search

eguil / Density_bining

Density bining code
2 stars 5 forks source link

Bug in integrated zonal mean #71

Closed eguil closed 5 years ago

eguil commented 5 years ago

These lines are wrong:

# Compute volume of isopycnals
            volBinz     = thickBinz  * areazt
            volBinza    = thickBinza * areazta
            volBinzp    = thickBinzp * areaztp
            volBinzi    = thickBinzi * areazti

Because the bathymetry in density coordinates cannot be taken into account (i.e. areazt does not depend on density). We then end up with a total volume of the ocean which is 1.7 e18 instead of 1.33 e18 m3. Solution to implement: areazt should be density dependent and computed for each time step. Also add under debug key the integral of isonvol to ensure 1.33 e18 is reached. Note that intensive quantities (T,S) and thickness are ok.

eguil commented 5 years ago

Solved with:

            # Compute volume of isopycnals
            # Create areai array with right dimensions to avoid loop
            areaitsig = npy.tile(npy.ma.reshape(areai,Nii*Nji), (N_s+1,1))
            areaitsig = npy.tile(npy.ma.reshape(areaitsig,(N_s+1)*Nii*Nji), (nyrtc,1))
            areaitsig = npy.ma.reshape(areaitsig,[nyrtc,N_s+1,Nji,Nii])
            # Create volume via zonal integral of thickness * area
            volBinz  = npy.ma.sum(thickBini *(1- thickBini.mask)*areaitsig, axis=3)
            volBinza = npy.ma.sum(thickBinia*(1-thickBinia.mask)*areaitsig, axis=3)
            volBinzp = npy.ma.sum(thickBinip*(1-thickBinip.mask)*areaitsig, axis=3)
            volBinzi = npy.ma.sum(thickBinii*(1-thickBinii.mask)*areaitsig, axis=3)

            voltoti = npy.ma.sum(volBinz)
            print '  Total volume in rho coordinates target grid (ref = 1.33 e+18)   : ', voltoti

Still a diff in total volume (1.28 instead of 1.37) using IPSL bug this may be due to the ESMF bug (see #54)

durack1 commented 5 years ago

@eguil I assume this is now fixed by https://github.com/eguil/Density_bining/commit/64987b47e976148bf94208404282b9680b14d2a0? I ask as I am now queuing the obs for a rerun against 64987b47e976148bf94208404282b9680b14d2a0

durack1 commented 5 years ago

@eguil @ysilvy it seems there is a new input required:

    densityBin(thetao,so,areacello,outfileDensity,debug=True,timeint='all')
TypeError: densityBin() takes at least 5 arguments (6 given)

For the obs files, where do I find the:

...
    - fileT(time,lev,lat,lon)   - 4D potential temperature array
    - fileS(time,lev,lat,lon)   - 4D salinity array

    - fileV(time,lev,lat,lon)   - 4D meridional velocity array <<-- THIS GUY?

    - fileFx(lat,lon)           - 2D array containing the cell area values
    - outFile(str)              - output file with full path specified.
...

I don't believe this will work for obs...

eguil commented 5 years ago

see #74