Changeset 1795 in lmdz_wrf

Timestamp:

Mar 12, 2018, 3:37:19 PM (7 years ago)

Author:

lfita

Message:

Adding:

`psl_ecmwf': sea-level pressure computation following ECMWF

Location:

trunk/tools

Files:

• diag_tools.py (modified) (1 diff)
• diagnostics.inf (modified) (1 diff)
• diagnostics.py (modified) (2 diffs)
• module_ForDiagnostics.f90 (modified) (2 diffs)
• module_ForDiagnosticsVars.f90 (modified) (2 diffs)

trunk/tools/diag_tools.py

@@ -833 +833 @@
 
     return mslpv, mslpdims, mslpvdims
+
+def Forcompute_psl_ecmwf(ps, hgt, ta1, pa2, unpa1, dimns, dimvns):
+    """ Function to compute the sea-level pressure following Mats Hamrud and Philippe Courtier [Pa]
+    Forcompute_psl_ptarget(ps, hgt, ta1, pa2, unpa1, dimns, dimvns)
+      [ps]= surface pressure values (assuming [[t],y,x]) [Pa]
+      [hgt]= opography (assuming [y,x]) [m]
+      [ta1]= air-temperature values at first half-mass level (assuming [[t],y,x]) [K]
+      [pa2]= pressure values at second full-mass levels (assuming [[t],y,x]) [Pa]
+      [unpa1]= pressure values at first half-mass levels (assuming [[t],y,x]) [Pa]
+      [dimns]= list of the name of the dimensions of [pa]
+      [dimvns]= list of the name of the variables with the values of the 
+        dimensions of [pa]
+    """
+    fname = 'Forcompute_psl_ecmwf'
+
+    vardims = dimns[:]
+    varvdims = dimvns[:]
+
+    if len(pa2.shape) == 3:
+        psl = np.zeros((pa2.shape[0],pa2.shape[1],pa2.shape[2]), dtype=np.float)
+        dx = pa2.shape[2]
+        dy = pa2.shape[1]
+        dt = pa2.shape[0]
+        pslt= fdin.module_fordiagnostics.compute_psl_ecmwf( ps=ps[:].transpose(),    \
+          hgt=hgt[:].transpose(), t=ta1[:].transpose(), press=pa2[:].transpose(),    \
+          unpress=unpa1[:].transpose(), d1=dx, d2=dy, d4=dt)
+        psl = pslt.transpose()
+    else:
+        print errormsg
+        print '  ' + fname + ': rank', len(pa2.shape), 'not ready !!'
+        print '  it only computes 3D [t,y,x] rank values'
+        quit(-1)
+
+    return psl, vardims, varvdims
 
 def Forcompute_psl_ptarget(pa, ps, ta, hgt, qv, target_pressure, dimns, dimvns):

trunk/tools/diagnostics.inf

@@ -24 +24 @@
 prw, WRFprw, WRFdens@QVAPOR
 psl, WRFmslp, WRFp@PSFC@HGT@WRFt@QVAPOR
+psl, WRFpsl_ecmwf, PSFC@HGT@WRFt@WRFp@ZNU@ZNW
 psl, WRFpsl_ptarget, WRFp@PSFC@WRFt@HGT@QVAPOR
 rvor, rvor, U@V

trunk/tools/diagnostics.py

@@ -79 +79 @@
   'OMEGAw', 'RAINTOT',                                                               \
   'rvors', 'td', 'turbulence', 'WRFcape_afwa', 'WRFclivi', 'WRFclwvl', 'WRFgeop',    \
-  'WRFmrso', 'WRFp',                                                                 \
+  'WRFmrso', 'WRFp', 'WRFpsl_ecmwf',                                                 \
   'WRFpsl_ptarget', 'WRFrvors', 'WRFslw', 'ws', 'wds', 'wss', 'WRFheight',           \
   'WRFheightrel', 'WRFua', 'WRFva', 'WRFzwind', 'WRFzwind_log', 'WRFzwindMO']
@@ -871 +871 @@
         ncvar.insert_variable(ncobj, 'zg', diagout, dnames, diagoutvd, newnc)
 
+# WRFmslp_ptarget sea-level pressure following ECMWF method as PSFC, HGT, WRFt, WRFp, ZNU, ZNW
+    elif diagn == 'WRFpsl_ecmwf':
+        var0 = ncobj.variables[depvars[0]][:]
+        var1 = ncobj.variables[depvars[1]][0,:,:]
+        var2 = WRFt[:,0,:,:]
+        var4 = WRFp[:,0,:,:]
+        var5 = ncobj.variables[depvars[4]][0,:]
+        var6 = ncobj.variables[depvars[5]][0,:]
+
+        # This is quite too appriximate!! passing pressure at half-levels to 2nd full 
+        #   level, using eta values at full (ZNW) and half (ZNU) mass levels
+        var3 = WRFp[:,0,:,:] + (var6[1] - var5[0])*(WRFp[:,1,:,:] - WRFp[:,0,:,:])/  \
+          (var5[1]-var5[0])
+
+        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
+        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
+
+        diagout = np.zeros(var0.shape, dtype=np.float)
+        diagout, diagoutd, diagoutvd = diag.Forcompute_psl_ecmwf(var0, var1, var2,   \
+          var3, var4, dnamesvar, dvnamesvar)
+
+        # Removing the nonChecking variable-dimensions from the initial list
+        varsadd = []
+        for nonvd in NONchkvardims:
+            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
+            varsadd.append(nonvd)
+
+        ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)
+
 # WRFmslp_ptarget sea-level pressure following ptarget method as WRFp, PSFC, WRFt, HGT, QVAPOR
     elif diagn == 'WRFpsl_ptarget':

trunk/tools/module_ForDiagnostics.f90

@@ -24 +24 @@
 ! compute_clt3D3: Computation of total cloudiness from a 3D CLDFRA being 3rd dimension the z-dim
 ! compute_clt: Computation of total cloudiness
+! compute_psl_ecmwf: Compute sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier [Pa]
 ! compute_massvertint1D: Subroutine to vertically integrate a 1D variable in eta vertical coordinates
 ! compute_vertint1D: Subroutine to vertically integrate a 1D variable in any vertical coordinates
@@ -602 +603 @@
   END SUBROUTINE compute_cape_afwa4D
 
+  SUBROUTINE compute_psl_ecmwf(ps, hgt, T, press, unpress, psl, d1, d2, d4)
+! Subroutine to compute sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier [Pa]
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: d1, d2, d4
+    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in)           :: ps, T, press, unpress
+    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: hgt
+    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out)          :: psl
+ 
+! Local
+    INTEGER                                              :: i, j, it
+
+!!!!!!! Variables
+! ps: surface pressure [Pa]
+! hgt: terrain height [m]
+! T: temperature at first half-mass level [K]
+! press: pressure at first full levels [Pa]
+! unpress: pressure at first mass (half) levels [Pa]
+! psl: sea-level pressure [Pa]
+
+    fname = 'compute_psl_ecmwf'
+
+    DO i=1, d1
+      DO j=1, d2
+        DO it=1, d4
+          CALL var_psl_ecmwf(ps(i,j,it), hgt(i,j), T(i,j,it), unpress(i,j,it), press(i,j,it),         &
+            psl(i,j,it))
+        END DO
+      END DO
+    END DO
+
+    RETURN
+
+  END SUBROUTINE compute_psl_ecmwf
+
   SUBROUTINE compute_zmla_generic4D(tpot, qratio, z, hgt, zmla3D, d1, d2, d3, d4)
 ! Subroutine to compute pbl-height following a generic method

trunk/tools/module_ForDiagnosticsVars.f90

@@ -28 +28 @@
 ! var_cllmh: low, medium, high-cloud [0,1]
 ! var_clt: total cloudiness [0,1]
+! var_psl_ecmwf: sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier [Pa]
 ! var_zmla_generic: Subroutine to compute pbl-height following a generic method
 ! var_zwind: Subroutine to extrapolate the wind at a given height following the 'power law' methodology
@@ -127 +128 @@
   END FUNCTION var_clt
 
+  SUBROUTINE var_psl_ecmwf(PRPRESS, hgt, PTB, PRESBH, PRESBF, psl)
+    ! Subroutine to compute sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier
+    !   method found in LMDZ in phylmd/pppmer.F90 in combination with phylmd/ctsar.F90
+
+!        IMPLICIT ARGUMENTS :  CONSTANTS FROM YOMCST,YOMGEM,YOMSTA.
+!        --------------------
+
+    IMPLICIT NONE
+
+    REAL, INTENT(in)                                     :: PRPRESS, hgt, PTB, PRESBH, PRESBF 
+    REAL, INTENT(out)                                    :: psl
+
+! Local
+    REAL                                                 :: ghgt, PTSTAR, PT0, ZTSTAR
+    REAL                                                 :: ZALPHA, POROG
+    REAL                                                 :: ZDTDZSG, ZOROG, ZT0, ZTX, ZTY, ZX, ZY, ZY2
+    REAL, PARAMETER                                      :: RDTDZ1 = -gammav
+
+!!!!!!! Variables
+! PRPRESS: Surface pressure [Pa]
+! hgt: Terrain height [m]
+! PTB: Temperature first half-level [K]
+! PRESBH: Pressure first half-level [Pa]
+! PRESBF: Pressure second full-level [Pa]
+! psl: sea-level pressure
+
+    fname = 'var_psl_ecmwf'
+    
+    ! Height by gravity
+    POROG = hgt*g
+
+    !* 1. COMPUTES SURFACE TEMPERATURE
+    !*   THEN STANDARD SURFACE TEMPERATURE.
+
+    ZDTDZSG=-RDTDZ1/g
+    ZALPHA=ZDTDZSG*r_d
+
+    PTSTAR=PTB*(1.0+ZALPHA*(PRESBH/PRESBF-1.0))
+    PT0=PTSTAR+ZDTDZSG*POROG
+
+    !* 2.    POST-PROCESS MSL PRESSURE.
+    !  --------------------------
+
+    !* 2.1   COMPUTATION OF MODIFIED ALPHA AND TSTAR.
+
+    ZTX=290.5
+    ZTY=255.0
+
+    IF (PTSTAR < ZTY) THEN
+      ZTSTAR=0.5*(ZTY+PTSTAR)
+    ELSEIF (PTSTAR < ZTX) THEN
+      ZTSTAR=PTSTAR
+    ELSE
+      ZTSTAR=0.5*(ZTX+PTSTAR)
+    ENDIF
+
+    ZT0=ZTSTAR+ZDTDZSG*POROG
+    IF (ZTX > ZTSTAR .AND. ZT0 > ZTX) THEN
+      ZT0=ZTX
+    ELSEIF (ZTX <= ZTSTAR .AND. ZT0 > ZTSTAR) THEN
+      ZT0=ZTSTAR
+    ELSE
+      ZT0=PT0
+    ENDIF
+
+    ZOROG=SIGN(MAX(1.0,ABS(POROG)),POROG)
+    ZALPHA=r_d*(ZT0-ZTSTAR)/ZOROG
+
+    !* 2.2   COMPUTATION OF MSL PRESSURE.
+
+    IF (ABS(POROG) >= 0.001) THEN
+      ZX=POROG/(r_d*ZTSTAR)
+      ZY=ZALPHA*ZX
+      ZY2=ZY*ZY
+
+      psl=PRPRESS*EXP(ZX*(1.0-0.5*ZY+1.0/3.*ZY2))
+    ELSE
+      psl=PRPRESS
+    ENDIF
+
+    RETURN
+
+  END SUBROUTINE var_psl_ecmwf
 
   SUBROUTINE compute_psl_ptarget4d2(press, ps, hgt, ta, qv, ptarget, psl, d1, d2, d3, d4)