Changeset 2675 in lmdz_wrf for trunk

Timestamp:

Jul 12, 2019, 10:31:58 PM (5 years ago)

Author:

lfita

Message:

Adding:

• 'frontogenesis': Calculation of the frontogenesis diagnostics

Location:

trunk/tools

Files:

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

trunk/tools/diag_tools.py

@@ -1600 +1600 @@
 
     if len(theta.shape) == 4:
-        dx = tas.shape[3]
-        dy = tas.shape[2]
-        dz = tas.shape[1]
-        dt = tas.shape[0]
+        dx = theta.shape[3]
+        dy = theta.shape[2]
+        dz = theta.shape[1]
+        dt = theta.shape[0]
 
         xdiab = np.zeros((dt,dz,dy,dx), dtype=np.float)
@@ -1623 +1623 @@
         dxst = dxs.transpose()
         dyst = dys.transpose()
-        dzst = dzs.transpose()
+
+        print '  ' + fname + ' Lluis shapes dzs', dzs.shape, 'dz:', dz
  
         pxdiab, pydiab, pzdiab, pxdef, pydef, pzdef, pxtilt, pytilt, pzdiv, pf =     \
           fdin.module_fordiagnostics.compute_frontogenesis(theta=thetat, ua=uat,     \
-          va=vat, wa=wat, press=presst, dsx=dxst, dsy=dyst, dsz=dzst, dst=dts, d1=dx,\
+          va=vat, wa=wat, press=presst, dsx=dxst, dsy=dyst, dsz=dzs, dst=dts, d1=dx, \
           d2=dy, d3=dz, d4=dt)
         xdiab = pxdiab.transpose()

trunk/tools/diagnostics.inf

@@ -18 +18 @@
 fog, fog_RUC, Q2@T2@PSFC
 front, front_R04, T2@U10@V10@WRFdx@WRFdy@params:1.:0.0001
+frontogenesis, frontogenesis, WRFt@UNua@UNva@UNwa@WRFp@WRFdx@WRFdy@WRFdz@WRFtime
 hur, LMDZrh, pres@t@r
 hur, WRFrh, WRFrh@T@P@PB

trunk/tools/diagnostics.py

@@ -106 +106 @@
   'WRFbils',  'WRFbnds',                                                             \
   'WRFclivi', 'WRFclwvi', 'WRFdens', 'WRFdx', 'WRFdxdy', 'WRFdxdywps', 'WRFdy',      \
-  'WRFgeop', 'WRFp', 'WRFtd',                                                        \
+  'WRFdz', 'WRFgeop', 'WRFp', 'WRFtd',                                               \
   'WRFpos', 'WRFprc', 'WRFprls', 'WRFrh', 'LMDZrh', 'LMDZrhs',                       \
   'WRFrhs', 'WRFrvors',                                                              \
@@ -173 +173 @@
 WRFdxdywps_compute = False
 LONLATdxdy_compute = False
+UNua_compute = False
+UNva_compute = False
+UNwa_compute = False
 
 # File creation
@@ -486 +489 @@
     print '  ' + main + ': Retrieving dz: real distance between grid points ' +      \
       'from WRF as dz=PHB(k+1)+PH(k+1)-(PHB(k)+PH(k))'
+    PH = ncobj.variables['PH'][:]
+    PHB = ncobj.variables['PHB'][:]
+
     dimv = ncobj.variables['PH'].shape
 
@@ -876 +882 @@
     elif diagn == 'frontogenesis':
             
-        var0 = ncobj.variables[depvars[0]][:]
+        if depvars[0] == 'WRFt': var0 = WRFt
+        else: var0 = ncobj.variables[depvars[0]][:]
+        dx = var0.shape[3]
+        dy = var0.shape[2]
+        dz = var0.shape[1]
+        dt = var0.shape[0]
         var1 = ncobj.variables[depvars[1]][:]
         var2 = ncobj.variables[depvars[2]][:]
-        if depvars[3] == 'WRFpress': var3 = WRFpress
-        else: var3 = ncobj.variables[depvars[3]]
-        if depvars[4] == 'WRFdx': var4 = WRFdx
+        var3 = ncobj.variables[depvars[3]][:]
+        if depvars[4] == 'WRFp': var4 = WRFp
         else: var4 = ncobj.variables[depvars[4]]
-        if depvars[5] == 'WRFdy': var5 = WRFdy
+        if depvars[5] == 'WRFdx': var5 = WRFdx
         else: var5 = ncobj.variables[depvars[5]]
-        if depvars[6] == 'WRFdz': var6 = WRFdz[0,:,0,0]
+        if depvars[6] == 'WRFdy': var6 = WRFdy
         else: var6 = ncobj.variables[depvars[6]]
-        if depvars[7] == 'WRFtime': var7 = WRFtime
+        if depvars[7] == 'WRFdz': var7 = WRFdz[0,0:dz,0,0]
         else: var7 = ncobj.variables[depvars[7]]
+        if depvars[8] == 'WRFtime': var08 = WRFtime
+        else: var08 = ncobj.variables[depvars[8]]
+
+        # Assuming monotonic time-axis...
+        var8 = var08[1] - var08[0]
 
         diag1, diag2, diag3, diag4, diag5, diag6, diag7, diag8, diag9, diag10,       \
           diagoutd, diagoutvd = diag.Forcompute_frontogenesis(var0, var1, var2, var3,\
-          var4, var5, var6, var7, dnames,dvnames)
+          var4, var5, var6, var7, var8, dnames, dvnames)
 
         ncvar.insert_variable(ncobj, 'diabh', diag1, diagoutd, diagoutvd, newnc,     \

trunk/tools/module_ForDiagnostics.f90

@@ -1596 +1596 @@
     INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
     REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: theta, ua, va, wa, press
-    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: dsx, dsy, dsz, dst
+    REAL(r_k), INTENT(in)                                :: dst
+    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: dsx, dsy
+    REAL(r_k), DIMENSION(d3), INTENT(in)                 :: dsz
     REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(out)       :: xdiab, ydiab, zdiab, xdef, ydef, zdef,    &
-      xtilt, ytilt, zdiv, f
+      xtilt, ytilt, zdiv
     REAL(r_k), DIMENSION(d1,d2,d3,d4,3), INTENT(out)     :: f
  
@@ -1623 +1625 @@
 
     DO it=1, d4
-      CALL var_Frontogenesis(dx, dy, dz, theta(:,:,:,it), ua(:,:,:,it), va(:,:,:,it), wa(:,:,:,it),   &
-        press(:,:,:,it), ddx, ddy, ddz, ddt, xdiab(:,:,:,it), ydiab(:,:,:,it), zdiab(:,:,:,it),       &
+      CALL var_Frontogenesis(d1, d2, d3, theta(:,:,:,it), ua(:,:,:,it), va(:,:,:,it), wa(:,:,:,it),   &
+        press(:,:,:,it), dsx, dsy, dsz, dst, xdiab(:,:,:,it), ydiab(:,:,:,it), zdiab(:,:,:,it),       &
         xdef(:,:,:,it), ydef(:,:,:,it), zdef(:,:,:,it), xtilt(:,:,:,it), ytilt(:,:,:,it),             &
         zdiv(:,:,:,it), f(:,:,:,it,:))
@@ -1631 +1633 @@
     RETURN
 
-  END SUBROUTINE compute_front_R04d3
+  END SUBROUTINE compute_frontogenesis
 
 END MODULE module_ForDiagnostics

trunk/tools/module_ForDiagnosticsVars.f90

@@ -2059 +2059 @@
     REAL(r_k), DIMENSION(dx,dy), INTENT(in)              :: ddx,ddy
     REAL(r_k), DIMENSION(dz), INTENT(in)                 :: ddz
-    REAL(r_k), DIMENSION(dx,dy,dz,), INTENT(in)          :: theta, ua, va, wa, press
-    REAL(r_k), DIMENSION(dx,dy,dz,), INTENT(out)         :: xdiab, ydiab, zdiab
-    REAL(r_k), DIMENSION(dx,dy,dz,), INTENT(out)         :: xdef, ydef, zdef
-    REAL(r_k), DIMENSION(dx,dy,dz,), INTENT(out)         :: xtilt, ytilt, zdiv
+    REAL(r_k), DIMENSION(dx,dy,dz), INTENT(in)           :: theta, ua, va, wa, press
+    REAL(r_k), DIMENSION(dx,dy,dz), INTENT(out)          :: xdiab, ydiab, zdiab
+    REAL(r_k), DIMENSION(dx,dy,dz), INTENT(out)          :: xdef, ydef, zdef
+    REAL(r_k), DIMENSION(dx,dy,dz), INTENT(out)          :: xtilt, ytilt, zdiv
     REAL(r_k), DIMENSION(dx,dy,dz,3), INTENT(out)        :: f
 
@@ -2110 +2110 @@
     zdiv(:,:,:) = thetatilt(:,:,:,3)
 
-    f(:,:,:,1) = 1./modthteagrad*thetagrad(:,:,:,1)*(xdiab(:,:,:) + xdef(:,:,:) +          &
-      xtilt(:,:,:))
+    f(:,:,:,1) = 1./modthetagrad*thetagrad(:,:,:,1)*(xdiab(:,:,:) + xdef(:,:,:) + xtilt(:,:,:))
     f(:,:,:,2) = thetagrad(:,:,:,2)*(ydiab(:,:,:) + ydef(:,:,:) + ytilt(:,:,:))
     f(:,:,:,3) = thetagrad(:,:,:,3)*(zdiab(:,:,:) + zdef(:,:,:) + zdiv(:,:,:))

trunk/tools/module_scientific.f90

@@ -7923 +7923 @@
   END FUNCTION vecmodule3D
 
-  SUBROUTINE FUNCTION matmodule3D(d1, d2, d3, mat, matmodule)
+  SUBROUTINE matmodule3D(d1, d2, d3, mat, matmodule)
   ! Subroutine to compute the module of a 3D matrix with 3 components
 
     IMPLICIT NONE
 
+    INTEGER, INTENT(in)                                  :: d1,d2,d3
     REAL(r_k), DIMENSION(d1,d2,d3,3), INTENT(in)         :: mat
     REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out)          :: matmodule
@@ -7933 +7934 @@
     fname = 'matmodule3D'
 
-    matmodule3D = SQRT(vec(:,:,:,1)*vec(:,:,:,1) + vec(:,:,:,2)*vec(:,:,:,2) +                        &
-      vec(:,:,:,3)*vec(:,:,:,3))
+    matmodule = SQRT(mat(:,:,:,1)*mat(:,:,:,1) + mat(:,:,:,2)*mat(:,:,:,2) + mat(:,:,:,3)*mat(:,:,:,3))
 
   END SUBROUTINE matmodule3D
@@ -7960 +7960 @@
     def = fillval64
 
-    def(:,:,:,1) = - uagrad(:,:,:,1)*vargrad(:,:,:1) - vagrad(:,:,:,1)*vargrad(:,:,:2)
-    def(:,:,:,2) = - uagrad(:,:,:,2)*vargrad(:,:,:1) - vagrad(:,:,:,2)*vargrad(:,:,:2)
-    def(:,:,:,3) = - uagrad(:,:,:,3)*vargrad(:,:,:1) - vagrad(:,:,:,3)*vargrad(:,:,:2)
+    def(:,:,:,1) = -uagrad(:,:,:,1)*vargrad(:,:,:,1) -vagrad(:,:,:,1)*vargrad(:,:,:,2)
+    def(:,:,:,2) = -uagrad(:,:,:,2)*vargrad(:,:,:,1) -vagrad(:,:,:,2)*vargrad(:,:,:,2)
+    def(:,:,:,3) = -uagrad(:,:,:,3)*vargrad(:,:,:,1) -vagrad(:,:,:,3)*vargrad(:,:,:,2)
 
   END SUBROUTINE deformation3D
@@ -7988 +7988 @@
     tilt = fillval64
 
-    tilt(:,:,:,1) = - uagrad(:,:,:,1)*vargrad(:,:,:3)
-    tilt(:,:,:,2) = - vagrad(:,:,:,2)*vargrad(:,:,:3)
-    tilt(:,:,:,3) = - wagrad(:,:,:,3)*vargrad(:,:,:3)
+    tilt(:,:,:,1) = -wagrad(:,:,:,1)*vargrad(:,:,:,3)
+    tilt(:,:,:,2) = -wagrad(:,:,:,2)*vargrad(:,:,:,3)
+    tilt(:,:,:,3) = -wagrad(:,:,:,3)*vargrad(:,:,:,3)
 
   END SUBROUTINE tilting3D
@@ -8000 +8000 @@
 
     INTEGER, INTENT(in)                                  :: dx, dy, dz, dt
-    REAL(r_k), DIMENSION(dx,dy,dz,dt), INTENT(in)        :: vargrad
+    REAL(r_k), INTENT(in)                                :: ddt
+    REAL(r_k), DIMENSION(dx,dy,dz,dt), INTENT(in)        :: var
     CHARACTER(LEN=*), INTENT(in)                         :: dtkind
     REAL(r_k), DIMENSION(dx,dy,dz,dt), INTENT(out)       :: vardt
@@ -8039 +8040 @@
     END IF    
 
-  END SUBROUTINE tilting3D
+  END SUBROUTINE deltat3D
 
 END MODULE module_scientific