Remove Fortran 2018 warnings

README.md

@@ -1,8 +1,6 @@
 # GFDL_atmos_cubed_sphere
 
-The source contained herein reflects the 201912 release of the Finite Volume Cubed-Sphere Dynamical Core (FV3) from GFDL
-
-The GFDL Microphysics is also available via this repository.
+The source contained herein reflects the 201912 release of the Finite Volume Cubed-Sphere Dynamical Core (FV3) for use in the current GFDL models (AM4/CM4/ESM4/SPEAR).
 
 # Where to find information
 
@@ -12,7 +10,6 @@ for more information.
 # Proper usage attribution
 
 Cite Putman and Lin (2007) and Harris and Lin (2013) when describing a model using the FV3 dynamical core.
-Cite Chen et al (2013) and Zhou et al (2019) when using the GFDL Microphysics.
 
 # What files are what
 
@@ -27,6 +24,13 @@ The top level directory structure groups source code and input files as follow:
 | ```tools/```         | contains source code of tools used within the core |
 | ```GFDL_tools/```    | contains source code of tools specific to GFDL models |
 
+# Dependencies
+
+The source code in this repository requires other NOAA-GFDL projects in
+order to compile.  As stated above, this branch is for use with the current
+NOAA-GFDL/AM4 and NOAA-GFDL/CM4 projects and is included as a submodule within
+those projects.
+
 # Disclaimer
 
 The United States Department of Commerce (DOC) GitHub project code is provided

RELEASE.md

@@ -1,31 +1,31 @@
-# RELEASE NOTES for FV3: Summary
+# RELEASE NOTES for GFDL  FV3: Summary
 
-FV3-201912-public --- 10 January 2020
-Lucas Harris, GFDL <[email protected]>
+2020.02 --- 22 April 2020
 
-This version has been tested against the current SHiELD (formerly fvGFS) physics
-and with FMS release candidate 2020.02 from https://github.com/NOAA-GFDL/FMS
+This version has been tested within current GFDL Models (AM4+, CM4+, ESM4+, SPEAR, etc.) and requires the 2020.02 release of the [FMS infrastructure](https://github.com/NOAA-GFDL/FMS).
 
-Includes all of the features of the GFDL Release to EMC, as well as:
+Includes all of the features from the [201912 Public Release](https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere/releases/tag/201912_public_release) which include:
 
-- Updated 2017 GFDL Microphysics (from S-J Lin and L Zhou included in GFSv15)
-- Updates for GFSv15 ICs (from T Black/J Abeles, EMC)
-- Updates to support new nesting capabilities in FMS (from Z Liang)
+- Updates to support new nesting capabilities in FMS
 - Re-written grid nesting code for efficiency and parallelization
 - Re-organized fv_eta for improved vertical level selection
-- 2018 Stand-alone regional capabilities (from T Black/J Abeles, EMC)
+- 2018 Stand-alone regional capabilities (from EMC)
 - Refactored model front-end (fv_control, fv_restart)
 - Support for point soundings
+- full non-hydrostatic capability is now included as a runtime option
 - And other updates
 
-# Interface changes
+# Directory structure changes
 
-drivers: renamed 'fvGFS' directory to SHiELD
+***drivers/***  (important for those moving from the GFDL internal project)
+  - renamed ***fvGFS*** to ***SHiELD***
+  - renamed ***coupled*** to ***GFDL***
 
-atmosphere.F90: 'mytile' is renamed 'mygrid'
+***model_nh_null/***
+  - has been removed
 
-The non-functional gfdl_cloud_microphys.F90 has been removed and replaced with the 2017 public release given to EMC. Also added a proper initialization routine, that includes the use of INTERNAL_FILE_NML and thereby requires the input_nml_file argument. If you do not define the compiler flag INTERNAL_FILE_NML then you can delete this argument.
+Update your build system as appropriate
 
-The namelist nggps_diag_nml has been eliminated. 'fdiag' is no longer handled by the dynamical core, and should be handled by the physics driver.
+# Documentation
 
 For a complete technical description see the NOAA Technical Memorandum OAR GFDL: https://repository.library.noaa.gov/view/noaa/23432

driver/GFDL/atmosphere.F90

@@ -1264,7 +1264,7 @@ subroutine fv_compute_p_z (npz, phis, pe, peln, delp, delz, pt, q_sph, &
     endif
     if (do_uni_zfull) then
        do k=1,npz
-       	  z_full(:,:,k)=0.5*(z_half(:,:,k)+z_half(:,:,k+1))
+          z_full(:,:,k)=0.5*(z_half(:,:,k)+z_half(:,:,k+1))
        enddo
     endif
   end subroutine fv_compute_p_z

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, &

model/fv_sg.F90

@@ -26,7 +26,8 @@ module fv_sg_mod
   use constants_mod,      only: rdgas, rvgas, cp_air, cp_vapor, hlv, hlf, kappa, grav
   use tracer_manager_mod, only: get_tracer_index
   use field_manager_mod,  only: MODEL_ATMOS
-  use gfdl_cloud_microphys_mod, only: wqs1, wqs2, wqsat2_moist
+!  use gfdl_cloud_microphys_mod, only: wqs1, wqs2, wqsat2_moist
+  use lin_cld_microphys_mod, only: wqs2, wqsat2_moist
   use fv_mp_mod,          only: mp_reduce_min, is_master
 
 implicit none

model/sw_core.F90

@@ -979,8 +979,10 @@ subroutine d_sw(delpc, delp,  ptc,   pt, u,  v, w, uc,vc, &
 !    endif
         call fv_tp_2d(pt, crx_adv,cry_adv, npx, npy, hord_tm, gx, gy,  &
                       xfx_adv,yfx_adv, gridstruct, bd, ra_x, ra_y,     &
-                      mfx=fx, mfy=fy, mass=delp, nord=nord_v, damp_c=damp_v)
-!                     mfx=fx, mfy=fy, mass=delp, nord=nord_t, damp_c=damp_t)
+! This is a difference from the publicly released version in github master
+! and is provided to ensure reproducibility with AM4/CM4 models from GFDL
+!                     mfx=fx, mfy=fy, mass=delp, nord=nord_v, damp_c=damp_v)
+                      mfx=fx, mfy=fy, mass=delp, nord=nord_t, damp_c=damp_t)
 #endif
 
      if ( inline_q ) then

tools/fv_diagnostics.F90

@@ -46,7 +46,8 @@ module fv_diagnostics_mod
  use sat_vapor_pres_mod, only: compute_qs, lookup_es
 
  use fv_arrays_mod, only: max_step
- use gfdl_cloud_microphys_mod, only: wqs1, qsmith_init
+! use gfdl_cloud_microphys_mod, only: wqs1, qsmith_init
+ use lin_cld_microphys_mod, only: wqs1, qsmith_init
 
  use column_diagnostics_mod, only:  column_diagnostics_init, &
                                     initialize_diagnostic_columns, &
@@ -707,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)
@@ -1576,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)
 
@@ -1731,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
 
@@ -4236,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
@@ -4919,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
@@ -4930,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
@@ -4941,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.
@@ -4952,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
 
@@ -4965,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
@@ -4985,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

tools/fv_mp_mod.F90

@@ -28,7 +28,7 @@ module fv_mp_mod
 ! !USES:
       use fms_mod,         only : fms_init, fms_end
       use mpp_mod,         only : FATAL, MPP_DEBUG, NOTE, MPP_CLOCK_SYNC,MPP_CLOCK_DETAILED, WARNING
-      use mpp_mod,         only : mpp_pe, mpp_npes, mpp_node, mpp_root_pe, mpp_error, mpp_set_warn_level
+      use mpp_mod,         only : mpp_pe, mpp_npes, mpp_root_pe, mpp_error, mpp_set_warn_level
       use mpp_mod,         only : mpp_declare_pelist, mpp_set_current_pelist, mpp_sync
       use mpp_mod,         only : mpp_clock_begin, mpp_clock_end, mpp_clock_id
       use mpp_mod,         only : mpp_chksum, stdout, stderr, mpp_broadcast