Merge pull request #47 from wfcooke/user/wfc/climate_nudging

User/wfc/climate nudging

model/fv_arrays.F90

@@ -60,6 +60,9 @@ module fv_arrays_mod
            id_dbz, id_maxdbz, id_basedbz, id_dbz4km, id_dbztop, id_dbz_m10C, &
            id_ctz, id_w1km, id_wmaxup, id_wmaxdn, id_cape, id_cin
 
+! Selected theta-level fields from 3D variables:
+ integer :: id_pv350K, id_pv550K
+
 ! Selected p-level fields from 3D variables:
  integer :: id_vort200, id_vort500, id_w500, id_w700
  integer :: id_vort850, id_w850, id_x850, id_srh25, &

tools/fv_diagnostics.F90

@@ -708,6 +708,10 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref)
           !--------------------
           idiag%id_pv = register_diag_field ( trim(field), 'pv', axes(1:3), Time,       &
                'potential vorticity', '1/s', missing_value=missing_value )
+          idiag%id_pv350K = register_diag_field ( trim(field), 'pv350K', axes(1:2), Time,       &
+                           '350-K potential vorticity; needs x350 scaling', '(K m**2) / (kg s)', missing_value=missing_value)
+          idiag%id_pv550K = register_diag_field ( trim(field), 'pv550K', axes(1:2), Time,       &
+                           '550-K potential vorticity; needs x550 scaling', '(K m**2) / (kg s)', missing_value=missing_value)
 
           ! -------------------
           ! Vertical flux correlation terms (good for averages)
@@ -1577,8 +1581,8 @@ subroutine fv_diag(Atm, zvir, Time, print_freq)
        endif
 
        if ( idiag%id_vort200>0 .or. idiag%id_vort500>0 .or. idiag%id_vort850>0 .or. idiag%id_vorts>0   &
-            .or. idiag%id_vort>0 .or. idiag%id_pv>0 .or. idiag%id_rh>0 .or. idiag%id_x850>0 .or.  &
-            idiag%id_uh03>0 .or. idiag%id_uh25>0) then
+            .or. idiag%id_vort>0 .or. idiag%id_pv>0 .or. idiag%id_pv350K>0 .or. idiag%id_pv550K>0 &
+            .or. idiag%id_rh>0 .or. idiag%id_x850>0 .or. idiag%id_uh03>0 .or. idiag%id_uh25>0) then
           call get_vorticity(isc, iec, jsc, jec, isd, ied, jsd, jed, npz, Atm(n)%u, Atm(n)%v, wk, &
                Atm(n)%gridstruct%dx, Atm(n)%gridstruct%dy, Atm(n)%gridstruct%rarea)
 
@@ -1732,11 +1736,36 @@ subroutine fv_diag(Atm, zvir, Time, print_freq)
           endif
 
 
-          if ( idiag%id_pv > 0 ) then
-! Note: this is expensive computation.
+          if ( idiag%id_pv > 0 .or. idiag%id_pv350K >0 .or. idiag%id_pv550K >0) then
+              if (allocated(a3)) deallocate(a3)
+              allocate ( a3(isc:iec,jsc:jec,npz+1) )
+            ! Modified pv_entropy to get potential temperature at layer interfaces (last variable)
+            ! The values are needed for interpolate_z
+            ! Note: this is expensive computation.
               call pv_entropy(isc, iec, jsc, jec, ngc, npz, wk,    &
-                              Atm(n)%gridstruct%f0, Atm(n)%pt, Atm(n)%pkz, Atm(n)%delp, grav)
-              used = send_data ( idiag%id_pv, wk, Time )
+                              Atm(n)%gridstruct%f0, Atm(n)%pt, Atm(n)%pkz, Atm(n)%delp, grav,a3)
+	      if ( idiag%id_pv > 0) then
+                used = send_data ( idiag%id_pv, wk, Time )
+	      endif
+	      if( idiag%id_pv350K > 0 .or. idiag%id_pv550K >0 ) then
+                !"pot temp" from pv_entropy is only semi-finished; needs p0^kappa (pk0)
+                do k=1,npz+1
+                   do j=jsc,jec
+                      do i=isc,iec
+                         a3(i,j,k) = a3(i,j,k) * pk0  
+                      enddo
+           	   enddo
+                enddo
+                !wk as input, which stores pv from pv_entropy;
+                !use z interpolation, both potential temp and z increase monotonically with height
+                !interpolate_vertical will apply log operation to 350, and also assumes a different vertical order
+		call interpolate_z(isc, iec, jsc, jec, npz, 350., a3, wk, a2)
+		used = send_data( idiag%id_pv350K, a2, Time)
+                !interpolate_vertical will apply log operation to 350, and also assumes a different vertical order
+		call interpolate_z(isc, iec, jsc, jec, npz, 550., a3, wk, a2)
+		used = send_data( idiag%id_pv550K, a2, Time)
+	      endif
+              deallocate ( a3 )
               if (prt_minmax) call prt_maxmin('PV', wk, isc, iec, jsc, jec, 0, 1, 1.)
           endif
 
@@ -4237,6 +4266,7 @@ subroutine get_height_given_pressure(is, ie, js, je, km, wz, kd, id, log_p, peln
                  go to 1000
              endif
           enddo
+       a2(i,j,n) =missing_value
 1000   continue
     enddo
  enddo
@@ -4920,7 +4950,7 @@ end subroutine updraft_helicity
 
 
 
- subroutine pv_entropy(is, ie, js, je, ng, km, vort, f_d, pt, pkz, delp, grav)
+ subroutine pv_entropy(is, ie, js, je, ng, km, vort, f_d, pt, pkz, delp, grav, te)
 
 ! !INPUT PARAMETERS:
    integer, intent(in)::  is, ie, js, je, ng, km
@@ -4931,7 +4961,9 @@ subroutine pv_entropy(is, ie, js, je, ng, km, vort, f_d, pt, pkz, delp, grav)
    real, intent(in):: f_d(is-ng:ie+ng,js-ng:je+ng)
 
 ! vort is relative vorticity as input. Becomes PV on output
-      real, intent(inout):: vort(is:ie,js:je,km)
+   real, intent(inout):: vort(is:ie,js:je,km)
+! output potential temperature at the interface so it can be used for diagnostics
+   real, intent(out):: te(is:ie,js:je,km+1) 
 
 ! !DESCRIPTION:
 !        EPV = 1/r * (vort+f_d) * d(S)/dz; where S is a conservative scalar
@@ -4942,7 +4974,7 @@ subroutine pv_entropy(is, ie, js, je, ng, km, vort, f_d, pt, pkz, delp, grav)
 ! z-surface is not that different from the hybrid sigma-p coordinate.
 ! See page 39, Pedlosky 1979: Geophysical Fluid Dynamics
 !
-! The follwoing simplified form is strictly correct only if vort is computed on
+! The following simplified form is strictly correct only if vort is computed on
 ! constant z surfaces. In addition hydrostatic approximation is made.
 !        EPV = - GRAV * (vort+f_d) / del(p) * del(pt) / pt
 ! where del() is the vertical difference operator.
@@ -4953,7 +4985,7 @@ subroutine pv_entropy(is, ie, js, je, ng, km, vort, f_d, pt, pkz, delp, grav)
 !---------------------------------------------------------------------
 !BOC
       real w3d(is:ie,js:je,km)
-      real te(is:ie,js:je,km+1), t2(is:ie,km), delp2(is:ie,km)
+      real t2(is:ie,km), delp2(is:ie,km)
       real te2(is:ie,km+1)
       integer i, j, k
 
@@ -4966,8 +4998,8 @@ subroutine pv_entropy(is, ie, js, je, ng, km, vort, f_d, pt, pkz, delp, grav)
         enddo
 #else
 ! Compute PT at layer edges.
-!$OMP parallel do default(none) shared(is,ie,js,je,km,pt,pkz,w3d,delp,te2,te) &
-!$OMP                          private(t2, delp2)
+!$OMP parallel do default(none) shared(is,ie,js,je,km,pt,pkz,w3d,delp,te) &
+!$OMP                          private(t2, te2, delp2)
      do j=js,je
         do k=1,km
           do i=is,ie
@@ -4986,7 +5018,8 @@ subroutine pv_entropy(is, ie, js, je, ng, km, vort, f_d, pt, pkz, delp, grav)
         enddo
      enddo
 
-!$OMP parallel do default(none) shared(is,ie,js,je,km,vort,f_d,te,w3d,delp,grav)
+!$OMP parallel do default(none) shared(is,ie,js,je,km,vort,te,w3d,delp,grav) &
+!$OMP                           private(f_d)
      do k=1,km
         do j=js,je
           do i=is,ie