From 1c6c1d831a1a1e65a70ec0bf6bb5a411f543afa8 Mon Sep 17 00:00:00 2001 From: Dusan Jovic Date: Thu, 24 Mar 2022 14:45:19 +0000 Subject: [PATCH 01/21] Resolve argument mismatch errors when using gfortran Switch from 'use mpi' to 'use mpi_f08' --- physics/GFS_debug.F90 | 4 +- physics/aerinterp.F90 | 33 ++++--------- physics/cires_tauamf_data.F90 | 9 +--- physics/h2ointerp.f90 | 10 +--- physics/iccninterp.F90 | 47 ++++++++---------- physics/module_mp_thompson.F90 | 9 ++-- physics/mp_fer_hires.meta | 2 +- physics/mp_thompson.F90 | 5 +- physics/mp_thompson.meta | 4 +- physics/mp_thompson_post.F90 | 3 +- physics/mp_thompson_post.meta | 2 +- physics/ozinterp.f90 | 10 +--- physics/rrtmgp_lw_cloud_optics.F90 | 5 +- physics/rrtmgp_lw_cloud_optics.meta | 2 +- physics/rrtmgp_lw_gas_optics.F90 | 5 +- physics/rrtmgp_lw_gas_optics.meta | 2 +- physics/rrtmgp_sw_cloud_optics.F90 | 5 +- physics/rrtmgp_sw_cloud_optics.meta | 2 +- physics/rrtmgp_sw_gas_optics.F90 | 5 +- physics/rrtmgp_sw_gas_optics.meta | 2 +- physics/sfcsub.F | 74 ++++------------------------- 21 files changed, 74 insertions(+), 166 deletions(-) diff --git a/physics/GFS_debug.F90 b/physics/GFS_debug.F90 index f01f25cbc..0b74bf851 100644 --- a/physics/GFS_debug.F90 +++ b/physics/GFS_debug.F90 @@ -388,7 +388,7 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, nthreads, blkno, errmsg, errflg) #ifdef MPI - use mpi + use mpi_f08 #endif #ifdef _OPENMP use omp_lib @@ -1041,7 +1041,7 @@ subroutine GFS_interstitialtoscreen_run (Model, Statein, Stateout, Sfcprop, Coup nthreads, blkno, errmsg, errflg) #ifdef MPI - use mpi + use mpi_f08 #endif #ifdef _OPENMP use omp_lib diff --git a/physics/aerinterp.F90 b/physics/aerinterp.F90 index 30ff97dff..b4673a7bd 100644 --- a/physics/aerinterp.F90 +++ b/physics/aerinterp.F90 @@ -112,8 +112,6 @@ SUBROUTINE read_aerdataf ( me, master, iflip, idate, FHOUR, errmsg, errflg) integer IDAT(8),JDAT(8) real(kind=kind_phys) RINC(5), rjday integer jdow, jdoy, jday - real(4) rinc4(5) - integer w3kindreal,w3kindint integer, allocatable :: invardims(:) ! @@ -131,13 +129,7 @@ SUBROUTINE read_aerdataf ( me, master, iflip, idate, FHOUR, errmsg, errflg) IDAT(5) = IDATE(1) RINC = 0. RINC(2) = FHOUR - call w3kind(w3kindreal,w3kindint) - if(w3kindreal == 4) then - rinc4 = rinc - CALL W3MOVDAT(RINC4,IDAT,JDAT) - else - CALL W3MOVDAT(RINC,IDAT,JDAT) - endif + CALL W3MOVDAT(RINC,IDAT,JDAT) ! jdow = 0 jdoy = 0 @@ -246,8 +238,6 @@ SUBROUTINE aerinterpol( me,master,nthrds,npts,IDATE,FHOUR,iflip, jindx1,jindx2, real(kind=kind_phys) prsl(npts,lev), aerpres(npts,levsaer) real(kind=kind_phys) RINC(5), rjday integer jdow, jdoy, jday - real(4) rinc4(5) - integer w3kindreal,w3kindint ! IDAT = 0 @@ -257,13 +247,7 @@ SUBROUTINE aerinterpol( me,master,nthrds,npts,IDATE,FHOUR,iflip, jindx1,jindx2, IDAT(5) = IDATE(1) RINC = 0. RINC(2) = FHOUR - call w3kind(w3kindreal,w3kindint) - if(w3kindreal == 4) then - rinc4 = rinc - CALL W3MOVDAT(RINC4,IDAT,JDAT) - else - CALL W3MOVDAT(RINC,IDAT,JDAT) - endif + CALL W3MOVDAT(RINC,IDAT,JDAT) ! jdow = 0 jdoy = 0 @@ -393,6 +377,7 @@ subroutine read_netfaer(nf, iflip,nt) use aerclm_def use netcdf integer, intent(in) :: iflip, nf, nt + integer :: ncerr integer :: ncid, varid, i,j,k,ii,klev character :: fname*50, mn*2, vname*10 real(kind=kind_io4),allocatable,dimension(:,:,:) :: buff @@ -406,11 +391,11 @@ subroutine read_netfaer(nf, iflip,nt) write(mn,'(i2.2)') nf fname=trim("aeroclim.m"//mn//".nc") - call nf_open(fname , nf90_NOWRITE, ncid) + ncerr = nf90_open(fname , NF90_NOWRITE, ncid) ! ====> construct 3-d pressure array (Pa) - call nf_inq_varid(ncid, "DELP", varid) - call nf_get_var(ncid, varid, buff) + ncerr = nf90_inq_varid(ncid, "DELP", varid) + ncerr = nf90_get_var(ncid, varid, buff) do j = jamin, jamax do i = iamin, iamax @@ -439,8 +424,8 @@ subroutine read_netfaer(nf, iflip,nt) ! for GFS, iflip 0: toa to sfc; 1: sfc to toa DO ii = 1, ntrcaerm vname=trim(specname(ii)) - call nf_inq_varid(ncid, vname, varid) - call nf_get_var(ncid, varid, buffx) + ncerr = nf90_inq_varid(ncid, vname, varid) + ncerr = nf90_get_var(ncid, varid, buffx) do j = jamin, jamax do k = 1, levsaer @@ -462,7 +447,7 @@ subroutine read_netfaer(nf, iflip,nt) ENDDO ! ii-loop (ntracaerm) ! close the file - call nf_close(ncid) + ncerr = nf90_close(ncid) deallocate (buff, pres_tmp) deallocate (buffx) return diff --git a/physics/cires_tauamf_data.F90 b/physics/cires_tauamf_data.F90 index e0d43e74e..52bfb67bd 100644 --- a/physics/cires_tauamf_data.F90 +++ b/physics/cires_tauamf_data.F90 @@ -179,7 +179,6 @@ subroutine gfs_idate_calendar(idate, fhour, ddd, fddd) real(kind=kind_phys) :: rinc(5), rjday integer :: jdow, jdoy, jday real(4) :: rinc4(5) - integer :: w3kindreal, w3kindint integer :: iw3jdn integer :: jd1, jddd @@ -195,13 +194,7 @@ subroutine gfs_idate_calendar(idate, fhour, ddd, fddd) rinc(1:5) = 0. rinc(2) = fhour ! - call w3kind(w3kindreal,w3kindint) - if(w3kindreal==4) then - rinc4 = rinc - call w3movdat(rinc4, idat,jdat) - else - call w3movdat(rinc, idat,jdat) - endif + call w3movdat(rinc, idat,jdat) ! jdate(8)- date and time (yr, mo, day, [tz], hr, min, sec) jdow = 0 jdoy = 0 diff --git a/physics/h2ointerp.f90 b/physics/h2ointerp.f90 index f26ae6c0c..de5eb855d 100644 --- a/physics/h2ointerp.f90 +++ b/physics/h2ointerp.f90 @@ -146,8 +146,6 @@ subroutine h2ointerpol(me,npts,idate,fhour,jindx1,jindx2,h2oplout,ddy) real(kind=kind_phys) h2oplout(npts,levh2o,h2o_coeff) real(kind=kind_phys) rinc(5), rjday integer jdow, jdoy, jday - real(4) rinc4(5) - integer w3kindreal, w3kindint ! idat = 0 idat(1) = idate(4) @@ -156,13 +154,7 @@ subroutine h2ointerpol(me,npts,idate,fhour,jindx1,jindx2,h2oplout,ddy) idat(5) = idate(1) rinc = 0. rinc(2) = fhour - call w3kind(w3kindreal,w3kindint) - if(w3kindreal==4) then - rinc4 = rinc - CALL W3MOVDAT(RINC4,IDAT,JDAT) - else - CALL W3MOVDAT(RINC,IDAT,JDAT) - endif + CALL W3MOVDAT(RINC,IDAT,JDAT) ! jdow = 0 jdoy = 0 diff --git a/physics/iccninterp.F90 b/physics/iccninterp.F90 index a3a08dee8..37e0bfb00 100644 --- a/physics/iccninterp.F90 +++ b/physics/iccninterp.F90 @@ -23,6 +23,7 @@ SUBROUTINE read_cidata (me, master) integer, intent(in) :: me integer, intent(in) :: master !--- locals + integer :: ncerr integer :: i, n, k, ncid, varid,j,it real(kind=kind_phys), allocatable, dimension(:) :: hyam,hybm real(kind=4), allocatable, dimension(:,:,:) :: ci_ps @@ -31,29 +32,29 @@ SUBROUTINE read_cidata (me, master) allocate (ciplin(lonscip,latscip,kcipl,timeci)) allocate (ccnin(lonscip,latscip,kcipl,timeci)) allocate (ci_pres(lonscip,latscip,kcipl,timeci)) - call nf_open("cam5_4_143_NAAI_monclimo2.nc", NF90_NOWRITE, ncid) - call nf_inq_varid(ncid, "lat", varid) - call nf_get_var(ncid, varid, ci_lat) - call nf_inq_varid(ncid, "lon", varid) - call nf_get_var(ncid, varid, ci_lon) - call nf_inq_varid(ncid, "PS", varid) - call nf_get_var(ncid, varid, ci_ps) - call nf_inq_varid(ncid, "hyam", varid) - call nf_get_var(ncid, varid, hyam) - call nf_inq_varid(ncid, "hybm", varid) - call nf_get_var(ncid, varid, hybm) - call nf_inq_varid(ncid, "NAAI", varid) - call nf_get_var(ncid, varid, ciplin) + ncerr = nf90_open("cam5_4_143_NAAI_monclimo2.nc", NF90_NOWRITE, ncid) + ncerr = nf90_inq_varid(ncid, "lat", varid) + ncerr = nf90_get_var(ncid, varid, ci_lat) + ncerr = nf90_inq_varid(ncid, "lon", varid) + ncerr = nf90_get_var(ncid, varid, ci_lon) + ncerr = nf90_inq_varid(ncid, "PS", varid) + ncerr = nf90_get_var(ncid, varid, ci_ps) + ncerr = nf90_inq_varid(ncid, "hyam", varid) + ncerr = nf90_get_var(ncid, varid, hyam) + ncerr = nf90_inq_varid(ncid, "hybm", varid) + ncerr = nf90_get_var(ncid, varid, hybm) + ncerr = nf90_inq_varid(ncid, "NAAI", varid) + ncerr = nf90_get_var(ncid, varid, ciplin) do it = 1,timeci do k=1, kcipl ci_pres(:,:,k,it)=hyam(k)*1.e5+hybm(k)*ci_ps(:,:,it) end do end do - call nf_close(ncid) - call nf_open("cam5_4_143_NPCCN_monclimo2.nc", NF90_NOWRITE, ncid) - call nf_inq_varid(ncid, "NPCCN", varid) - call nf_get_var(ncid, varid, ccnin) - call nf_close(ncid) + ncerr = nf90_close(ncid) + ncerr = nf90_open("cam5_4_143_NPCCN_monclimo2.nc", NF90_NOWRITE, ncid) + ncerr = nf90_inq_varid(ncid, "NPCCN", varid) + ncerr = nf90_get_var(ncid, varid, ccnin) + ncerr = nf90_close(ncid) !--- deallocate (hyam, hybm, ci_ps) if (me == master) then @@ -145,8 +146,6 @@ SUBROUTINE ciinterpol(me,npts,IDATE,FHOUR,jindx1,jindx2,ddy, & real(kind=kind_phys) cipres(npts,kcipl), prsl(npts,lev) real(kind=kind_phys) RINC(5), rjday integer jdow, jdoy, jday - real(4) rinc4(5) - integer w3kindreal,w3kindint ! IDAT=0 IDAT(1)=IDATE(4) @@ -155,13 +154,7 @@ SUBROUTINE ciinterpol(me,npts,IDATE,FHOUR,jindx1,jindx2,ddy, & IDAT(5)=IDATE(1) RINC=0. RINC(2)=FHOUR - call w3kind(w3kindreal,w3kindint) - if(w3kindreal==4) then - rinc4=rinc - CALL W3MOVDAT(RINC4,IDAT,JDAT) - else - CALL W3MOVDAT(RINC,IDAT,JDAT) - endif + CALL W3MOVDAT(RINC,IDAT,JDAT) ! jdow = 0 jdoy = 0 diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index c23b6d1d8..75f199305 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -64,7 +64,7 @@ MODULE module_mp_thompson USE module_mp_radar #ifdef MPI - use mpi + use mpi_f08 #endif IMPLICIT NONE @@ -419,7 +419,7 @@ MODULE module_mp_thompson REAL:: t1_qs_me, t2_qs_me, t1_qg_me, t2_qg_me !..MPI communicator - INTEGER:: mpi_communicator + TYPE(MPI_Comm):: mpi_communicator !..Write tables with master MPI task after computing them in thompson_init LOGICAL:: thompson_table_writer @@ -444,7 +444,8 @@ SUBROUTINE thompson_init(is_aerosol_aware_in, & IMPLICIT NONE LOGICAL, INTENT(IN) :: is_aerosol_aware_in - INTEGER, INTENT(IN) :: mpicomm, mpirank, mpiroot + TYPE(MPI_Comm), INTENT(IN) :: mpicomm + INTEGER, INTENT(IN) :: mpirank, mpiroot INTEGER, INTENT(IN) :: threads CHARACTER(len=*), INTENT(INOUT) :: errmsg INTEGER, INTENT(INOUT) :: errflg @@ -1810,7 +1811,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & qgten1, qiten1, niten1, nrten1, ncten1, qcten1) #ifdef MPI - use mpi + use mpi_f08 #endif implicit none diff --git a/physics/mp_fer_hires.meta b/physics/mp_fer_hires.meta index 9f7c63d4d..058154856 100644 --- a/physics/mp_fer_hires.meta +++ b/physics/mp_fer_hires.meta @@ -55,7 +55,7 @@ long_name = MPI communicator units = index dimensions = () - type = integer + type = MPI_Comm intent = in [mpirank] standard_name = mpi_rank diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 7c76ea933..eb6166b3d 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -6,6 +6,7 @@ !! This module contains the aerosol-aware Thompson microphysics scheme. module mp_thompson + use mpi_f08 use machine, only : kind_phys use module_mp_thompson, only : thompson_init, mp_gt_driver, thompson_finalize, calc_effectRad @@ -72,7 +73,7 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & real(kind_phys), intent(in ) :: phil(:,:) real(kind_phys), intent(in ) :: area(:) ! MPI information - integer, intent(in ) :: mpicomm + type(MPI_Comm), intent(in ) :: mpicomm integer, intent(in ) :: mpirank integer, intent(in ) :: mpiroot ! Threading/blocking information @@ -362,7 +363,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & integer, intent(in) :: decfl ! MPI and block information integer, intent(in) :: blkno - integer, intent(in) :: mpicomm + type(MPI_Comm), intent(in) :: mpicomm integer, intent(in) :: mpirank integer, intent(in) :: mpiroot ! Extended diagnostic output diff --git a/physics/mp_thompson.meta b/physics/mp_thompson.meta index a3bc20615..c728cd14a 100644 --- a/physics/mp_thompson.meta +++ b/physics/mp_thompson.meta @@ -221,7 +221,7 @@ long_name = MPI communicator units = index dimensions = () - type = integer + type = MPI_Comm intent = in [mpirank] standard_name = mpi_rank @@ -593,7 +593,7 @@ long_name = MPI communicator units = index dimensions = () - type = integer + type = MPI_Comm intent = in [mpirank] standard_name = mpi_rank diff --git a/physics/mp_thompson_post.F90 b/physics/mp_thompson_post.F90 index c53f61b0c..9ab841928 100644 --- a/physics/mp_thompson_post.F90 +++ b/physics/mp_thompson_post.F90 @@ -1,5 +1,6 @@ module mp_thompson_post + use mpi_f08 use machine, only : kind_phys implicit none @@ -66,7 +67,7 @@ subroutine mp_thompson_post_run(ncol, nlev, tgrs_save, tgrs, prslk, dtp, ttendli real(kind_phys), intent(in) :: ttendlim integer, intent(in) :: kdt ! MPI information - integer, intent(in ) :: mpicomm + type(MPI_Comm), intent(in ) :: mpicomm integer, intent(in ) :: mpirank integer, intent(in ) :: mpiroot ! CCPP error handling diff --git a/physics/mp_thompson_post.meta b/physics/mp_thompson_post.meta index 82b035e99..7b843000e 100644 --- a/physics/mp_thompson_post.meta +++ b/physics/mp_thompson_post.meta @@ -101,7 +101,7 @@ long_name = MPI communicator units = index dimensions = () - type = integer + type = MPI_Comm intent = in [mpirank] standard_name = mpi_rank diff --git a/physics/ozinterp.f90 b/physics/ozinterp.f90 index 6fe86c8e1..2f0a28128 100644 --- a/physics/ozinterp.f90 +++ b/physics/ozinterp.f90 @@ -149,8 +149,6 @@ SUBROUTINE ozinterpol(me,npts,IDATE,FHOUR,jindx1,jindx2,ozplout,ddy) real(kind=kind_phys) ozplout(npts,levozp,oz_coeff) real(kind=kind_phys) RINC(5), rjday integer jdow, jdoy, jday - real(4) rinc4(5) - integer w3kindreal,w3kindint ! IDAT=0 IDAT(1)=IDATE(4) @@ -159,13 +157,7 @@ SUBROUTINE ozinterpol(me,npts,IDATE,FHOUR,jindx1,jindx2,ozplout,ddy) IDAT(5)=IDATE(1) RINC=0. RINC(2)=FHOUR - call w3kind(w3kindreal,w3kindint) - if(w3kindreal==4) then - rinc4=rinc - CALL W3MOVDAT(RINC4,IDAT,JDAT) - else - CALL W3MOVDAT(RINC,IDAT,JDAT) - endif + CALL W3MOVDAT(RINC,IDAT,JDAT) ! jdow = 0 jdoy = 0 diff --git a/physics/rrtmgp_lw_cloud_optics.F90 b/physics/rrtmgp_lw_cloud_optics.F90 index 5ddcec078..28f1f5fb4 100644 --- a/physics/rrtmgp_lw_cloud_optics.F90 +++ b/physics/rrtmgp_lw_cloud_optics.F90 @@ -8,7 +8,7 @@ module rrtmgp_lw_cloud_optics use radiation_tools, only: check_error_msg use netcdf #ifdef MPI - use mpi + use mpi_f08 #endif implicit none @@ -81,8 +81,9 @@ subroutine rrtmgp_lw_cloud_optics_init(nrghice, mpicomm, mpirank, mpiroot, doGP_cldoptics_LUT ! Use RRTMGP cloud-optics: LUTs? integer, intent(inout) :: & nrghice ! Number of ice-roughness categories + type(MPI_Comm), intent(in) :: & + mpicomm ! MPI communicator integer, intent(in) :: & - mpicomm, & ! MPI communicator mpirank, & ! Current MPI rank mpiroot ! Master MPI rank character(len=128),intent(in) :: & diff --git a/physics/rrtmgp_lw_cloud_optics.meta b/physics/rrtmgp_lw_cloud_optics.meta index 35e27979e..3d34884cb 100644 --- a/physics/rrtmgp_lw_cloud_optics.meta +++ b/physics/rrtmgp_lw_cloud_optics.meta @@ -70,7 +70,7 @@ long_name = MPI communicator units = index dimensions = () - type = integer + type = MPI_Comm intent = in [errmsg] standard_name = ccpp_error_message diff --git a/physics/rrtmgp_lw_gas_optics.F90 b/physics/rrtmgp_lw_gas_optics.F90 index 67a888911..af47d86ab 100644 --- a/physics/rrtmgp_lw_gas_optics.F90 +++ b/physics/rrtmgp_lw_gas_optics.F90 @@ -8,7 +8,7 @@ module rrtmgp_lw_gas_optics use radiation_tools, only: check_error_msg use netcdf #ifdef MPI - use mpi + use mpi_f08 #endif implicit none @@ -82,8 +82,9 @@ subroutine rrtmgp_lw_gas_optics_init(rrtmgp_root_dir, rrtmgp_lw_file_gas, mpicom character(len=128),intent(in) :: & rrtmgp_root_dir, & ! RTE-RRTMGP root directory rrtmgp_lw_file_gas ! RRTMGP file containing coefficients used to compute gaseous optical properties + type(MPI_Comm),intent(in) :: & + mpicomm ! MPI communicator integer,intent(in) :: & - mpicomm, & ! MPI communicator mpirank, & ! Current MPI rank mpiroot ! Master MPI rank diff --git a/physics/rrtmgp_lw_gas_optics.meta b/physics/rrtmgp_lw_gas_optics.meta index 0b484b6ac..37b44f502 100644 --- a/physics/rrtmgp_lw_gas_optics.meta +++ b/physics/rrtmgp_lw_gas_optics.meta @@ -42,7 +42,7 @@ long_name = MPI communicator units = index dimensions = () - type = integer + type = MPI_Comm intent = in [minGPpres] standard_name = minimum_pressure_in_RRTMGP diff --git a/physics/rrtmgp_sw_cloud_optics.F90 b/physics/rrtmgp_sw_cloud_optics.F90 index f80440522..2a175a0e6 100644 --- a/physics/rrtmgp_sw_cloud_optics.F90 +++ b/physics/rrtmgp_sw_cloud_optics.F90 @@ -8,7 +8,7 @@ module rrtmgp_sw_cloud_optics use radiation_tools, only: check_error_msg use netcdf #ifdef MPI - use mpi + use mpi_f08 #endif implicit none @@ -81,8 +81,9 @@ subroutine rrtmgp_sw_cloud_optics_init(doG_cldoptics, doGP_cldoptics_PADE, doGP_cldoptics_LUT ! Use RRTMGP cloud-optics: LUTs? integer, intent(inout) :: & nrghice ! Number of ice-roughness categories + type(MPI_Comm), intent(in) :: & + mpicomm ! MPI communicator integer, intent(in) :: & - mpicomm, & ! MPI communicator mpirank, & ! Current MPI rank mpiroot ! Master MPI rank character(len=128),intent(in) :: & diff --git a/physics/rrtmgp_sw_cloud_optics.meta b/physics/rrtmgp_sw_cloud_optics.meta index d73258cb2..f852cf7cb 100644 --- a/physics/rrtmgp_sw_cloud_optics.meta +++ b/physics/rrtmgp_sw_cloud_optics.meta @@ -70,7 +70,7 @@ long_name = MPI communicator units = index dimensions = () - type = integer + type = MPI_Comm intent = in [errmsg] standard_name = ccpp_error_message diff --git a/physics/rrtmgp_sw_gas_optics.F90 b/physics/rrtmgp_sw_gas_optics.F90 index 260f65fe7..9f2078a51 100644 --- a/physics/rrtmgp_sw_gas_optics.F90 +++ b/physics/rrtmgp_sw_gas_optics.F90 @@ -7,7 +7,7 @@ module rrtmgp_sw_gas_optics use mo_optical_props, only: ty_optical_props_2str use netcdf #ifdef MPI - use mpi + use mpi_f08 #endif implicit none @@ -86,8 +86,9 @@ subroutine rrtmgp_sw_gas_optics_init(rrtmgp_root_dir, rrtmgp_sw_file_gas, character(len=128),intent(in) :: & rrtmgp_root_dir, & ! RTE-RRTMGP root directory rrtmgp_sw_file_gas ! RRTMGP file containing coefficients used to compute gaseous optical properties + type(MPI_Comm),intent(in) :: & + mpicomm ! MPI communicator integer,intent(in) :: & - mpicomm, & ! MPI communicator mpirank, & ! Current MPI rank mpiroot ! Master MPI rank character(len=*), dimension(:), intent(in) :: & diff --git a/physics/rrtmgp_sw_gas_optics.meta b/physics/rrtmgp_sw_gas_optics.meta index 1fdbc946b..8a018172d 100644 --- a/physics/rrtmgp_sw_gas_optics.meta +++ b/physics/rrtmgp_sw_gas_optics.meta @@ -50,7 +50,7 @@ long_name = MPI communicator units = index dimensions = () - type = integer + type = MPI_Comm intent = in [errmsg] standard_name = ccpp_error_message diff --git a/physics/sfcsub.F b/physics/sfcsub.F index e8b61f083..321e87a98 100644 --- a/physics/sfcsub.F +++ b/physics/sfcsub.F @@ -2740,7 +2740,7 @@ subroutine fixrdg(lugb,idim,jdim,fngrib, & use sfccyc_module, only : mdata implicit none integer lgrib,n,lskip,jret,j,ndata,lugi,jdim,idim,lugb, - & iret, me,kpds5,kdata,i,w3kindreal,w3kindint + & iret, me,kpds5,kdata,i ! character*(*) fngrib ! @@ -2748,7 +2748,6 @@ subroutine fixrdg(lugb,idim,jdim,fngrib, & logical gaus real (kind=kind_io8) blno,blto real (kind=kind_io8), allocatable :: data8(:) - real (kind=kind_io4), allocatable :: data4(:) ! logical*1, allocatable :: lbms(:) ! @@ -2807,20 +2806,8 @@ subroutine fixrdg(lugb,idim,jdim,fngrib, & jpds = kpds0 lskip = -1 kdata=idim*jdim - call w3kind(w3kindreal,w3kindint) - if (w3kindreal == 8) then call getgb(lugb,lugi,kdata,lskip,jpds,jgds,ndata,lskip, & kpds,kgds,lbms,data8,jret) - else if (w3kindreal == 4) then - allocate(data4(1:idim*jdim)) - call getgb(lugb,lugi,kdata,lskip,jpds,jgds,ndata,lskip, - & kpds,kgds,lbms,data4,jret) - data8 = real(data4, kind=kind_io8) - deallocate(data4) - else - write(0,*)' FATAL ERROR: Invalid w3kindreal' - call abort - endif ! if(jret == 0) then if(ndata.eq.0) then @@ -7041,8 +7028,6 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil,slmskl,slmskw, & & 196.5,227.5,258.0,288.5,319.0,349.5,380.5/ ! real (kind=kind_io8) fha(5) - real(4) fha4(5) - integer w3kindreal,w3kindint integer ida(8),jda(8),ivtyp, kpd7 ! real (kind=kind_io8), allocatable :: tsf(:,:),sno(:,:), @@ -7134,13 +7119,7 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil,slmskl,slmskw, & ida(2) = im ida(3) = id ida(5) = ih - call w3kind(w3kindreal,w3kindint) - if(w3kindreal == 4) then - fha4 = fha - call w3movdat(fha4,ida,jda) - else - call w3movdat(fha,ida,jda) - endif + call w3movdat(fha,ida,jda) jy = jda(1) jm = jda(2) jd = jda(3) @@ -7210,13 +7189,7 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil,slmskl,slmskw, & ida(2) = im ida(3) = id ida(5) = ih - call w3kind(w3kindreal,w3kindint) - if(w3kindreal == 4) then - fha4 = fha - call w3movdat(fha4,ida,jda) - else - call w3movdat(fha,ida,jda) - endif + call w3movdat(fha,ida,jda) jy = jda(1) jm = jda(2) jd = jda(3) @@ -8310,7 +8283,7 @@ subroutine fixrdc(lugb,fngrib,kpds5,kpds7,mon,slmask, & implicit none integer imax,jmax,ijmax,i,j,n,jret,inttyp,iret,imsk, & & jmsk,len,lugb,kpds5,mon,lskip,lgrib,ndata,lugi,me,kmami & - &, jj,w3kindreal,w3kindint + &, jj real (kind=kind_io8) wlon,elon,rnlat,dlat,dlon,rslat,blno,blto ! ! @@ -8322,7 +8295,6 @@ subroutine fixrdc(lugb,fngrib,kpds5,kpds7,mon,slmask, & real (kind=kind_io8) gdata(len), slmask(len) real (kind=kind_io8), allocatable :: data(:,:), rslmsk(:,:) real (kind=kind_io8), allocatable :: data8(:) - real (kind=kind_io4), allocatable :: data4(:) real (kind=kind_io8), allocatable :: rlngrb(:), rltgrb(:) ! logical lmask, yr2kc, gaus, ijordr @@ -8390,17 +8362,8 @@ subroutine fixrdc(lugb,fngrib,kpds5,kpds7,mon,slmask, & jpds = kpds0 jpds(9) = mon if(jpds(9).eq.13) jpds(9) = 1 - call w3kind(w3kindreal,w3kindint) - if (w3kindreal==8) then - call getgb(lugb,lugi,mdata,lskip,jpds,jgds,ndata,lskip, - & kpds,kgds,lbms,data8,jret) - else if (w3kindreal==4) then - allocate(data4(1:mdata)) - call getgb(lugb,lugi,mdata,lskip,jpds,jgds,ndata,lskip, - & kpds,kgds,lbms,data4,jret) - data8 = real(data4, kind=kind_io8) - deallocate(data4) - endif + call getgb(lugb,lugi,mdata,lskip,jpds,jgds,ndata,lskip, + & kpds,kgds,lbms,data8,jret) if (me .eq. 0) write(6,*) ' input grib file dates=', & (kpds(i),i=8,11) if(jret.eq.0) then @@ -8485,7 +8448,7 @@ subroutine fixrda(lugb,fngrib,kpds5,slmask, & integer nrepmx,nvalid,imo,iyr,idy,jret,ihr,nrept,lskip,lugi, & & lgrib,j,ndata,i,inttyp,jmax,imax,ijmax,ij,jday,len,iret, & & jmsk,imsk,ih,kpds5,lugb,iy,id,im,jh,jd,jdoy,jdow,jm,me, & - & monend,jy,iy4,kmami,iret2,jj,w3kindreal,w3kindint + & monend,jy,iy4,kmami,iret2,jj real (kind=kind_io8) rnlat,rslat,wlon,elon,dlon,dlat,fh,blno, & & rjday,blto ! @@ -8507,7 +8470,6 @@ subroutine fixrda(lugb,fngrib,kpds5,slmask, & real (kind=kind_io8) gdata(len), slmask(len) real (kind=kind_io8), allocatable :: data(:,:),rslmsk(:,:) real (kind=kind_io8), allocatable :: data8(:) - real (kind=kind_io4), allocatable :: data4(:) real (kind=kind_io8), allocatable :: rlngrb(:), rltgrb(:) ! logical lmask, yr2kc, gaus, ijordr @@ -8529,7 +8491,6 @@ subroutine fixrda(lugb,fngrib,kpds5,slmask, & data mjday/31,28,31,30,31,30,31,31,30,31,30,31/ ! real (kind=kind_io8) fha(5) - real(4) fha4(5) integer ida(8),jda(8) ! allocate(data8(1:mdata)) @@ -8548,13 +8509,7 @@ subroutine fixrda(lugb,fngrib,kpds5,slmask, & ida(2)=im ida(3)=id ida(5)=ih - call w3kind(w3kindreal,w3kindint) - if(w3kindreal==4) then - fha4=fha - call w3movdat(fha4,ida,jda) - else - call w3movdat(fha,ida,jda) - endif + call w3movdat(fha,ida,jda) jy=jda(1) jm=jda(2) jd=jda(3) @@ -8637,17 +8592,8 @@ subroutine fixrda(lugb,fngrib,kpds5,slmask, & jpds(10)=idy ! jpds(11)=ihr jpds(21)=(iyr-1)/100+1 - call w3kind(w3kindreal,w3kindint) - if (w3kindreal == 8) then - call getgb(lugb,lugi,mdata,lskip,jpds,jgds,ndata,lskip, - & kpds,kgds,lbms,data8,jret) - elseif (w3kindreal == 4) then - allocate (data4(1:mdata)) - call getgb(lugb,lugi,mdata,lskip,jpds,jgds,ndata,lskip, - & kpds,kgds,lbms,data4,jret) - data8 = real(data4, kind=kind_io8) - deallocate(data4) - endif + call getgb(lugb,lugi,mdata,lskip,jpds,jgds,ndata,lskip, + & kpds,kgds,lbms,data8,jret) if (me .eq. 0) write(6,*) ' input grib file dates=', & (kpds(i),i=8,11) if(jret.eq.0) then From 16ef7f241efb568da6b860e27a7aac3e9b423036 Mon Sep 17 00:00:00 2001 From: Dusan Jovic Date: Tue, 19 Jul 2022 21:14:44 +0000 Subject: [PATCH 02/21] Declare variables passed to w3nco library as double precission --- physics/aerinterp.F90 | 7 ++++--- physics/h2ointerp.f90 | 3 ++- physics/iccninterp.F90 | 3 ++- physics/ozinterp.f90 | 3 ++- physics/sfcsub.F | 12 ++++++------ 5 files changed, 16 insertions(+), 12 deletions(-) diff --git a/physics/aerinterp.F90 b/physics/aerinterp.F90 index 90765e4d5..790aa0096 100644 --- a/physics/aerinterp.F90 +++ b/physics/aerinterp.F90 @@ -98,7 +98,7 @@ END SUBROUTINE read_aerdata ! !********************************************************************** SUBROUTINE read_aerdataf ( me, master, iflip, idate, FHOUR, errmsg, errflg) - use machine, only: kind_phys, kind_io4, kind_io8 + use machine, only: kind_phys, kind_dbl_prec use aerclm_def !--- in/out @@ -111,7 +111,7 @@ SUBROUTINE read_aerdataf ( me, master, iflip, idate, FHOUR, errmsg, errflg) logical :: file_exist integer IDAT(8),JDAT(8) real(kind=kind_phys) rjday - real(8) RINC(5) + real(kind=kind_dbl_prec) rinc(5) integer jdow, jdoy, jday integer, allocatable :: invardims(:) @@ -218,7 +218,7 @@ END SUBROUTINE setindxaer SUBROUTINE aerinterpol( me,master,nthrds,npts,IDATE,FHOUR,iflip, jindx1,jindx2, & ddy,iindx1,iindx2,ddx,lev,prsl,aerout) ! - use machine, only: kind_phys, kind_io4, kind_io8 + use machine, only: kind_phys, kind_dbl_prec use aerclm_def implicit none @@ -238,6 +238,7 @@ SUBROUTINE aerinterpol( me,master,nthrds,npts,IDATE,FHOUR,iflip, jindx1,jindx2, real(kind=kind_phys) aerpm(npts,levsaer,ntrcaer) real(kind=kind_phys) prsl(npts,lev), aerpres(npts,levsaer) real(kind=kind_phys) rjday + real(kind=kind_dbl_prec) rinc(5) integer jdow, jdoy, jday ! IDAT = 0 diff --git a/physics/h2ointerp.f90 b/physics/h2ointerp.f90 index 2d16decc0..f5a1f36c6 100644 --- a/physics/h2ointerp.f90 +++ b/physics/h2ointerp.f90 @@ -131,7 +131,7 @@ subroutine h2ointerpol(me,npts,idate,fhour,jindx1,jindx2,h2oplout,ddy) ! ! May 2015 Shrinivas Moorthi - Prepare for H2O interpolation ! - use machine , only : kind_phys + use machine , only : kind_phys, kind_dbl_prec use h2o_def implicit none integer j,j1,j2,l,npts,nc,n1,n2 @@ -145,6 +145,7 @@ subroutine h2ointerpol(me,npts,idate,fhour,jindx1,jindx2,h2oplout,ddy) real(kind=kind_phys) ddy(npts) real(kind=kind_phys) h2oplout(npts,levh2o,h2o_coeff) real(kind=kind_phys) rjday + real(kind=kind_dbl_prec) rinc(5) integer jdow, jdoy, jday ! idat = 0 diff --git a/physics/iccninterp.F90 b/physics/iccninterp.F90 index 1d613c53c..dd752d9b8 100644 --- a/physics/iccninterp.F90 +++ b/physics/iccninterp.F90 @@ -129,7 +129,7 @@ END SUBROUTINE setindxci SUBROUTINE ciinterpol(me,npts,IDATE,FHOUR,jindx1,jindx2,ddy, & iindx1,iindx2,ddx,lev, prsl, ciplout,ccnout) ! - USE MACHINE, ONLY : kind_phys + USE MACHINE, ONLY : kind_phys, kind_dbl_prec use iccn_def implicit none integer i1,i2, iday,j,j1,j2,l,npts,nc,n1,n2,lev,k,i @@ -145,6 +145,7 @@ SUBROUTINE ciinterpol(me,npts,IDATE,FHOUR,jindx1,jindx2,ddy, & real(kind=kind_phys) ccnout(npts,lev),ccnpm(npts,kcipl) real(kind=kind_phys) cipres(npts,kcipl), prsl(npts,lev) real(kind=kind_phys) rjday + real(kind=kind_dbl_prec) rinc(5) integer jdow, jdoy, jday ! IDAT=0 diff --git a/physics/ozinterp.f90 b/physics/ozinterp.f90 index f92140eb1..e26d5a56e 100644 --- a/physics/ozinterp.f90 +++ b/physics/ozinterp.f90 @@ -135,7 +135,7 @@ END SUBROUTINE setindxoz ! SUBROUTINE ozinterpol(me,npts,IDATE,FHOUR,jindx1,jindx2,ozplout,ddy) ! - USE MACHINE, ONLY : kind_phys + USE MACHINE, ONLY : kind_phys, kind_dbl_prec USE OZNE_DEF implicit none integer iday,j,j1,j2,l,npts,nc,n1,n2 @@ -148,6 +148,7 @@ SUBROUTINE ozinterpol(me,npts,IDATE,FHOUR,jindx1,jindx2,ozplout,ddy) real(kind=kind_phys) DDY(npts) real(kind=kind_phys) ozplout(npts,levozp,oz_coeff) real(kind=kind_phys) rjday + real(kind=kind_dbl_prec) rinc(5) integer jdow, jdoy, jday ! IDAT=0 diff --git a/physics/sfcsub.F b/physics/sfcsub.F index d5a1d997d..4c7168b90 100644 --- a/physics/sfcsub.F +++ b/physics/sfcsub.F @@ -2747,7 +2747,7 @@ subroutine fixrdg(lugb,idim,jdim,fngrib, & real (kind=kind_io8) gdata(idim*jdim) logical gaus real (kind=kind_io8) blno,blto - real (kind=kind_io8), allocatable :: data8(:) + real (kind=kind_dbl_prec), allocatable :: data8(:) ! logical*1, allocatable :: lbms(:) ! @@ -6969,7 +6969,7 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil,slmskl,slmskw, & &, gaus, blno, blto, me,lprnt,iprnt, fnalbc2, ialb & &, tile_num_ch, i_index, j_index) ! - use machine , only : kind_io8,kind_io4 + use machine , only : kind_io8,kind_io4, kind_dbl_prec implicit none character(len=*), intent(in) :: tile_num_ch integer, intent(in) :: i_index(len), j_index(len) @@ -7029,7 +7029,7 @@ subroutine clima(lugb,iy,im,id,ih,fh,len,lsoil,slmskl,slmskw, & data dayhf/ 15.5, 45.0, 74.5,105.0,135.5,166.0, & 196.5,227.5,258.0,288.5,319.0,349.5,380.5/ ! - real (kind=kind_io8) fha(5) + real (kind=kind_dbl_prec) fha(5) integer ida(8),jda(8),ivtyp, kpd7 ! real (kind=kind_io8), allocatable :: tsf(:,:),sno(:,:), @@ -8296,7 +8296,7 @@ subroutine fixrdc(lugb,fngrib,kpds5,kpds7,mon,slmask, & ! real (kind=kind_io8) gdata(len), slmask(len) real (kind=kind_io8), allocatable :: data(:,:), rslmsk(:,:) - real (kind=kind_io8), allocatable :: data8(:) + real (kind=kind_dbl_prec), allocatable :: data8(:) real (kind=kind_io8), allocatable :: rlngrb(:), rltgrb(:) ! logical lmask, yr2kc, gaus, ijordr @@ -8471,7 +8471,7 @@ subroutine fixrda(lugb,fngrib,kpds5,slmask, & ! real (kind=kind_io8) gdata(len), slmask(len) real (kind=kind_io8), allocatable :: data(:,:),rslmsk(:,:) - real (kind=kind_io8), allocatable :: data8(:) + real (kind=kind_dbl_prec), allocatable :: data8(:) real (kind=kind_io8), allocatable :: rlngrb(:), rltgrb(:) ! logical lmask, yr2kc, gaus, ijordr @@ -8492,7 +8492,7 @@ subroutine fixrda(lugb,fngrib,kpds5,slmask, & integer mjday(12) data mjday/31,28,31,30,31,30,31,31,30,31,30,31/ ! - real (kind=kind_io8) fha(5) + real (kind=kind_dbl_prec) fha(5) integer ida(8),jda(8) ! allocate(data8(1:mdata)) From 696c0ea32d86259836016977621d05c316db0d73 Mon Sep 17 00:00:00 2001 From: Dusan Jovic Date: Wed, 20 Jul 2022 11:48:43 -0500 Subject: [PATCH 03/21] More argument mismatch fixes --- physics/GFS_time_vary_pre.fv3.F90 | 18 +++--------------- physics/GFS_time_vary_pre.scm.F90 | 18 +++--------------- 2 files changed, 6 insertions(+), 30 deletions(-) diff --git a/physics/GFS_time_vary_pre.fv3.F90 b/physics/GFS_time_vary_pre.fv3.F90 index 86dfe3f40..ba661f578 100644 --- a/physics/GFS_time_vary_pre.fv3.F90 +++ b/physics/GFS_time_vary_pre.fv3.F90 @@ -92,10 +92,8 @@ subroutine GFS_time_vary_pre_timestep_init (jdat, idat, dtp, nsswr, real(kind=kind_phys), parameter :: con_24 = 24.0_kind_phys real(kind=kind_phys), parameter :: con_hr = 3600.0_kind_phys - real(kind=kind_sngl_prec) :: rinc4(5) real(kind=kind_dbl_prec) :: rinc8(5) - integer :: w3kindreal,w3kindint integer :: iw3jdn integer :: jd0, jd1 real :: fjd @@ -113,19 +111,9 @@ subroutine GFS_time_vary_pre_timestep_init (jdat, idat, dtp, nsswr, !--- jdat is being updated directly inside of FV3GFS_cap.F90 !--- update calendars and triggers - call w3kind(w3kindreal,w3kindint) - if (w3kindreal == 8) then - rinc8(1:5) = 0 - call w3difdat(jdat,idat,4,rinc8) - sec = rinc8(4) - else if (w3kindreal == 4) then - rinc4(1:5) = 0 - call w3difdat(jdat,idat,4,rinc4) - sec = rinc4(4) - else - write(0,*)' FATAL ERROR: Invalid w3kindreal' - call abort - endif + rinc8(1:5) = 0 + call w3difdat(jdat,idat,4,rinc8) + sec = rinc8(4) phour = sec/con_hr !--- set current bucket hour zhour = phour diff --git a/physics/GFS_time_vary_pre.scm.F90 b/physics/GFS_time_vary_pre.scm.F90 index 2bb6b3ceb..3293e09e4 100644 --- a/physics/GFS_time_vary_pre.scm.F90 +++ b/physics/GFS_time_vary_pre.scm.F90 @@ -91,10 +91,8 @@ subroutine GFS_time_vary_pre_timestep_init (jdat, idat, dtp, nsswr, & real(kind=kind_phys), parameter :: con_24 = 24.0_kind_phys real(kind=kind_phys), parameter :: con_hr = 3600.0_kind_phys - real(kind=kind_sngl_prec) :: rinc4(5) real(kind=kind_dbl_prec) :: rinc8(5) - integer :: w3kindreal,w3kindint integer :: iw3jdn integer :: jd0, jd1 real :: fjd @@ -114,19 +112,9 @@ subroutine GFS_time_vary_pre_timestep_init (jdat, idat, dtp, nsswr, & !--- jdat is being updated directly inside of the time integration !--- loop of scm.F90 !--- update calendars and triggers - call w3kind(w3kindreal,w3kindint) - if (w3kindreal == 8) then - rinc8(1:5) = 0 - call w3difdat(jdat,idat,4,rinc8) - sec = rinc8(4) - else if (w3kindreal == 4) then - rinc4(1:5) = 0 - call w3difdat(jdat,idat,4,rinc4) - sec = rina4c(4) - else - write(0,*)' FATAL ERROR: Invalid w3kindreal' - call abort - endif + rinc8(1:5) = 0 + call w3difdat(jdat,idat,4,rinc8) + sec = rinc8(4) phour = sec/con_hr !--- set current bucket hour zhour = phour From 92b8c857bd211770e53b2b3b80ecdf5115ba7e57 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 7 Nov 2022 11:35:12 -0700 Subject: [PATCH 04/21] If MPI is used, find package --- CMakeLists.txt | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index d14778b06..4de7aa24e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -8,6 +8,12 @@ project(ccpp_physics set(PACKAGE "ccpp-physics") set(AUTHORS "Grant Firl" "Dom Heinzeller" "Man Zhang" "Mike Kavulich" "Chunxi Zhang") +#------------------------------------------------------------------------------ +# Set MPI flags for C/C++/Fortran +if (MPI) + find_package(MPI REQUIRED C Fortran) +endif (OPENMP) + #------------------------------------------------------------------------------ # Set OpenMP flags for C/C++/Fortran if (OPENMP) From f943d204fa9bb2200ce8b5ae545919c62c003664 Mon Sep 17 00:00:00 2001 From: Dusan Jovic Date: Fri, 18 Nov 2022 20:59:10 +0000 Subject: [PATCH 05/21] Fixed cmake if/endif mismatch --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 4de7aa24e..8a88fadac 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -12,7 +12,7 @@ set(AUTHORS "Grant Firl" "Dom Heinzeller" "Man Zhang" "Mike Kavulich" "Chunxi Zh # Set MPI flags for C/C++/Fortran if (MPI) find_package(MPI REQUIRED C Fortran) -endif (OPENMP) +endif() #------------------------------------------------------------------------------ # Set OpenMP flags for C/C++/Fortran From 347d74520d1bb8312ad4cceba1ec6c6038b8c650 Mon Sep 17 00:00:00 2001 From: Dusan Jovic Date: Thu, 16 Mar 2023 17:46:04 +0000 Subject: [PATCH 06/21] Use mpi_f08 in rrtmgp_lw_main and rrtmgp_sw_main --- physics/rrtmgp_lw_main.F90 | 4 +++- physics/rrtmgp_lw_main.meta | 4 ++-- physics/rrtmgp_sw_main.F90 | 4 +++- physics/rrtmgp_sw_main.meta | 2 +- 4 files changed, 9 insertions(+), 5 deletions(-) diff --git a/physics/rrtmgp_lw_main.F90 b/physics/rrtmgp_lw_main.F90 index c0bc99d35..f3a2f5ba6 100644 --- a/physics/rrtmgp_lw_main.F90 +++ b/physics/rrtmgp_lw_main.F90 @@ -7,6 +7,7 @@ !! ! ########################################################################################### module rrtmgp_lw_main + use mpi_f08 use machine, only: kind_phys, kind_dbl_prec use mo_optical_props, only: ty_optical_props_1scl, ty_optical_props_2str use mo_cloud_optics, only: ty_cloud_optics @@ -68,8 +69,9 @@ subroutine rrtmgp_lw_main_init(rrtmgp_root_dir, rrtmgp_lw_file_gas, rrtmgp_lw_fi doGP_sgs_cnv ! Flag to include sgs convective clouds integer, intent(inout) :: & nrghice ! Number of ice-roughness categories + type(MPI_Comm),intent(in) :: & + mpicomm ! MPI communicator integer,intent(in) :: & - mpicomm, & ! MPI communicator mpirank, & ! Current MPI rank mpiroot, & ! Master MPI rank rrtmgp_phys_blksz, & ! Number of horizontal points to process at once. diff --git a/physics/rrtmgp_lw_main.meta b/physics/rrtmgp_lw_main.meta index fd96eb14b..ebff48750 100644 --- a/physics/rrtmgp_lw_main.meta +++ b/physics/rrtmgp_lw_main.meta @@ -90,7 +90,7 @@ long_name = MPI communicator units = index dimensions = () - type = integer + type = MPI_Comm intent = in [rrtmgp_phys_blksz] standard_name = number_of_columns_per_RRTMGP_LW_block @@ -638,4 +638,4 @@ units = 1 dimensions = () type = integer - intent = out \ No newline at end of file + intent = out diff --git a/physics/rrtmgp_sw_main.F90 b/physics/rrtmgp_sw_main.F90 index b25e093e7..b1a8e7cc4 100644 --- a/physics/rrtmgp_sw_main.F90 +++ b/physics/rrtmgp_sw_main.F90 @@ -1,6 +1,7 @@ ! ########################################################################################### ! ########################################################################################### module rrtmgp_sw_main + use mpi_f08 use machine, only: kind_phys, kind_dbl_prec use mo_optical_props, only: ty_optical_props_2str use mo_cloud_optics, only: ty_cloud_optics @@ -55,8 +56,9 @@ subroutine rrtmgp_sw_main_init(rrtmgp_root_dir, rrtmgp_sw_file_gas, rrtmgp_sw_fi doGP_sgs_cnv ! Flag to include sgs convective clouds integer, intent(inout) :: & nrghice ! Number of ice-roughness categories + type(MPI_Comm),intent(in) :: & + mpicomm ! MPI communicator integer,intent(in) :: & - mpicomm, & ! MPI communicator mpirank, & ! Current MPI rank mpiroot, & ! Master MPI rank rrtmgp_phys_blksz, & ! Number of horizontal points to process at once. diff --git a/physics/rrtmgp_sw_main.meta b/physics/rrtmgp_sw_main.meta index dbb93a5df..c2b6555b5 100644 --- a/physics/rrtmgp_sw_main.meta +++ b/physics/rrtmgp_sw_main.meta @@ -104,7 +104,7 @@ long_name = MPI communicator units = index dimensions = () - type = integer + type = MPI_Comm intent = in [active_gases_array] standard_name = list_of_active_gases_used_by_RRTMGP From 360a612a4ff606096acb32c09c3e7276c4075224 Mon Sep 17 00:00:00 2001 From: Dusan Jovic Date: Fri, 3 Nov 2023 13:48:20 -0500 Subject: [PATCH 07/21] Update GFS_phys_time_vary.fv3.F90 to compile without '-fallow-argument-mismatch' --- physics/GFS_phys_time_vary.fv3.F90 | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/physics/GFS_phys_time_vary.fv3.F90 b/physics/GFS_phys_time_vary.fv3.F90 index 4b6909f74..2d0b35d79 100644 --- a/physics/GFS_phys_time_vary.fv3.F90 +++ b/physics/GFS_phys_time_vary.fv3.F90 @@ -830,7 +830,6 @@ subroutine GFS_phys_time_vary_timestep_init ( real(kind_phys) :: rannie(cny) real(kind_phys) :: rndval(cnx*cny*nrcm) real(kind_dbl_prec) :: rinc(5) - real(kind_sngl_prec) :: rinc4(5) ! Initialize CCPP error handling variables errmsg = '' @@ -850,7 +849,7 @@ subroutine GFS_phys_time_vary_timestep_init ( !$OMP shared(ozpl,ddy_o3,h2o_phys,jindx1_h,jindx2_h,h2opl,ddy_h,iaerclm,master) & !$OMP shared(levs,prsl,iccn,jindx1_ci,jindx2_ci,ddy_ci,iindx1_ci,iindx2_ci) & !$OMP shared(ddx_ci,in_nm,ccn_nm,do_ugwp_v1,jindx1_tau,jindx2_tau,ddy_j1tau) & -!$OMP shared(ddy_j2tau,tau_amf,iflip,ozphys,rjday,n1,n2,idat,jdat,rinc,rinc4) & +!$OMP shared(ddy_j2tau,tau_amf,iflip,ozphys,rjday,n1,n2,idat,jdat,rinc) & !$OMP shared(w3kindreal,w3kindint,jdow,jdoy,jday) & !$OMP private(iseed,iskip,i,j,k) @@ -910,13 +909,7 @@ subroutine GFS_phys_time_vary_timestep_init ( idat(5)=idate(1) rinc=0. rinc(2)=fhour - call w3kind(w3kindreal,w3kindint) - if(w3kindreal==4) then - rinc4=rinc - CALL w3movdat(rinc4,idat,jdat) - else - CALL w3movdat(rinc,idat,jdat) - endif + CALL w3movdat(rinc,idat,jdat) jdow = 0 jdoy = 0 jday = 0 From bac19946862a08b453e94678727c868f1b60c450 Mon Sep 17 00:00:00 2001 From: Dusan Jovic Date: Wed, 7 Feb 2024 16:47:36 +0000 Subject: [PATCH 08/21] Change the type of mpi communicator in few more files --- physics/MP/Ferrier_Aligo/mp_fer_hires.F90 | 3 ++- physics/MP/NSSL/mp_nssl.F90 | 4 ++-- physics/MP/NSSL/mp_nssl.meta | 2 +- physics/smoke_dust/rrfs_smoke_wrapper.F90 | 3 ++- physics/smoke_dust/rrfs_smoke_wrapper.meta | 2 +- 5 files changed, 8 insertions(+), 6 deletions(-) diff --git a/physics/MP/Ferrier_Aligo/mp_fer_hires.F90 b/physics/MP/Ferrier_Aligo/mp_fer_hires.F90 index 938beae5d..47e80e9d9 100644 --- a/physics/MP/Ferrier_Aligo/mp_fer_hires.F90 +++ b/physics/MP/Ferrier_Aligo/mp_fer_hires.F90 @@ -36,6 +36,7 @@ subroutine mp_fer_hires_init(ncol, nlev, dtp, imp_physics, & mpicomm, mpirank,mpiroot, & threads, errmsg, errflg) + USE mpi_f08 USE machine, ONLY : kind_phys USE MODULE_MP_FER_HIRES, ONLY : FERRIER_INIT_HR implicit none @@ -45,7 +46,7 @@ subroutine mp_fer_hires_init(ncol, nlev, dtp, imp_physics, & real(kind_phys), intent(in) :: dtp integer, intent(in) :: imp_physics integer, intent(in) :: imp_physics_fer_hires - integer, intent(in) :: mpicomm + type(MPI_Comm), intent(in) :: mpicomm integer, intent(in) :: mpirank integer, intent(in) :: mpiroot integer, intent(in) :: threads diff --git a/physics/MP/NSSL/mp_nssl.F90 b/physics/MP/NSSL/mp_nssl.F90 index 0b111f7cd..5c47bce73 100644 --- a/physics/MP/NSSL/mp_nssl.F90 +++ b/physics/MP/NSSL/mp_nssl.F90 @@ -40,7 +40,7 @@ subroutine mp_nssl_init(ncol, nlev, errflg, errmsg, threads, restart, & use module_mp_nssl_2mom, only: nssl_2mom_init, nssl_2mom_init_const #ifdef MPI - use mpi + use mpi_f08 #endif implicit none @@ -56,7 +56,7 @@ subroutine mp_nssl_init(ncol, nlev, errflg, errmsg, threads, restart, & integer, intent(in) :: mpirank integer, intent(in) :: mpiroot - integer, intent(in) :: mpicomm + type(MPI_Comm), intent(in) :: mpicomm integer, intent(in) :: imp_physics integer, intent(in) :: imp_physics_nssl real(kind_phys), intent(in) :: nssl_cccn, nssl_alphah, nssl_alphahl diff --git a/physics/MP/NSSL/mp_nssl.meta b/physics/MP/NSSL/mp_nssl.meta index 8449f26cf..deb96569d 100644 --- a/physics/MP/NSSL/mp_nssl.meta +++ b/physics/MP/NSSL/mp_nssl.meta @@ -68,7 +68,7 @@ long_name = MPI communicator units = index dimensions = () - type = integer + type = MPI_Comm intent = in [qc] standard_name = cloud_liquid_water_mixing_ratio diff --git a/physics/smoke_dust/rrfs_smoke_wrapper.F90 b/physics/smoke_dust/rrfs_smoke_wrapper.F90 index 3842cba54..0ac285b23 100755 --- a/physics/smoke_dust/rrfs_smoke_wrapper.F90 +++ b/physics/smoke_dust/rrfs_smoke_wrapper.F90 @@ -4,6 +4,7 @@ module rrfs_smoke_wrapper + use mpi_f90 use machine , only : kind_phys use rrfs_smoke_config, only : kemit, dust_opt, seas_opt, do_plumerise, & addsmoke_flag, plumerisefire_frq, wetdep_ls_opt, & @@ -225,7 +226,7 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, integer :: i, j, k, kp, n ! MPI variables integer :: mpiid - integer, intent(in) :: mpicomm + type(MPI_comm), intent(in) :: mpicomm integer, intent(in) :: mpirank integer, intent(in) :: mpiroot diff --git a/physics/smoke_dust/rrfs_smoke_wrapper.meta b/physics/smoke_dust/rrfs_smoke_wrapper.meta index 271d2dd36..1f26e5a93 100755 --- a/physics/smoke_dust/rrfs_smoke_wrapper.meta +++ b/physics/smoke_dust/rrfs_smoke_wrapper.meta @@ -774,7 +774,7 @@ long_name = MPI communicator units = index dimensions = () - type = integer + type = MPI_Comm intent = in [mpirank] standard_name = mpi_rank From 8c68c2f9b63ae1456678cf31a74fe0005a8bcb7c Mon Sep 17 00:00:00 2001 From: Dusan Jovic Date: Wed, 7 Feb 2024 14:52:01 -0600 Subject: [PATCH 09/21] Fix mpi use statement --- physics/smoke_dust/rrfs_smoke_wrapper.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/smoke_dust/rrfs_smoke_wrapper.F90 b/physics/smoke_dust/rrfs_smoke_wrapper.F90 index 0ac285b23..37ffccb35 100755 --- a/physics/smoke_dust/rrfs_smoke_wrapper.F90 +++ b/physics/smoke_dust/rrfs_smoke_wrapper.F90 @@ -4,7 +4,7 @@ module rrfs_smoke_wrapper - use mpi_f90 + use mpi_f08 use machine , only : kind_phys use rrfs_smoke_config, only : kemit, dust_opt, seas_opt, do_plumerise, & addsmoke_flag, plumerisefire_frq, wetdep_ls_opt, & From bb9b90fefdadb82a97d68e9fee2c6d7326e6020f Mon Sep 17 00:00:00 2001 From: Dusan Jovic Date: Mon, 26 Feb 2024 13:44:20 +0000 Subject: [PATCH 10/21] Update MPI find_package and check if MPI F08 module is supported --- CMakeLists.txt | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index d0d265437..b8cf73f7c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required(VERSION 3.3) +cmake_minimum_required(VERSION 3.10) project(ccpp_physics VERSION 5.0.0 @@ -9,9 +9,10 @@ set(PACKAGE "ccpp-physics") set(AUTHORS "Grant Firl" "Dustin Swales" "Man Zhang" "Mike Kavulich" ) #------------------------------------------------------------------------------ -# Set MPI flags for C/C++/Fortran -if (MPI) - find_package(MPI REQUIRED C Fortran) +# Set MPI flags for C/C++/Fortran with MPI F08 interface +find_package(MPI REQUIRED C Fortran) +if(NOT MPI_Fortran_HAVE_F08_MODULE) + message(FATAL_ERROR "MPI implementation does not support the Fortran 2008 mpi_f08 interface") endif() #------------------------------------------------------------------------------ From 0cb5d97355d347f5448568bba1b64409f6e1ddf5 Mon Sep 17 00:00:00 2001 From: Dusan Jovic Date: Tue, 27 Feb 2024 00:01:41 +0000 Subject: [PATCH 11/21] Remove C from MPI find_package --- CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index b8cf73f7c..c56070123 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,8 +9,8 @@ set(PACKAGE "ccpp-physics") set(AUTHORS "Grant Firl" "Dustin Swales" "Man Zhang" "Mike Kavulich" ) #------------------------------------------------------------------------------ -# Set MPI flags for C/C++/Fortran with MPI F08 interface -find_package(MPI REQUIRED C Fortran) +# Set MPI flags for Fortran with MPI F08 interface +find_package(MPI REQUIRED Fortran) if(NOT MPI_Fortran_HAVE_F08_MODULE) message(FATAL_ERROR "MPI implementation does not support the Fortran 2008 mpi_f08 interface") endif() From df5d9d4366cfd83f1af09cb3af3cfa10ffb31226 Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Thu, 7 Mar 2024 12:47:52 -0500 Subject: [PATCH 12/21] add cs_conv V2 changes on top of latest ufs/dev code --- physics/CONV/Chikira_Sugiyama/cs_conv.F90 | 1468 ++++++++++++--------- 1 file changed, 820 insertions(+), 648 deletions(-) diff --git a/physics/CONV/Chikira_Sugiyama/cs_conv.F90 b/physics/CONV/Chikira_Sugiyama/cs_conv.F90 index ab7388df8..94dadba87 100644 --- a/physics/CONV/Chikira_Sugiyama/cs_conv.F90 +++ b/physics/CONV/Chikira_Sugiyama/cs_conv.F90 @@ -60,6 +60,7 @@ module cs_conv !DD and precipitation. Decrease for more precip real(kind_phys), public :: precz0, preczh, clmd, clmp, clmdpa + real(kind_phys), public, parameter :: c0t=0.002, d0t=0.002 ! ! Private data ! @@ -225,15 +226,17 @@ subroutine cs_conv_run( IJSDIM , KMAX , ntracp1 , NN, & ! ! output arguments of CS_CUMLUS ! - real(kind_phys), dimension(IJSDIM,KMAX,nctp) :: vverti + real(kind_phys), dimension(IJSDIM,KMAX+1,nctp) :: vverti, sigmai real(kind_phys) GTT(IJSDIM,KMAX) !< temperature tendency [K/s] real(kind_phys) GTQ(IJSDIM,KMAX,NTR) !< tracer tendency [kg/kg/s] real(kind_phys) GTU(IJSDIM,KMAX) !< zonal velocity tendency [m/s2] real(kind_phys) GTV(IJSDIM,KMAX) !< meridional velocity tendency [m/s2] - real(kind_phys) GTPRP(IJSDIM,KMAX) !< precipitation (including snowfall) flux at interfaces [kg/m2/s] - real(kind_phys) GSNWP(IJSDIM,KMAX) !< snowfall flux at interfaces [kg/m2/s] - + real(kind_phys) CMDET(IJSDIM,KMAX) !< detrainment mass flux [kg/m2/s] + real(kind_phys) GTPRP(IJSDIM,KMAX+1) !< precipitation (including snowfall) flux at interfaces [kg/m2/s] + real(kind_phys) GSNWP(IJSDIM,KMAX+1) !< snowfall flux at interfaces [kg/m2/s] + real(kind_phys) GMFX0(IJSDIM,KMAX+1) !< updraft mass flux [kg/m2/s] + real(kind_phys) GMFX1(IJSDIM,KMAX+1) !< downdraft mass flux [kg/m2/s] integer KT(IJSDIM,nctp) !< cloud top index for each cloud type real(kind_phys) :: cape(IJSDIM) !< convective available potential energy (J/kg) @@ -377,13 +380,14 @@ subroutine cs_conv_run( IJSDIM , KMAX , ntracp1 , NN, & !> -# Initialize the sigma diagnostics do n=1,nctp - do k=1,kmax + do k=1,kmax+1 do i=ists,iens vverti(i,k,n) = zero + sigmai(i,k,n) = zero enddo enddo enddo - do k=1,kmax + do k=1,kmax+1 do i=ists,iens sigma(i,k) = zero enddo @@ -394,9 +398,9 @@ subroutine cs_conv_run( IJSDIM , KMAX , ntracp1 , NN, & otspt(1:ntr,1), otspt(1:ntr,2), & lprnt , ipr , & GTT , GTQ , GTU , GTV , & ! output - dt_mf , & ! output - GTPRP , GSNWP , ud_mf , & ! output - dd_mf , cape , KT , & ! output + CMDET , & ! output + GTPRP , GSNWP , GMFX0 , & ! output + GMFX1 , cape , KT , & ! output CBMFX , & ! modified GDT , GDQ , GDU , GDV , & ! input GDTM , & ! input @@ -404,7 +408,7 @@ subroutine cs_conv_run( IJSDIM , KMAX , ntracp1 , NN, & delp , delpi , & DELTA , DELTI , ISTS , IENS, mype,& ! input fscav, fswtr, wcbmaxm, nctp, & - sigma, vverti, & ! input/output !DDsigma + sigmai, sigma, vverti, & ! input/output !DDsigma do_aw, do_awdd, flx_form) ! ! @@ -432,6 +436,10 @@ subroutine cs_conv_run( IJSDIM , KMAX , ntracp1 , NN, & t(i,k) = GDT(i,k) + GTT(i,k) * delta u(i,k) = GDU(i,k) + GTU(i,k) * delta v(i,k) = GDV(i,k) + GTV(i,k) * delta +! Set the mass fluxes. + ud_mf (i,k) = GMFX0(i,k) + dd_mf (i,k) = GMFX1(i,k) + dt_mf (i,k) = CMDET(i,k) enddo enddo @@ -458,8 +466,8 @@ subroutine cs_conv_run( IJSDIM , KMAX , ntracp1 , NN, & ! CNV_PRC3(i,k) = 0.0 CNV_NDROP(i,k) = 0.0 CNV_NICE(i,k) = 0.0 - cf_upi(i,k) = max(0.0, min(1.0, 0.5*(sigma(i,k)+sigma(i,kp1)))) - CLCN(i,k) = cf_upi(i,k) !downdraft is below updraft +! cf_upi(i,k) = max(0.0,min(0.01*log(1.0+500*ud_mf(i,k)),0.1)) +! CLCN(i,k) = cf_upi(i,k) !downdraft is below updraft !! clcn(i,k) = max(0.0,min(0.01*log(1.0+500*ud_mf(i,k)/delta),0.25)) w_upi(i,k) = 0.0 @@ -492,9 +500,9 @@ subroutine cs_conv_run( IJSDIM , KMAX , ntracp1 , NN, & ! CNV_PRC3(i,k) = 0.0 CNV_NDROP(i,k) = 0.0 CNV_NICE(i,k) = 0.0 - cf_upi(i,k) = max(0.0,min(0.01*log(1.0+500*ud_mf(i,k)),0.25)) +! cf_upi(i,k) = max(0.0,min(0.01*log(1.0+500*ud_mf(i,k)),0.1)) ! & 500*ud_mf(i,k)),0.60)) - CLCN(i,k) = cf_upi(i,k) !downdraft is below updraft +! CLCN(i,k) = cf_upi(i,k) !downdraft is below updraft w_upi(i,k) = ud_mf(i,k)*(t(i,k)+epsvt*gdq(i,k,1)) * rair & / (max(cf_upi(i,k),1.e-12)*gdp(i,k)) @@ -586,11 +594,11 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions GDT , GDQ , GDU , GDV , & ! input GDTM , & ! input GDP , GDPM , GDZ , GDZM , & ! input - delp , delpi , & + delp , delpinv , & DELTA , DELTI , ISTS , IENS, mype,& ! input fscav, fswtr, wcbmaxm, nctp, & ! - sigma, vverti, & ! input/output !DDsigma - do_aw, do_awdd, flx_form ) + sigmai, sigma, vverti, & ! input/output !DDsigma + do_aw, do_awdd, flx_form) ! IMPLICIT NONE @@ -598,6 +606,8 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions INTEGER, INTENT(IN) :: im, IJSDIM, KMAX, NTR, mype, nctp, ipr !! DD, for GFS, pass in logical, intent(in) :: do_aw, do_awdd, flx_form ! switch to apply Arakawa-Wu to the tendencies logical, intent(in) :: otspt1(ntr), otspt2(ntr), lprnt + REAL(kind_phys),intent(in) :: DELP (IJSDIM, KMAX) + REAL(kind_phys),intent(in) :: DELPINV (IJSDIM, KMAX) ! ! [OUTPUT] REAL(kind_phys), INTENT(OUT) :: GTT (IJSDIM, KMAX ) ! heating rate @@ -605,35 +615,35 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions REAL(kind_phys), INTENT(OUT) :: GTU (IJSDIM, KMAX ) ! tendency of u REAL(kind_phys), INTENT(OUT) :: GTV (IJSDIM, KMAX ) ! tendency of v REAL(kind_phys), INTENT(OUT) :: CMDET (IJSDIM, KMAX ) ! detrainment mass flux - + REAL(kind_phys) :: GTLDET( IJSDIM, KMAX ) ! cloud liquid tendency by detrainment + REAL(kind_phys) :: GTIDET( IJSDIM, KMAX ) ! cloud ice tendency by detrainment ! assuming there is no flux at the top of the atmospherea - Moorthi - REAL(kind_phys), INTENT(OUT) :: GTPRP (IJSDIM, KMAX ) ! rain+snow flux - REAL(kind_phys), INTENT(OUT) :: GSNWP (IJSDIM, KMAX ) ! snowfall flux - REAL(kind_phys), INTENT(OUT) :: GMFX0 (IJSDIM, KMAX ) ! updraft mass flux - REAL(kind_phys), INTENT(OUT) :: GMFX1 (IJSDIM, KMAX ) ! downdraft mass flux + REAL(kind_phys), INTENT(OUT) :: GTPRP (IJSDIM, KMAX+1 ) ! rain+snow flux + REAL(kind_phys), INTENT(OUT) :: GSNWP (IJSDIM, KMAX+1 ) ! snowfall flux + REAL(kind_phys), INTENT(OUT) :: GMFX0 (IJSDIM, KMAX+1 ) ! updraft mass flux + REAL(kind_phys), INTENT(OUT) :: GMFX1 (IJSDIM, KMAX+1 ) ! downdraft mass flux REAL(kind_phys), INTENT(OUT) :: CAPE (IJSDIM ) INTEGER , INTENT(OUT) :: KT (IJSDIM, NCTP ) ! cloud top ! ! [MODIFIED] - REAL(kind_phys), INTENT(INOUT) :: CBMFX (IM, NCTP) ! cloud base mass flux - -!DDsigma - output added for AW sigma diagnostics -! sigma and vert. velocity as a function of cloud type (1==sfc) - real(kind_phys), intent(out), dimension(IM,KMAX) :: sigma !sigma totaled over cloud type - on interfaces (1=sfc) - real(kind_phys), intent(out), dimension(IM,KMAX,nctp) :: vverti + REAL(kind_phys), INTENT(INOUT) :: CBMFX ( IM, NCTP ) !! cloud base mass flux + !DDsigma - output added for AW sigma diagnostics + real(kind_phys), intent(out) :: sigmai(IM,KMAX+1,nctp) !DDsigma sigma by cloud type - on interfaces (1=sfc) + real(kind_phys), intent(out) :: vverti(IM,KMAX+1,nctp) !DDsigma vert. vel. by cloud type - on interfaces (1=sfc) + real(kind_phys), intent(out) :: sigma(IM,KMAX+1) !DDsigma sigma totaled over cloud type - on interfaces (1=sfc) + ! for computing AW flux form of tendencies -! The tendencies are summed over all cloud types -! real(kind_phys), intent(out), dimension(IM,KMAX) :: & !DDsigmadiag - real(kind_phys), allocatable, dimension(:,:) :: sfluxterm, qvfluxterm,& ! tendencies of DSE and water vapor due to eddy mass flux - qlfluxterm, qifluxterm,& ! tendencies of cloud water and cloud ice due to eddy mass flux +! real(kind_phys), dimension(IM,KMAX) :: & !DDsigmadiag +! sfluxterm, qvfluxterm +! real(kind_phys), dimension(IM,KMAX) :: & !DDsigmadiag +! qlfluxterm, qifluxterm +! real(kind_phys), dimension(ijsdim,kmax,ntrq:ntr) :: trfluxterm ! tendencies of tracers due to eddy mass flux + real(kind_phys), dimension(IM,KMAX) :: & !DDsigmadiag + condtermt, condtermq, frzterm, prectermq, prectermfrz + !DDsigma -! The fluxes are for an individual cloud type and reused. -! condtermt, condtermq are eddy flux of temperature and water vapor - condtermt, condtermq, frzterm, & - prectermq, prectermfrz - real(kind_phys), allocatable, dimension(:,:,:) :: trfluxterm ! tendencies of tracers due to eddy mass flux ! ! [INPUT] REAL(kind_phys), INTENT(IN) :: GDT (IJSDIM, KMAX ) ! temperature T @@ -653,167 +663,190 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions ! ! [INTERNAL WORK] REAL(kind_phys), allocatable :: GPRCC (:, :) ! rainfall -! REAL(kind_phys) GPRCC (IJSDIM, NTR) ! rainfall -! REAL(kind_phys) GSNWC (IJSDIM) ! snowfall -! REAL(kind_phys) CUMCLW(IJSDIM, KMAX) ! cloud water in cumulus -! REAL(kind_phys) CUMFRC(IJSDIM) ! cumulus cloud fraction -! -! REAL(kind_phys) GTCFRC(IJSDIM, KMAX) ! change in cloud fraction -! REAL(kind_phys) FLIQC (IJSDIM, KMAX) ! liquid ratio in cumulus -! -! REAL(kind_phys) GDCFRC(IJSDIM, KMAX) ! cloud fraction -! - REAL(kind_phys) GDW (IJSDIM, KMAX) ! total water - REAL(kind_phys) DELP (IJSDIM, KMAX) - REAL(kind_phys) DELPI (IJSDIM, KMAX) - REAL(kind_phys) GDQS (IJSDIM, KMAX) ! saturate moisture - REAL(kind_phys) FDQS (IJSDIM, KMAX) - REAL(kind_phys) GAM (IJSDIM, KMAX) - REAL(kind_phys) GDS (IJSDIM, KMAX) ! dry static energy - REAL(kind_phys) GDH (IJSDIM, KMAX) ! moist static energy - REAL(kind_phys) GDHS (IJSDIM, KMAX) ! saturate MSE -! - REAL(kind_phys) GCYM (IJSDIM, KMAX, NCTP)! norm. mass flux (half lev) - REAL(kind_phys) GCHB (IJSDIM) ! cloud base MSE-Li*Qi - REAL(kind_phys) GCWB (IJSDIM) ! cloud base total water - REAL(kind_phys) GCUB (IJSDIM) ! cloud base U - REAL(kind_phys) GCVB (IJSDIM) ! cloud base V - REAL(kind_phys) GCIB (IJSDIM) ! cloud base ice - REAL(kind_phys) GCtrB (IJSDIM,ntrq:ntr) ! cloud base tracer - REAL(kind_phys) GCYT (IJSDIM, NCTP) ! norm. mass flux @top - REAL(kind_phys) GCHT (IJSDIM, NCTP) ! cloud top MSE - REAL(kind_phys) GCQT (IJSDIM, NCTP) ! cloud top q - REAL(kind_phys) GCwT (IJSDIM) ! cloud top total water - REAL(kind_phys) GCUT (IJSDIM, NCTP) ! cloud top U - REAL(kind_phys) GCVT (IJSDIM, NCTP) ! cloud top V - REAL(kind_phys) GCLT (IJSDIM, NCTP) ! cloud top cloud water - REAL(kind_phys) GCIT (IJSDIM, NCTP) ! cloud top cloud ice + REAL(kind_phys) GSNWC ( IJSDIM ) !! snowfall + REAL(kind_phys) CUMCLW( IJSDIM, KMAX ) !! cloud water in cumulus + REAL(kind_phys) CUMFRC( IJSDIM ) !! cumulus cloud fraction +!COSP + REAL(kind_phys) QLIQC ( IJSDIM, KMAX ) !! cumulus cloud liquid water [kg/kg] + REAL(kind_phys) QICEC ( IJSDIM, KMAX ) !! cumulus cloud ice [kg/kg] + REAL(kind_phys) GPRCPF( IJSDIM, KMAX ) !! rainfall flux at full level + REAL(kind_phys) GSNWPF( IJSDIM, KMAX ) !! snowfall flux at full level +! + REAL(kind_phys) GTCFRC( IJSDIM, KMAX ) !! change in cloud fraction + REAL(kind_phys) FLIQC ( IJSDIM, KMAX ) !! liquid ratio in cumulus +! +!#ifdef OPT_CHASER +! REAL(kind_phys) RFXC ( IJSDIM, KMAX+1 ) !! precipi. flx [kg/m2/s] +! REAL(kind_phys) SFXC ( IJSDIM, KMAX+1 ) !! ice/snow flx [kg/m2/s] +! INTEGER LEVCUM( IJSDIM, KMAX ) !! flag for cum. cloud top +! REAL(kind_phys) LNFRC ( IJSDIM, KMAX ) !! areal rates of clouds +! REAL(kind_phys) REVC ( IJSDIM, KMAX ) !! evaporation rates +!#endif +! + REAL(kind_phys) GDCFRC( IJSDIM, KMAX ) !! cloud fraction +! +! REAL(kind_phys) GTQL ( IJSDIM, KMAX ) !! tendency of cloud liquid +! + REAL(kind_phys) GDW ( IJSDIM, KMAX ) !! total water + REAL(kind_phys) GDQS ( IJSDIM, KMAX ) !! saturate moisture + REAL(kind_phys) FDQS ( IJSDIM, KMAX ) + REAL(kind_phys) GAM ( IJSDIM, KMAX ) + REAL(kind_phys) GDS ( IJSDIM, KMAX ) !! dry static energy + REAL(kind_phys) GDH ( IJSDIM, KMAX ) !! moist static energy + REAL(kind_phys) GDHS ( IJSDIM, KMAX ) !! saturate MSE +! + REAL(kind_phys) GCYM ( IJSDIM, KMAX, NCTP ) !! norm. mass flux (half lev) + REAL(kind_phys) GCHB ( IJSDIM ) !! cloud base MSE-Li*Qi + REAL(kind_phys) GCWB ( IJSDIM ) !! cloud base total water + REAL(kind_phys) GCtrB ( IJSDIM, ntrq:ntr ) !! cloud base water vapor tracer + REAL(kind_phys) GCUB ( IJSDIM ) !! cloud base U + REAL(kind_phys) GCVB ( IJSDIM ) !! cloud base V + REAL(kind_phys) GCIB ( IJSDIM ) !! cloud base ice + REAL(kind_phys) ELAM ( IJSDIM, KMAX, NCTP ) !! entrainment (rate*massflux) + REAL(kind_phys) GCYT ( IJSDIM, NCTP ) !! norm. mass flux @top + REAL(kind_phys) GCHT ( IJSDIM, NCTP ) !! cloud top MSE + REAL(kind_phys) GCQT ( IJSDIM, NCTP ) !! cloud top q + REAL(kind_phys) GCwT ( IJSDIM ) !! cloud top total water + REAL(kind_phys) GCUT ( IJSDIM, NCTP ) !! cloud top U + REAL(kind_phys) GCVT ( IJSDIM, NCTP ) !! cloud top V + REAL(kind_phys) GCLT ( IJSDIM, NCTP ) !! cloud top cloud water + REAL(kind_phys) GCIT ( IJSDIM, NCTP ) !! cloud top cloud ice REAL(kind_phys) GCtrT (IJSDIM, ntrq:ntr, NCTP) ! cloud top tracer - REAL(kind_phys) GTPRT (IJSDIM, NCTP) ! precipitation/M - REAL(kind_phys) GCLZ (IJSDIM, KMAX) ! cloud liquid for each CTP - REAL(kind_phys) GCIZ (IJSDIM, KMAX) ! cloud ice for each CTP - -! REAL(kind_phys) ACWF (IJSDIM, NCTP) ! cloud work function - REAL(kind_phys) ACWF (IJSDIM ) ! cloud work function - REAL(kind_phys) GPRCIZ(IJSDIM, KMAX) ! precipitation - REAL(kind_phys) GSNWIZ(IJSDIM, KMAX) ! snowfall - REAL(kind_phys) GTPRC0(IJSDIM) ! precip. before evap. - - REAL(kind_phys) GMFLX (IJSDIM, KMAX) ! mass flux (updraft+downdraft) - REAL(kind_phys) QLIQ (IJSDIM, KMAX) ! total cloud liquid - REAL(kind_phys) QICE (IJSDIM, KMAX) ! total cloud ice - REAL(kind_phys) GPRCI (IJSDIM, KMAX) ! rainfall generation - REAL(kind_phys) GSNWI (IJSDIM, KMAX) ! snowfall generation - - REAL(kind_phys) GPRCP (IJSDIM, KMAX) ! rainfall flux -! - REAL(kind_phys) GTEVP (IJSDIM, KMAX) ! evaporation+sublimation - REAL(kind_phys) GMDD (IJSDIM, KMAX) ! downdraft mass flux - -! REAL(kind_phys) CUMHGT(IJSDIM, NCTP) ! cloud top height -! REAL(kind_phys) CTOPP (IJSDIM) ! cloud top pressure - - REAL(kind_phys) GDZTR (IJSDIM) ! tropopause height -! REAL(kind_phys) FLIQOU(IJSDIM, KMAX) ! liquid ratio in cumulus - INTEGER KB (IJSDIM) - INTEGER KSTRT (IJSDIM) ! tropopause level - REAL(kind_phys) GAMX - REAL(kind_phys) CIN (IJSDIM) - INTEGER JBUOY (IJSDIM) - REAL(kind_phys) DELZ, BUOY, DELWC, DELER - REAL(kind_phys) WCBX (IJSDIM) -! REAL(kind_phys) ERMR (NCTP) ! entrainment rate (ASMODE) -! SAVE ERMR - INTEGER KTMX (NCTP) ! max of cloud top - INTEGER KTMXT ! max of cloud top -! REAL(kind_phys) TIMED - REAL(kind_phys) GDCLDX, GDMU2X, GDMU3X -! -! REAL(kind_phys) HBGT (IJSDIM) ! imbalance in column heat -! REAL(kind_phys) WBGT (IJSDIM) ! imbalance in column water + REAL(kind_phys) GTPRT ( IJSDIM, NCTP ) !! precipitation/M + REAL(kind_phys) GCLZ ( IJSDIM, KMAX ) !! cloud liquid for each CTP + REAL(kind_phys) GCIZ ( IJSDIM, KMAX ) !! cloud ice for each CTP + + REAL(kind_phys) ACWF ( IJSDIM ) !! cloud work function + REAL(kind_phys) GPRCIZ( IJSDIM, KMAX+1, NCTP ) !! precipitation + REAL(kind_phys) GSNWIZ( IJSDIM, KMAX+1, NCTP ) !! snowfall + REAL(kind_phys) GTPRC0( IJSDIM ) !! precip. before evap. + + REAL(kind_phys) GMFLX ( IJSDIM, KMAX+1 ) !! mass flux (updraft+downdraft) + REAL(kind_phys) QLIQ ( IJSDIM, KMAX ) !! total cloud liquid + REAL(kind_phys) QICE ( IJSDIM, KMAX ) !! total cloud ice + REAL(kind_phys) GPRCI ( IJSDIM, KMAX ) !! rainfall generation + REAL(kind_phys) GSNWI ( IJSDIM, KMAX ) !! snowfall generation + + REAL(kind_phys) GPRCP ( IJSDIM, KMAX+1 ) !! rainfall flux +! + REAL(kind_phys) GTEVP ( IJSDIM, KMAX ) !! evaporation+sublimation + REAL(kind_phys) GMDD ( IJSDIM, KMAX+1 ) !! downdraft mass flux + + REAL(kind_phys) CUMHGT( IJSDIM, NCTP ) !! cloud top height + REAL(kind_phys) CTOPP ( IJSDIM ) !! cloud top pressure + + REAL(kind_phys) GDZTR ( IJSDIM ) !! tropopause height + REAL(kind_phys) FLIQOU( IJSDIM, KMAX ) !! liquid ratio in cumulus +!#ifdef OPT_CHASER +! REAL(kind_phys) TOPFLX( IJSDIM, NCTP ) !! flux at each cloud top +!#endif + INTEGER KB ( IJSDIM ) + INTEGER KSTRT ( IJSDIM ) !! tropopause level + REAL(kind_phys) GAMX + REAL(kind_phys) CIN ( IJSDIM ) + INTEGER JBUOY ( IJSDIM ) + REAL(kind_phys) DELZ, BUOY, DELWC, DELER +!M REAL(kind_phys) WCB ( NCTP ) !! updraft velocity**2 @base +!M SAVE WCB + REAL(kind_phys) WCBX (IJSDIM) +! REAL(kind_phys) ERMR ( NCTP ) !! entrainment rate (ASMODE) +! SAVE ERMR + INTEGER KTMX ( NCTP ) !! max of cloud top + INTEGER KTMXT !! max of cloud top + REAL(kind_phys) TIMED + REAL(kind_phys) GDCLDX, GDMU2X, GDMU3X +! + LOGICAL OOUT1, OOUT2 + INTEGER KBMX, I, K, CTP, ierr, n, kp1, l, l1, kk, kbi, kmi, km1 + real(kind_phys) tem1, tem2, tem3, cbmfl, mflx_e, teme, tems + + REAL(kind_phys) HBGT ( IJSDIM ) !! imbalance in column heat + REAL(kind_phys) WBGT ( IJSDIM ) !! imbalance in column water -!DDsigma begin local work variables - all on model interfaces (sfc=1) - REAL(kind_phys) lamdai ! lamda for cloud type ctp - REAL(kind_phys) gdqm, gdlm, gdim ! water vapor + !DDsigma begin local work variables - all on model interfaces (sfc=1) + REAL(kind_phys) lamdai( IJSDIM, KMAX+1, nctp ) !! lamda for cloud type ctp + REAL(kind_phys) lamdaprod( IJSDIM, KMAX+1 ) !! product of (1+lamda) through cloud type ctp + REAL(kind_phys) gdrhom !! density + REAL(kind_phys) gdtvm !! virtual temperature + REAL(kind_phys) gdqm, gdwm,gdlm, gdim !! water vaper REAL(kind_phys) gdtrm(ntrq:ntr) ! tracer - -! the following are new arguments to cumup to get them out for AW - REAL(kind_phys) wcv (IJSDIM, KMAX) ! in-cloud vertical velocity - REAL(kind_phys) GCTM (IJSDIM, KMAX) ! cloud T (half lev) !DDsigmadiag make output - REAL(kind_phys) GCQM (IJSDIM, KMAX) ! cloud q (half lev) !DDsigmadiag make output - REAL(kind_phys) GCwM (IJSDIM, KMAX) ! cloud q (half lev) !DDsigmadiag make output - REAL(kind_phys) GCiM (IJSDIM, KMAX) ! cloud q (half lev) !DDsigmadiag make output - REAL(kind_phys) GClM (IJSDIM, KMAX) ! cloud q (half lev) !DDsigmadiag make output - REAL(kind_phys) GChM (IJSDIM, KMAX) ! cloud q (half lev) !DDsigmadiag make output + character(len=4) :: cproc !DDsigmadiag + + ! the following are new arguments to cumup to get them out + REAL(kind_phys) wcv( IJSDIM, KMAX+1, nctp) !! in-cloud vertical velocity + REAL(kind_phys) GCTM ( IJSDIM, KMAX+1 ) !! cloud T (half lev) !DDsigmadiag make output + REAL(kind_phys) GCQM ( IJSDIM, KMAX+1, nctp ) !! cloud q (half lev) !DDsigmadiag make output + REAL(kind_phys) GCwM ( IJSDIM, KMAX+1, nctp ) !! cloud q (half lev) !DDsigmadiag make output + REAL(kind_phys) GCiM ( IJSDIM, KMAX+1 ) !! cloud q (half lev) !DDsigmadiag make output + REAL(kind_phys) GClM ( IJSDIM, KMAX+1 ) !! cloud q (half lev) !DDsigmadiag make output + REAL(kind_phys) GChM ( IJSDIM, KMAX+1, nctp ) !! cloud q (half lev) !DDsigmadiag make output REAL(kind_phys) GCtrM (IJSDIM, KMAX, ntrq:ntr) ! cloud tracer (half lev) !DDsigmadiag make output - -! eddy flux profiles for dse, water vapor, cloud water, cloud ice - REAL(kind_phys), dimension(Kmax+1) :: sfluxtem, qvfluxtem, qlfluxtem, qifluxtem - REAL(kind_phys), dimension(Kmax+1,ntrq:ntr) :: trfluxtem ! tracer - -! tendency profiles - condensation heating, condensation moistening, heating due to -! freezing, total precip production, frozen precip production - REAL(kind_phys), dimension(ijsdim,Kmax) :: dtcondtem, dqcondtem, dtfrztem, dqprectem,& ! Moorthi - dfrzprectem, lamdaprod !< product of (1+lamda) through cloud type ctp - REAL(kind_phys), dimension(ijsdim,Kmax) :: dtevap, dqevap, dtmelt, dtsubl - -! factor to modify precip rate to force conservation of water. With bug fixes it's -! not doing anything now. - REAL(kind_phys), dimension(ijsdim) :: moistening_aw - real(kind_phys), dimension(ijsdim,kmax) :: gctbl, gcqbl,gcwbl, gcqlbl, gcqibl, & !DDsigmadiag updraft profiles below cloud Base - sigmad ! downdraft area fraction + +! these are the fluxes at the interfaces - AW will operate on them + REAL(kind_phys), dimension(ijsdim,Kmax+1,nctp) :: sfluxtem, qvfluxtem, qlfluxtem, qifluxtem + REAL(kind_phys), dimension(ijsdim,Kmax+1,ntrq:ntr,nctp) :: trfluxtem ! tracer + + REAL(kind_phys), dimension(ijsdim,Kmax+1) :: dtcondtem, dqcondtem, dtfrztem, dqprectem,dfrzprectem + REAL(kind_phys), dimension(ijsdim,Kmax) :: dtevap, dqevap, dtmelt, dtsubl + REAL(kind_phys), dimension(ijsdim) :: moistening_aw + real(kind_phys) rhs_q, rhs_h, sftem, qftem, qlftem, qiftem + real(kind_phys), dimension(ijsdim,kmax+1) :: gctbl, gcqbl,gcwbl, gcqlbl, gcqibl !DDsigmadiag updraft profiles below cloud Base real(kind_phys), dimension(ijsdim,kmax,ntrq:ntr) :: gctrbl !DDsigmadiag tracer updraft profiles below cloud Base -! rhs_q, rhs_h are residuals of condensed water, MSE budgets to compute condensation, -! and heating due to freezing - real(kind_phys) :: rhs_q, rhs_h, fsigma, sigmai, delpinv -! real(kind_phys) :: rhs_q, rhs_h, sftem, qftem, qlftem, qiftem, & -! fsigma ! factor to reduce mass flux terms (1-sigma**2) for AW -!DDsigma end local work variables -! -! profiles of heating due to precip evaporation, melting and sublimation, and the -! evap, melting and sublimation rates. - - REAL(kind_phys), allocatable, dimension(:,:) :: dtdwn, & ! t tendency downdraft detrainment - dqvdwn, & ! qv tendency downdraft detrainment - dqldwn, & ! ql tendency downdraft detrainment - dqidwn ! qi tendency downdraft detrainment - REAL(kind_phys), allocatable, dimension(:,:,:) :: dtrdwn ! tracer tendency downdraft detrainment - + real(kind_phys), dimension(ijsdim,kmax+1) :: sigmad + real(kind_phys) :: fsigma( IJSDIM, KMAX+1 ) ! factor to reduce mass flux terms (1-sigma**2) for AW + real(kind_phys) :: lamdamax ! for sorting lamda values + integer loclamdamax + real(kind_phys) :: pr_tot, pr_ice, pr_liq !DDsigma end local work variables ! ! [INTERNAL PARM] - REAL(kind_phys), parameter :: WCBMIN = zero ! min. of updraft velocity at cloud base - -!M REAL(kind_phys) :: WCBMAX = 1.4_kind_phys ! max. of updraft velocity at cloud base + REAL(kind_phys) :: WCBMIN = 0._kind_phys !! min. of updraft velocity at cloud base +!M REAL(kind_phys) :: WCBMAX = 1.4_kind_phys !! max. of updraft velocity at cloud base !M wcbas commented by Moorthi since it is not used -!M REAL(kind_phys) :: WCBAS = 2._kind_phys ! updraft velocity**2 at cloud base (ASMODE) -!M REAL(kind_phys) :: ERAMIN = 1.e-5_kind_phys ! min. of entrainment rate - ! used only in OPT_ASMODE -!M REAL(kind_phys) :: ERAMAX = 2.e-3_kind_phys ! max. of entrainment rate - ! used only in OPT_ASMODE - LOGICAL :: OINICB = .false. ! set 0.d0 to CBMFX when .true. - -! REAL(kind_phys) :: VARMIN = 1.e-13_kind_phys ! minimum of PDF variance -! REAL(kind_phys) :: VARMAX = 5.e-7_kind_phys ! maximum of PDF variance -! REAL(kind_phys) :: SKWMAX = 0.566_kind_phys ! maximum of PDF skewness +!M REAL(kind_phys) :: WCBAS = 2._kind_phys !! updraft velocity**2 at cloud base (ASMODE) +!M REAL(kind_phys) :: ERAMIN = 1.e-5_kind_phys !! min. of entrainment rate + !! used only in OPT_ASMODE +!M REAL(kind_phys) :: ERAMAX = 2.e-3_kind_phys !! max. of entrainment rate + !! used only in OPT_ASMODE +! downdraft mass flux terms now slot nctp+1 in the *fluxterm arrays + REAL(kind_phys) dtdwn ( IJSDIM, KMAX ) !! t tendency downdraft detrainment + REAL(kind_phys) dqvdwn ( IJSDIM, KMAX ) !! qv tendency downdraft detrainment + REAL(kind_phys) dqldwn ( IJSDIM, KMAX ) !! ql tendency downdraft detrainment + REAL(kind_phys) dqidwn ( IJSDIM, KMAX ) !! qi tendency downdraft detrainment + REAL(kind_phys), dimension(ijsdim,kmax,ntrq:ntr) :: dtrdwn ! tracer tendency downdraft detrainment + + LOGICAL :: OINICB = .false. !! set 0.d0 to CBMFX + + REAL(kind_phys) :: VARMIN = 1.e-13_kind_phys !! minimum of PDF variance + REAL(kind_phys) :: VARMAX = 5.e-7_kind_phys !! maximum of PDF variance + REAL(kind_phys) :: SKWMAX = 0.566_kind_phys !! maximum of PDF skewness + + REAL(kind_phys) :: PSTRMX = 400.e2_kind_phys !! max P of tropopause + REAL(kind_phys) :: PSTRMN = 50.e2_kind_phys !! min P of tropopause + REAL(kind_phys) :: GCRSTR = 1.e-4_kind_phys !! crit. dT/dz tropopause - REAL(kind_phys) :: PSTRMX = 400.e2_kind_phys ! max P of tropopause - REAL(kind_phys) :: PSTRMN = 50.e2_kind_phys ! min P of tropopause - REAL(kind_phys) :: GCRSTR = 1.e-4_kind_phys ! crit. dT/dz tropopause - - real(kind=kind_phys) :: tem, esat, mflx_e, cbmfl, tem1, tem2, tem3 - INTEGER :: KBMX, I, K, CTP, ierr, n, kp1, km1, kk, kbi, l, l1 + ! 0: mass fixer is not applied + ! tracers which may become negative values + ! e.g. subgrid-PDFs + ! 1: mass fixer is applied, total mass may change through cumulus scheme + ! e.g. moisture, liquid cloud, ice cloud, aerosols + ! 2: mass fixer is applied, total mass never change through cumulus scheme + ! e.g. CO2 + real(kind=kind_phys), parameter :: zero=0.0, one=1.0 + real(kind=kind_phys) :: tem, esat ! - LOGICAL, SAVE :: OFIRST = .TRUE. ! called first time? + LOGICAL, SAVE :: OFIRST = .TRUE. !! called first time? ! + IF ( OFIRST ) THEN - IF (OFIRST) THEN OFIRST = .FALSE. IF (OINICB) THEN CBMFX = zero ENDIF ENDIF ! + + kp1 = kmax + 1 do n=1,ntr do k=1,kmax do i=1,ijsdim @@ -821,65 +854,82 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions enddo enddo enddo - + do k=1,kmax+1 + do i=1,ijsdim + gmflx(i,k) = zero + gmfx0(i,k) = zero + enddo + enddo do k=1,kmax do i=1,ijsdim - gtt(i,k) = zero - gtu(i,k) = zero - gtv(i,k) = zero - gmflx(i,k) = zero - gmfx0(i,k) = zero - gprci(i,k) = zero - gsnwi(i,k) = zero - qliq(i,k) = zero - qice(i,k) = zero -! gtcfrc(i,k) = zero -! cumclw(i,k) = zero -! fliqc(i,k) = zero - sigma(i,k) = zero + gtt(i,k) = zero + gtu(i,k) = zero + gtv(i,k) = zero + gprci(i,k) = zero + gsnwi(i,k) = zero + qliq(i,k) = zero + qice(i,k) = zero +! gtcfrc(i,k) = zero +! cumclw(i,k) = zero +! fliqc(i,k) = zero + fliqou(i,k) = zero + gprcpf(i,k) = zero + gsnwpf(i,k) = zero + cmdet(i,k) = zero enddo enddo if (flx_form) then - allocate(sfluxterm(ijsdim,kmax), qvfluxterm(ijsdim,kmax), qlfluxterm(ijsdim,kmax), & - qifluxterm(ijsdim,kmax), condtermt(ijsdim,kmax), condtermq(ijsdim,kmax), & - frzterm(ijsdim,kmax), prectermq(ijsdim,kmax), prectermfrz(ijsdim,kmax), & - dtdwn(ijsdim,kmax), dqvdwn(ijsdim,kmax), dqldwn(ijsdim,kmax), & - dqidwn(ijsdim,kmax), trfluxterm(ijsdim,kmax,ntrq:ntr), & - dtrdwn(ijsdim,kmax,ntrq:ntr)) - do k=1,kmax - do i=1,ijsdim - sfluxterm(i,k) = zero - qvfluxterm(i,k) = zero - qlfluxterm(i,k) = zero - qifluxterm(i,k) = zero - condtermt(i,k) = zero - condtermq(i,k) = zero - frzterm(i,k) = zero - prectermq(i,k) = zero - prectermfrz(i,k) = zero - dtdwn(i,k) = zero - dqvdwn(i,k) = zero - dqldwn(i,k) = zero - dqidwn(i,k) = zero - cmdet(i,k) = zero + do ctp = 1,nctp + do k=1,kp1 + do i=1,ijsdim + sfluxtem(i,k,ctp) = zero + qvfluxtem(i,k,ctp) = zero + qlfluxtem(i,k,ctp) = zero + qifluxtem(i,k,ctp) = zero + enddo + enddo + do n = ntrq,ntr + do k=1,kp1 + do i=1,ijsdim + trfluxtem(i,k,n,ctp) = zero + enddo + enddo enddo enddo - do n = ntrq,ntr do k=1,kmax do i=1,ijsdim - trfluxterm(i,k,n) = zero - dtrdwn(i,k,n) = zero + condtermt(i,k) = zero + condtermq(i,k) = zero + frzterm(i,k) = zero + prectermq(i,k) = zero + prectermfrz(i,k) = zero enddo enddo - enddo + do k=1,kmax + do i=1,ijsdim + dtdwn(i,k) = zero + dqvdwn(i,k) = zero + dqldwn(i,k) = zero + dqidwn(i,k) = zero + enddo + enddo + do n = ntrq,ntr + do k=1,kmax + do i=1,ijsdim + dtrdwn(i,k,n) = zero + enddo + enddo + enddo endif do i=1,ijsdim -! gprcc(i,:) = zero - gtprc0(i) = zero -! hbgt(i) = zero -! wbgt(i) = zero - gdztr(i) = zero - kstrt(i) = kmax +! gprcc(i,:) = zero +! gmflx(i,kp1) = zero + gmfx0(i,kp1) = zero + gtprc0(i) = zero +! hbgt(i) = zero +! wbgt(i) = zero + gdztr(i) = zero + kstrt(i) = kmax enddo do k=1,kmax @@ -907,9 +957,8 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions ! !> -# Compute tropopause height (GDZTR) DO K=1,KMAX - kp1 = k + 1 DO I=ISTS,IENS - GAMX = (GDTM(I,KP1)-GDTM(I,K)) / (GDZM(I,KP1)-GDZM(I,K)) + GAMX = (GDTM(I,K+1)-GDTM(I,K)) / (GDZM(I,K+1)-GDZM(I,K)) IF ((GDP(I,K) < PSTRMX .AND. GAMX > GCRSTR) .OR. GDP(I,K) < PSTRMN) THEN KSTRT(I) = MIN(K, KSTRT(I)) ENDIF @@ -925,12 +974,12 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions !> -# Call cumbas() to compute cloud base properties CALL CUMBAS(IJSDIM, KMAX , & !DD dimensions - KB , GCYM(1,1,1) , KBMX , & ! output + KB , GCYM(:,:,1) , KBMX , & ! output ntr , ntrq , & GCHB , GCWB , GCUB , GCVB , & ! output GCIB , gctrb, & ! output GDH , GDW , GDHS , GDQS , & ! input - GDQ(1,1,iti) , GDU , GDV , GDZM , & ! input + GDQ(:,:,iti) , GDU , GDV , GDZM , & ! input GDPM , FDQS , GAM , & ! input lprnt, ipr, & ISTS , IENS , & !) ! input @@ -955,7 +1004,7 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions CAPE(I) = CAPE(I) + BUOY * GRAV * (GDZM(I,K+1) - GDZM(I,K)) JBUOY(I) = 2 ELSEIF (BUOY < zero .AND. JBUOY(I) /= 2) THEN - CIN(I) = CIN(I) - BUOY * GRAV * (GDZM(I,K+1) - GDZM(I,K)) + CIN(I) = CIN(I) + BUOY * GRAV * (GDZM(I,K+1) - GDZM(I,K)) JBUOY(I) = -1 ENDIF endif @@ -968,12 +1017,25 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions !DDsigma some initialization before summing over cloud type !> -# Initialize variables before summing over cloud types - do k=1,kmax ! Moorthi + if(flx_form) then + do k=1,kp1 ! Moorthi do i=1,ijsdim lamdaprod(i,k) = one + sigma(i,k) = 0.0 enddo enddo + do ctp=1,nctp + do k=1,kp1 + do i=1,ijsdim + lamdai(i,k,ctp) = zero + sigmai(i,k,ctp) = zero + vverti(i,k,ctp) = zero + enddo + enddo + enddo + endif + do ctp=2,nctp do k=1,kmax do i=1,ijsdim @@ -990,15 +1052,6 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions WCBX(I) = DELWC * DELWC enddo - do k=1,kmax ! Moorthi - do i=1,ijsdim - dqcondtem(i,k) = zero - dqprectem(i,k) = zero - dfrzprectem(i,k) = zero - dtfrztem(i,k) = zero - dtcondtem(i,k) = zero - enddo - enddo ! getting more incloud profiles of variables to compute eddy flux tendencies ! and condensation rates @@ -1010,51 +1063,48 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions !> -# Call cumup() to compute in-cloud properties CALL CUMUP(IJSDIM, KMAX, NTR, ntrq, & !DD dimensions ACWF , & ! output - GCLZ , GCIZ , GPRCIZ , GSNWIZ, & ! output - GCYT(1,CTP) , GCHT(1,CTP) , GCQT (1,CTP), & ! output - GCLT(1,CTP) , GCIT(1,CTP) , GTPRT(1,CTP), & ! output - GCUT(1,CTP) , GCVT(1,CTP) , gctrt(1,ntrq:ntr,ctp), & ! output - KT (1,CTP) , KTMX(CTP) , & ! output - GCYM(1,1,CTP) , & ! modified - wcv , & ! !DD-sigma new output + GCLZ , GCIZ , GPRCIZ(:,:,CTP), GSNWIZ(:,:,CTP), & ! output + GCYT(:,CTP) , GCHT(:,CTP) , GCQT (:,CTP), & ! output + GCLT(:,CTP) , GCIT(:,CTP) , GTPRT(:,CTP), & ! output + GCUT(:,CTP) , GCVT(:,CTP) , gctrt(:,ntrq:ntr,ctp), & ! output + KT (:,CTP) , KTMX(CTP) , & ! output + GCYM(:,:,CTP) , & ! modified + wcv(:,:,CTP) , & ! !DD-sigma new output GCHB , GCWB , GCUB , GCVB , & ! input !DDsigmadiag GCIB , gctrb , & ! input GDU , GDV , GDH , GDW , & ! input GDHS , GDQS , GDT , GDTM , & ! input - GDQ , GDQ(1,1,iti) , GDZ , GDZM , & ! input + GDQ , GDQ(:,:,iti) , GDZ , GDZM , & ! input GDPM , FDQS , GAM , GDZTR , & ! input CPRES , WCBX , & ! input KB , CTP , ISTS , IENS , & ! input - gctm , gcqm, gcwm, gchm, gcwt, gclm, gcim, gctrm, & ! additional incloud profiles and cloud top total water + gctm , gcqm(:,:,CTP), gcwm(:,:,CTP), gchm(:,:,CTP),& + gcwt, gclm, gcim, gctrm, & ! additional incloud profiles and cloud top total water lprnt , ipr ) ! !> -# Call cumbmx() to compute cloud base mass flux CALL CUMBMX(IJSDIM, KMAX, & !DD dimensions - CBMFX(1,CTP), & ! modified - ACWF , GCYT(1,CTP), GDZM , & ! input + CBMFX(:,CTP), & ! modified + ACWF , GCYT(:,CTP), GDZM , & ! input GDW , GDQS , DELP , & ! input - KT (1,CTP), KTMX(CTP) , KB , & ! input + KT (:,CTP), KTMX(CTP) , KB , & ! input DELTI , ISTS , IENS ) !DDsigma - begin sigma computation ! At this point cbmfx is updated and we have everything we need to compute sigma - do i=ISTS,IENS - if (flx_form) then -!> -# Initialize eddy fluxes for cloud types - do k=1,kmax+1 - sfluxtem(k) = zero - qvfluxtem(k) = zero - qlfluxtem(k) = zero - qifluxtem(k) = zero - enddo - do n=ntrq,ntr ! tracers - do k=1,kmax+1 - trfluxtem(k,n) = zero - enddo + if (flx_form) then + do k=1,kmax + 1 ! Moorthi + do i=1,ijsdim + dqcondtem(i,k) = zero + dqprectem(i,k) = zero + dfrzprectem(i,k) = zero + dtfrztem(i,k) = zero + dtcondtem(i,k) = zero enddo - endif + enddo + do i=ISTS,IENS cbmfl = cbmfx(i,ctp) kk = kt(i,ctp) ! cloud top index @@ -1062,56 +1112,54 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions kbi = kb(i) ! cloud base index do k=kbi,kk ! loop from cloud base to cloud top km1 = k - 1 - rhs_h = zero - rhs_q = zero -!> -# Interpolate environment variables to layer interface +! get environment variables interpolated to layer interface GDQM = half * (GDQ(I,K,1) + GDQ(I,KM1,1)) ! as computed in cumup ! GDwM = half * (GDw(I,K) + GDw(I,KM1 )) - GDlM = half * (GDQ(I,K,3) + GDQ(I,KM1,3)) - GDiM = half * (GDQ(I,K,2) + GDQ(I,KM1,2)) + GDlM = half * (GDQ(I,K,itl) + GDQ(I,KM1,itl)) + GDiM = half * (GDQ(I,K,iti) + GDQ(I,KM1,iti)) do n = ntrq,NTR GDtrM(n) = half * (GDQ(I,K,n) + GDQ(I,KM1,n)) ! as computed in cumup enddo mflx_e = gcym(i,k,ctp) * cbmfl ! mass flux at level k for cloud ctp - if (do_aw) then !> -# Compute lamda for a cloud type and then updraft area fraction !! (sigmai) following Equations 23 and 12 of !! Arakawa and Wu (2013) \cite arakawa_and_wu_2013 , respectively - lamdai = mflx_e * rair * gdtm(i,k)*(one+epsvt*gdqm) & - / (gdpm(i,k)*wcv(i,k)) - lamdaprod(i,k) = lamdaprod(i,k) * (one+lamdai) - -! vverti(i,k,ctp) = wcv(i,k) -! sigmai(i,k,ctp) = lamdai / lamdaprod(i,k) -! sigma(i,k) = max(zero, min(one, sigma(i,k) + sigmai(i,k,ctp))) - - sigmai = lamdai / lamdaprod(i,k) - sigma(i,k) = max(zero, min(one, sigma(i,k) + sigmai)) - vverti(i,k,ctp) = sigmai * wcv(i,k) - else - sigma(i,k) = 0.0 - endif + lamdai(i,k,ctp) = mflx_e * rair * gdtm(i,k)*(one+epsvt*gdqm) & + / (gdpm(i,k)*wcv(i,k,ctp)) + +! just compute lamdai here, we will compute sigma, sigmai, and vverti outside +! the cloud type loop after we can sort lamdai +! lamdaprod(i,k) = lamdaprod(i,k) * (one+lamdai(i,k,ctp)) +! +!! vverti(i,k,ctp) = wcv(i,k) +!! sigmai(i,k,ctp) = lamdai / lamdaprod(i,k) +!! sigma(i,k) = max(zero, min(one, sigma(i,k) + sigmai(i,k,ctp))) +! +! sigmai(i,k,ctp) = lamdai(i,k,ctp) / lamdaprod(i,k) +! sigma(i,k) = max(zero, min(one, sigma(i,k) + sigmai(i,k,ctp))) +! vverti(i,k,ctp) = sigmai(i,k,ctp) * wcv(i,k,ctp) - if (flx_form) then +! sigma effect won't be applied until later, when lamda is sorted ! fsigma = 1.0 ! no aw effect, comment following lines to undo AW - fsigma = one - sigma(i,k) +! fsigma = one - sigma(i,k) !> -# Compute tendencies based on mass flux and condensation ! fsigma is the AW reduction of flux tendencies if(k == kbi) then do l=2,kbi ! compute eddy fluxes below cloud base - tem = - fsigma * gcym(i,l,ctp) * cbmfl +! tem = - fsigma * gcym(i,l,ctp) * cbmfl + tem = - gcym(i,l,ctp) * cbmfl ! first get environment variables at layer interface l1 = l - 1 GDQM = half * (GDQ(I,l,1) + GDQ(I,l1,1)) - GDlM = half * (GDQ(I,l,3) + GDQ(I,l1,3)) - GDiM = half * (GDQ(I,l,2) + GDQ(I,l1,2)) + GDlM = half * (GDQ(I,l,itl) + GDQ(I,l1,itl)) + GDiM = half * (GDQ(I,l,iti) + GDQ(I,l1,iti)) !! GDwM = half * (GDw(I,l) + GDw(I,l1)) do n = ntrq,NTR GDtrM(n) = half * (GDQ(I,l,n) + GDQ(I,l1,n)) ! as computed in cumup @@ -1119,12 +1167,12 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions ! flux = mass flux * (updraft variable minus environment variable) !centered differences - sfluxtem(l) = tem * (gdtm(i,l)-gctbl(i,l)) - qvfluxtem(l) = tem * (gdqm-gcqbl(i,l)) - qlfluxtem(l) = tem * (gdlm-gcqlbl(i,l)) - qifluxtem(l) = tem * (gdim-gcqibl(i,l)) + sfluxtem(i,l,ctp) = tem * (gdtm(i,l)-gctbl(i,l)) + qvfluxtem(i,l,ctp) = tem * (gdqm-gcqbl(i,l)) + qlfluxtem(i,l,ctp) = tem * (gdlm-gcqlbl(i,l)) + qifluxtem(i,l,ctp) = tem * (gdim-gcqibl(i,l)) do n = ntrq,NTR - trfluxtem(l,n) = tem * (gdtrm(n)-gctrbl(i,l,n)) + trfluxtem(i,l,n,ctp) = tem * (gdtrm(n)-gctrbl(i,l,n)) enddo ! if(lprnt .and. i == ipr) write(0,*)' l=',l,' kbi=',kbi,' tem =', tem,' trfluxtem=',trfluxtem(l,ntr),& ! ' gdtrm=',gdtrm(ntr),' gctrbl=',gctrbl(i,l,ntr),' gq=',GDQ(I,l,ntr),GDQ(I,l1,ntr),' l1=',l1,' ctp=',ctp,& @@ -1146,14 +1194,15 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions else ! flux = mass flux * (updraft variable minus environment variable) - tem = - fsigma * mflx_e +! tem = - fsigma * mflx_e + tem = - mflx_e !centered - sfluxtem(k) = tem * (gdtm(i,k)+gocp*gdzm(i,k)-gctm(i,k)) - qvfluxtem(k) = tem * (gdqm-gcqm(i,k)) - qlfluxtem(k) = tem * (gdlm-gclm(i,k)) - qifluxtem(k) = tem * (gdim-gcim(i,k)) + sfluxtem(i,k,ctp) = tem * (gdtm(i,k)+gocp*gdzm(i,k)-gctm(i,k)) + qvfluxtem(i,k,ctp) = tem * (gdqm-gcqm(i,k,ctp)) + qlfluxtem(i,k,ctp) = tem * (gdlm-gclm(i,k)) + qifluxtem(i,k,ctp) = tem * (gdim-gcim(i,k)) do n = ntrq,NTR - trfluxtem(k,n) = tem * (gdtrm(n)-gctrm(i,k,n)) + trfluxtem(i,k,n,ctp) = tem * (gdtrm(n)-gctrm(i,k,n)) enddo !upstream - This better matches what the original CS tendencies do @@ -1185,117 +1234,57 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions ! ' fsigma=',fsigma,' mflx_e=',mflx_e,' trfluxtemk=',trfluxtem(k,ntr),' sigma=',sigma(i,k) -! the condensation terms - these come from the MSE and condensed water budgets for -! an entraining updraft -! if(k > kb(i)) then ! comment for test -! if(k <= kk) then ! Moorthi -! if(k < kt(i,ctp)) then -! rhs_h = cbmfl*(gcym(i,k)*gchm(i,k) - (gcym(i,km1)*gchm(i,km1) & -! + GDH(I,Km1 )*(gcym(i,k)-gcym(i,km1))) ) -! rhs_q = cbmfl*(gcym(i,k)*(gcwm(i,k)-gcqm(i,k)) & -! - (gcym(i,km1)*(gcwm(i,km1)-gcqm(i,km1)) & -! + (GDw( I,Km1 )-gdq(i,km1,1))*(gcym(i,k)-gcym(i,km1))) ) -! tem = cbmfl * (one - sigma(i,k)) - tem = cbmfl * (one - 0.5*(sigma(i,k)+sigma(i,km1))) - tem1 = gcym(i,k,ctp) * (one - sigma(i,k)) - tem2 = gcym(i,km1,ctp) * (one - sigma(i,km1)) - rhs_h = cbmfl * (tem1*gchm(i,k) - (tem2*gchm(i,km1) & - + GDH(I,Km1)*(tem1-tem2)) ) - rhs_q = cbmfl * (tem1*(gcwm(i,k)-gcqm(i,k)) & - - (tem2*(gcwm(i,km1)-gcqm(i,km1)) & - + (GDw(I,Km1)-gdq(i,km1,1))*(tem1-tem2)) ) - -! ELSE -! rhs_h = cbmfl*(gcht(i,ctp) - (gcym(i,k-1)*gchm(i,k-1) + GDH( I,K-1 )*(gcyt(i,ctp)-gcym(i,k-1))) ) -! rhs_q = cbmfl*((gcwt(i)-gcqt(i,ctp)) - (gcym(i,k-1)*(gcwm(i,k-1)-gcqm(i,k-1)) + (GDw( I,K-1 )-gdq(i,k-1,1))*(gcyt(i,ctp)-gcym(i,k-1))) ) -! endif - -!> -# Compute condensation, total precipitation production, frozen precipitation production, -!! heating due to freezing, and total temperature tendency due to in-cloud microphysics - dqcondtem(i,km1) = -rhs_q ! condensation -! dqprectem(i,km1) = cbmfl * (GPRCIZ(i,k) + GSNWIZ(i,k)) - dqprectem(i,km1) = tem * (GPRCIZ(i,k) + GSNWIZ(i,k)) ! total precip production -! dfrzprectem(i,km1) = cbmfl * GSNWIZ(i,k) - dfrzprectem(i,km1) = tem * GSNWIZ(i,k) ! production of frozen precip - dtfrztem(i,km1) = rhs_h*oneocp ! heating due to freezing - dtcondtem(i,km1) = - elocp * dqcondtem(i,km1) + dtfrztem(i,km1) - endif ! if(k > kbi) then - endif ! if (flx_form) enddo ! end of k=kbi,kk loop endif ! end of if(cbmfl > zero) -! get tendencies by difference of fluxes, sum over cloud type - - if (flx_form) then - do k = 1,kk - delpinv = delpi(i,k) -!> -# Sum single cloud microphysical tendencies over all cloud types - condtermt(i,k) = condtermt(i,k) + dtcondtem(i,k) * delpinv - condtermq(i,k) = condtermq(i,k) + dqcondtem(i,k) * delpinv - prectermq(i,k) = prectermq(i,k) + dqprectem(i,k) * delpinv - prectermfrz(i,k) = prectermfrz(i,k) + dfrzprectem(i,k) * delpinv - frzterm(i,k) = frzterm(i,k) + dtfrztem(i,k) * delpinv - -!> -# Compute flux tendencies and vertical flux divergence - sfluxterm(i,k) = sfluxterm(i,k) - (sfluxtem(k+1) - sfluxtem(k)) * delpinv - qvfluxterm(i,k) = qvfluxterm(i,k) - (qvfluxtem(k+1) - qvfluxtem(k)) * delpinv - qlfluxterm(i,k) = qlfluxterm(i,k) - (qlfluxtem(k+1) - qlfluxtem(k)) * delpinv - qifluxterm(i,k) = qifluxterm(i,k) - (qifluxtem(k+1) - qifluxtem(k)) * delpinv - do n = ntrq,ntr - trfluxterm(i,k,n) = trfluxterm(i,k,n) - (trfluxtem(k+1,n) - trfluxtem(k,n)) * delpinv - enddo -! if (lprnt .and. i == ipr) write(0,*)' k=',k,' trfluxtem=',trfluxtem(k+1,ntr),trfluxtem(k,ntr),& -! ' ctp=',ctp,' trfluxterm=',trfluxterm(i,k,ntr) - enddo - endif ! if (flx_form) enddo ! end of i loop -! - do i=ists,iens - if (cbmfx(i,ctp) > zero) then - tem = one - sigma(i,kt(i,ctp)) - gcyt(i,ctp) = tem * gcyt(i,ctp) - gtprt(i,ctp) = tem * gtprt(i,ctp) - gclt(i,ctp) = tem * gclt(i,ctp) - gcht(i,ctp) = tem * gcht(i,ctp) - gcqt(i,ctp) = tem * gcqt(i,ctp) - gcit(i,ctp) = tem * gcit(i,ctp) - if (.not. flx_form) then - do n = ntrq,ntr - gctrt(i,n,ctp) = tem * gctrt(i,n,ctp) - enddo - end if - gcut(i,ctp) = tem * gcut(i,ctp) - gcvt(i,ctp) = tem * gcvt(i,ctp) - do k=1,kmax - kk = kb(i) - if (k < kk) then - tem = one - sigma(i,kk) - tem1 = tem - else - tem = one - sigma(i,k) - tem1 = one - 0.5*(sigma(i,k)+sigma(i,k-1)) - endif - gcym(i,k,ctp) = tem * gcym(i,k,ctp) - gprciz(i,k) = tem1 * gprciz(i,k) - gsnwiz(i,k) = tem1 * gsnwiz(i,k) - gclz(i,k) = tem1 * gclz(i,k) - gciz(i,k) = tem1 * gciz(i,k) - enddo - endif - enddo + endif ! if (flx_form) +! +! we don't reduce these values in AW, just the tendencies due to fluxes +! do i=ists,iens +! if (cbmfx(i,ctp) > zero) then +! tem = one - sigma(i,kt(i,ctp)) +! gcyt(i,ctp) = tem * gcyt(i,ctp) +! gtprt(i,ctp) = tem * gtprt(i,ctp) +! gclt(i,ctp) = tem * gclt(i,ctp) +! gcht(i,ctp) = tem * gcht(i,ctp) +! gcqt(i,ctp) = tem * gcqt(i,ctp) +! gcit(i,ctp) = tem * gcit(i,ctp) +! do n = ntrq,ntr +! gctrt(i,n,ctp) = tem * gctrt(i,n,ctp) +! enddo +! gcut(i,ctp) = tem * gcut(i,ctp) +! gcvt(i,ctp) = tem * gcvt(i,ctp) +! do k=1,kmax +! kk = kb(i) +! if (k < kk) then +! tem = one - sigma(i,kk) +! tem1 = tem +! else +! tem = one - sigma(i,k) +! tem1 = one - 0.5*(sigma(i,k)+sigma(i,k-1)) +! endif +! gcym(i,k,ctp) = tem * gcym(i,k,ctp) +! gprciz(i,k) = tem1 * gprciz(i,k) +! gsnwiz(i,k) = tem1 * gsnwiz(i,k) +! gclz(i,k) = tem1 * gclz(i,k) +! gciz(i,k) = tem1 * gciz(i,k) +! enddo +! endif +! enddo ! !> -# Call cumflx() to compute cloud mass flux and precipitation CALL CUMFLX(IM , IJSDIM, KMAX , & !DD dimensions GMFX0 , GPRCI , GSNWI , CMDET, & ! output QLIQ , QICE , GTPRC0, & ! output - CBMFX(1,CTP) , GCYM(1,1,ctp), GPRCIZ , GSNWIZ , & ! input - GTPRT(1,CTP) , GCLZ , GCIZ , GCYT(1,ctp),& ! input - KB , KT(1,CTP) , KTMX(CTP) , & ! input + CBMFX(:,CTP) , GCYM(:,:,ctp), GPRCIZ(:,:,CTP), GSNWIZ(:,:,CTP) , & ! input + GTPRT(:,CTP) , GCLZ , GCIZ , GCYT(:,ctp),& ! input + KB , KT(:,CTP) , KTMX(CTP) , & ! input ISTS , IENS ) ! input ENDDO ! end of cloud type ctp loop @@ -1333,46 +1322,127 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions GDH , GDQ , GDU , GDV , & ! input ! GTT , GTQ , GTCFRC, GTU , GTV , & ! modified ! GDH , GDQ , GDCFRC, GDU , GDV , & ! input - CBMFX , GCYT , DELPI , GCHT , GCQT , & ! input - GCLT , GCIT , GCUT , GCVT , GDQ(1,1,iti),& ! input + CBMFX , GCYT , DELPInv , GCHT , GCQT , & ! input + GCLT , GCIT , GCUT , GCVT , GDQ(:,:,iti),& ! input gctrt , & KT , ISTS , IENS, nctp ) ! input endif !for now area fraction of the downdraft is zero, it will be computed -! within cumdwn and applied there -! Get AW downdraft eddy flux and microphysical tendencies out of downdraft code. +! within cumdwn and applied there. So we will get the total sigma now before calling it, +! and apply to the diabatic terms from the updrafts. - do k=1,kmax - do i=ists,iens - sigmad(i,k) = zero - enddo - enddo +! if (do_aw.and.flx_form) then + if (flx_form) then + do k=1,kp1 + do i=ists,iens + lamdamax = maxval(lamdai(i,k,:)) + do while (lamdamax > zero) + loclamdamax = maxloc(lamdai(i,k,:),dim=1) + lamdaprod(i,k) = lamdaprod(i,k) * (one+lamdai(i,k,loclamdamax)) + sigmai(i,k,loclamdamax) = lamdai(i,k,loclamdamax) / lamdaprod(i,k) + sigma(i,k) = max(zero, min(one, sigma(i,k) + sigmai(i,k,loclamdamax))) + vverti(i,k,loclamdamax) = sigmai(i,k,loclamdamax) * wcv(i,k,loclamdamax) + + ! make this lamdai negative so it won't be counted again + lamdai(i,k,loclamdamax) = -lamdai(i,k,loclamdamax) + ! get new lamdamax + lamdamax = maxval(lamdai(i,k,:)) + enddo + ! restore original values of lamdai + lamdai(i,k,:) = abs(lamdai(i,k,:)) +! write(6,'(i2,14f7.4)') k,sigmai(i,k,:) + enddo + enddo + endif + +! the condensation terms - these come from the MSE and condensed water budgets for +! an entraining updraft + if(flx_form) then + DO CTP=1,NCTP ! loop over cloud types + dtcondtem(:,:) = zero + dqcondtem(:,:) = zero + dqprectem(:,:) = zero + dfrzprectem(:,:) = zero + dtfrztem(:,:) = zero + do i=ISTS,IENS + cbmfl = cbmfx(i,ctp) + kk = kt(i,ctp) ! cloud top index + if(cbmfl > zero) then ! this should avoid zero wcv in the denominator + kbi = kb(i) ! cloud base index + do k=kbi,kk ! loop from cloud base to cloud top + km1 = k - 1 + rhs_h = zero + rhs_q = zero + if(k > kbi) then +! tem = cbmfl * (one - sigma(i,k)) + tem = cbmfl * (one - 0.5*(sigma(i,k)+sigma(i,km1))) + tem1 = gcym(i,k,ctp) * (one - sigma(i,k)) + tem2 = gcym(i,km1,ctp) * (one - sigma(i,km1)) + rhs_h = cbmfl * (tem1*gchm(i,k,ctp) - (tem2*gchm(i,km1,ctp) & + + GDH(I,Km1)*(tem1-tem2)) ) + rhs_q = cbmfl * (tem1*(gcwm(i,k,ctp)-gcqm(i,k,ctp)) & + - (tem2*(gcwm(i,km1,ctp)-gcqm(i,km1,ctp)) & + + (GDw(I,Km1)-gdq(i,km1,1))*(tem1-tem2)) ) +! + dqcondtem(i,km1) = -rhs_q ! condensation + dqprectem(i,km1) = tem * (GPRCIZ(i,k,ctp) + GSNWIZ(i,k,ctp)) ! total precip production + dfrzprectem(i,km1) = tem * GSNWIZ(i,k,ctp) ! production of frozen precip + dtfrztem(i,km1) = rhs_h*oneocp ! heating due to freezing +! total temperature tendency due to in cloud microphysics + dtcondtem(i,km1) = - elocp * dqcondtem(i,km1) + dtfrztem(i,km1) + + endif ! if(k > kbi) then + enddo ! end of k=kbi,kk loop + + endif ! end of if(cbmfl > zero) + + +! get tendencies by difference of fluxes, sum over cloud type + + do k = 1,kk +! sum single cloud microphysical tendencies over all cloud types + condtermt(i,k) = condtermt(i,k) + dtcondtem(i,k) * delpinv(i,k) + condtermq(i,k) = condtermq(i,k) + dqcondtem(i,k) * delpinv(i,k) + prectermq(i,k) = prectermq(i,k) + dqprectem(i,k) * delpinv(i,k) + prectermfrz(i,k) = prectermfrz(i,k) + dfrzprectem(i,k) * delpinv(i,k) + frzterm(i,k) = frzterm(i,k) + dtfrztem(i,k) * delpinv(i,k) + +! if (lprnt .and. i == ipr) write(0,*)' k=',k,' trfluxtem=',trfluxtem(k+1,ntr),trfluxtem(k,ntr),& +! ' ctp=',ctp,' trfluxterm=',trfluxterm(i,k,ntr) + enddo + + enddo ! end of i loop + enddo ! end of nctp loop + endif +!downdraft sigma and mass-flux tendency terms are now put into +! the nctp+1 slot of the cloud-type dimensiond variables + + do k=1,kmax + do i=ists,iens + sigmad(i,k) = zero + enddo + enddo !> -# Call cumdwn() to compute cumulus downdraft and assocated melt, freeze !! and evaporation - CALL CUMDWN(IM , IJSDIM, KMAX , NTR , ntrq , & ! DD dimensions + CALL CUMDWN(IM, IJSDIM, KMAX, NTR, ntrq, nctp, & ! DD dimensions GTT , GTQ , GTU , GTV , & ! modified GMFLX , & ! modified updraft+downdraft flux GPRCP , GSNWP , GTEVP , GMDD , & ! output GPRCI , GSNWI , & ! input - GDH , GDW , GDQ , GDQ(1,1,iti) , & ! input + GDH , GDW , GDQ , GDQ(:,:,iti) , & ! input GDQS , GDS , GDHS , GDT , & ! input GDU , GDV , GDZ , & ! input - GDZM , FDQS , DELP , DELPI , & ! input + GDZM , FDQS , DELP , DELPInv , & ! input sigmad, do_aw , do_awdd, flx_form, & ! DDsigma input dtmelt, dtevap, dtsubl, & ! DDsigma input dtdwn , dqvdwn, dqldwn, dqidwn, & ! DDsigma input dtrdwn, & KB , KTMXT , ISTS , IENS ) ! input -! here we substitute the AW tendencies into tendencies to be passed out -! if (do_aw) then -! do k=1,kmax -! do i=ists,iens -! sigma(i,k) = sigma(i,k) + sigmad(i,k) -! enddo -! enddo + +! sigma = sigma + sigmad !> -# Call cumsbw() to compute cloud subsidence heating if (.not. flx_form) then @@ -1381,20 +1451,20 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions CALL CUMSBH(IM , IJSDIM, KMAX , NTR , ntrq , & !DD dimensions GTT , GTQ , & ! modified GTU , GTV , & ! modified - GDH , GDQ , GDQ(1,1,iti) , & ! input + GDH , GDQ , GDQ(:,:,iti) , & ! input GDU , GDV , & ! input - DELPI , GMFLX , GMFX0 , & ! input + DELPINV , GMFLX , GMFX0 , & ! input KTMXT , CPRES , kb, ISTS , IENS ) ! input else CALL CUMSBW(IM , IJSDIM, KMAX , & !DD dimensions GTU , GTV , & ! modified GDU , GDV , & ! input - DELPI , GMFLX , GMFX0 , & ! input + DELPINV , GMFLX , GMFX0 , & ! input KTMXT , CPRES , kb, ISTS , IENS ) ! input endif ! -! for now the following routines appear to be of no consequence to AW - DD +! for now the following routines appear to be of no consequence - DD ! if (.not. flx_form) then ! Tracer Updraft properties @@ -1411,20 +1481,20 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions GCYM , GCYT , GCQT , GCLT , GCIT , & ! input GTPRT , GTEVP , GTPRC0, & ! input KB , KBMX , KT , KTMX , KTMXT , & ! input - DELPI , OTSPT1, ISTS , IENS, & ! input + DELPInv , OTSPT1, ISTS , IENS, & ! input fscav , fswtr, nctp) ! ! Tracer Change due to Downdraft ! --------------- CALL CUMDNR(im ,IJSDIM , KMAX , NTR , & !DD dimensions GTQ , & ! modified - GDQ , GMDD , DELPI , & ! input + GDQ , GMDD , DELPInv , & ! input KTMXT , OTSPT1, ISTS , IENS ) ! input !! !! Tracer change due to Subsidence !! --------------- !! This will be done by cumsbh, now DD 20170907 -! CALL CUMSBR(im , IJSDIM, KMAX , NTR , & !DD dimensions +! CALL CUMSBR(im , IJSDIM, KMAX , NTR ,NCTP, & !DD dimensions ! GTQ , & ! modified ! GDQ , DELPI , & ! input ! GMFLX , KTMXT , OTSPT2, & ! input @@ -1447,6 +1517,60 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions !> -# Compute AW tendencies of T, ql and qi if(flx_form) then ! compute AW tendencies ! AW lump all heating together, compute qv term + +! sigma interpolated to the layer for condensation, etc. terms, precipitation + if(do_aw) then + do k=1,kmax + kp1 = k+1 + do i=1,ijsdim + fsigma(i,k) = one - half*(sigma(i,k)+sigma(i,kp1)) + enddo + enddo + else + do k=1,kmax+1 + do i=1,ijsdim + fsigma(i,k) = one + enddo + enddo + endif + +! AW adjustment of precip fluxes from downdraft model + if(do_aw) then + kp1 = kmax+1 + DO I=ISTS,IENS + GSNWP( I,kp1 ) = zero + GPRCP( I,kp1 ) = zero + ENDDO + tem1 = cpoemelt/grav + tem2 = cpoel/grav + tem3 = cpoesub/grav + DO K=KMAX,1,-1 + kp1 = k+1 + DO I=ISTS,IENS + tem = -dtmelt(i,k) * delp(i,k) * tem1 + teme = -dtevap(i,k) * delp(i,k) * tem2 + tems = -dtsubl(i,k) * delp(i,k) * tem3 + GSNWP(I,k) = GSNWP(I,kp1) + fsigma(i,k) * (GSNWI(i,k) - tem - tems) + GPRCP(I,k) = GPRCP(I,kp1) + fsigma(i,k) * (GPRCI(i,k) + tem - teme) + ENDDO + ENDDO + endif + + +! some of the above routines have set the tendencies and they need to be +! reinitialized, gtt not needed, but gtq needed Anning 5/25/2020 + do n=1,ntr + do k=1,kmax + do i=1,ijsdim + gtq(i,k,n) = zero + enddo + enddo + enddo +! do k=1,kmax +! do i=1,ijsdim +! gtt(i,k) = zero +! enddo +! enddo do k=1,kmax do i=ists,iens dqevap(i,k) = - dtevap(i,k)*cpoel - dtsubl(i,k)*cpoesub @@ -1454,25 +1578,70 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions dtsubl(i,k) = zero enddo enddo - do i=1,ijsdim - moistening_aw(i) = zero - enddo - tem2 = one / delta + + +! diabatic terms from updraft and downdraft models DO K=1,KMAX DO I=ISTS,IENS tem = frzterm(i,k)*cpoEMELT - prectermfrz(i,k) - gtt(i,k) = dtdwn(i,k) + sfluxterm(i,k) + condtermt(i,k) & - + dtmelt(i,k) + dtevap(i,k) - gtq(i,k,1) = dqvdwn(i,k) + qvfluxterm(i,k) + condtermq(i,k) & - + dqevap(i,k) - gtq(i,k,itl) = dqldwn(i,k) + qlfluxterm(i,k) - condtermq(i,k) & +! gtt(i,k) = gtt(i,k) + fsigma(i,k)*(dtmelt(i,k) + dtevap(i,k)) + condtermt(i,k) +! gtq(i,k,1) = gtq(i,k,1) + fsigma(i,k)*dqevap(i,k) + condtermq(i,k) +! gtq(i,k,itl) = gtq(i,k,itl) - (condtermq(i,k) + prectermq(i,k) + tem) +! gtq(i,k,iti) = gtq(i,k,iti) + tem + gtt(i,k) = dtdwn(i,k) + condtermt(i,k) & + + fsigma(i,k)*(dtmelt(i,k) + dtevap(i,k)) + gtq(i,k,1) = dqvdwn(i,k) + condtermq(i,k) & + + fsigma(i,k) * dqevap(i,k) + gtq(i,k,itl) = dqldwn(i,k) - condtermq(i,k) & - prectermq(i,k) - tem - gtq(i,k,iti) = dqidwn(i,k) + qifluxterm(i,k) + tem + gtq(i,k,iti) = dqidwn(i,k) + tem + ! detrainment terms get zeroed ! gtldet(i,k) = zero ! gtidet(i,k) = zero + ENDDO + ENDDO +!! flux tendencies - compute the vertical flux divergence + DO ctp =1,nctp + DO I=ISTS,IENS + cbmfl = cbmfx(i,ctp) + kk = kt(i,ctp) ! cloud top index + if(cbmfl > zero) then ! this should avoid zero wcv in the denominator + DO K=1,kk + kp1 = k+1 + gtt(i,k) = gtt(i,k) - (fsigma(i,kp1)*sfluxtem(i,kp1,ctp) & + - fsigma(i,k)*sfluxtem(i,k,ctp)) * delpinv(i,k) + gtq(i,k,1) = gtq(i,k,1) - (fsigma(i,kp1)*qvfluxtem(i,kp1,ctp) & + - fsigma(i,k)*qvfluxtem(i,k,ctp)) * delpinv(i,k) + gtq(i,k,itl) = gtq(i,k,itl) - (fsigma(i,kp1)*qlfluxtem(i,kp1,ctp) & + - fsigma(i,k)*qlfluxtem(i,k,ctp)) * delpinv(i,k) + gtq(i,k,iti) = gtq(i,k,iti) - (fsigma(i,kp1)*qifluxtem(i,kp1,ctp) & + - fsigma(i,k)*qifluxtem(i,k,ctp)) * delpinv(i,k) + ENDDO +! replace tracer tendency only if to be advected. + DO n = ntrq,NTR + if (OTSPT2(n)) then + DO K=1,kk + kp1 = k+1 + gtq(i,k,n) = - (fsigma(i,kp1)*trfluxtem(i,kp1,n,ctp) & + - fsigma(i,k)*trfluxtem(i,k,n,ctp)) * delpinv(i,k) + ENDDO + endif + ENDDO + end if + ENDDO + ENDDO + +! if(kdt>4) stop 1000 + DO I=ISTS,IENS + moistening_aw(i) = zero + enddo +! adjust tendencies that will lead to negative water mixing ratios + tem2 = one / delta + DO K=1,KMAX + DO I=ISTS,IENS tem1 = - gdq(i,k,itl)*tem2 if (gtq(i,k,itl) < tem1) then tem3 = gtq(i,k,itl) - tem1 @@ -1504,7 +1673,7 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions if (OTSPT2(n)) then DO K=1,KMAX DO I=ISTS,IENS - gtq(i,k,n) = dtrdwn(i,k,n) + trfluxterm(i,k,n) + gtq(i,k,n) = gtq(i,k,n) + dtrdwn(i,k,n) ENDDO ENDDO endif @@ -1597,46 +1766,74 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions !> -# Ensures conservation of water. !In fact, no adjustment of the precip ! is occuring now which is a good sign! DD - if(flx_form .and. adjustp) then + if(flx_form) then DO I = ISTS, IENS if(gprcp(i,1)+gsnwp(i,1) > 1.e-12_kind_phys) then - moistening_aw(i) = - moistening_aw(i) / (gprcp(i,1)+gsnwp(i,1)) - else - moistening_aw(i) = 1.0 + moistening_aw(i) = -moistening_aw(i) / (gprcp(i,1)+gsnwp(i,1)) +! print*,'moistening_aw',moistening_aw(i) + gprcp(i,:) = gprcp(i,:) * moistening_aw(i) + gsnwp(i,:) = gsnwp(i,:) * moistening_aw(i) endif - ENDDO - do k=1,kmax - DO I = ISTS, IENS - gprcp(i,k) = max(0.0, gprcp(i,k) * moistening_aw(i)) - gsnwp(i,k) = max(0.0, gsnwp(i,k) * moistening_aw(i)) - ENDDO - enddo - + END DO endif + +! second method of determining sfc precip only +! if(flx_form) then +! DO I = ISTS, IENS +! pr_tot = zero +! pr_liq = zero +! pr_ice = zero +! do k = 1,kmax +! pr_tot = pr_tot - (gtq(i,k,1)+gtq(i,k,itl)+gtq(i,k,iti)) * delp(i,k) * gravi +! pr_ice = pr_ice + ( cp*gtt(i,k) + el*gtq(i,k,1) - emelt*gtq(i,k,iti) ) & +! * delp(i,k)*gravi +! enddo + !pr_ice = max( min(pr_tot, pr_ice / (emelt)),zero) +! pr_ice = pr_ice / emelt +! pr_liq = pr_tot - pr_ice +! END DO +! print *,'precip1',pr_tot*86400.,gprcp(1,1)*86400.,gsnwp(1,1)*86400. +! print *,'precip2',pr_tot*86400.,pr_liq*86400.,pr_ice*86400. +! endif + + DO K = 1, KMAX + DO I = ISTS, IENS + GPRCPF( I,K ) = 0.5*( GPRCP( I,K )+GPRCP( I,K+1 ) ) + GSNWPF( I,K ) = 0.5*( GSNWP( I,K )+GSNWP( I,K+1 ) ) + END DO + END DO + ! -! do i=ISTS,IENS -! GPRCC(I,1) = GPRCP(I,1) -! GSNWC(I ) = GSNWP(I,1) -! enddo - do k=1,kmax +! do i=ISTS,IENS +! GPRCC( I,1 ) = GPRCP( I,1 ) +! GSNWC( I ) = GSNWP( I,1 ) +! enddo + +! adjust sfc precip consistently with vertically integrated +! temperature and moisture tendencies + + do k=1,kmax+1 do i=ISTS,IENS GTPRP(I,k) = GPRCP(I,k) + GSNWP(I,k) enddo enddo ! !DD provide GFS with a separate downdraft mass flux - DO K=1,KMAX - DO I=ISTS,IENS - GMFX1(I,K) = GMFX0(I,K) - GMFLX(I,K) - ENDDO - ENDDO -! - if (flx_form) then - deallocate(sfluxterm, qvfluxterm, qlfluxterm, qifluxterm,& - condtermt, condtermq, frzterm, prectermq, & - prectermfrz, dtdwn, dqvdwn, dqldwn, & - dqidwn, trfluxterm, dtrdwn) - endif + if(do_aw) then + DO K = 1, KMAX+1 + DO I = ISTS, IENS + fsigma(i,k) = one - sigma(i,k) + GMFX0( I,K ) = GMFX0( I,K ) * fsigma(i,k) + GMFLX( I,K ) = GMFLX( I,K ) * fsigma(i,k) + END DO + END DO + endif + DO K = 1, KMAX+1 + DO I = ISTS, IENS + GMFX1( I,K ) = GMFX0( I,K ) - GMFLX( I,K ) + END DO + END DO + if (allocated(gprcc)) deallocate(gprcc) ! @@ -1748,28 +1945,6 @@ SUBROUTINE CUMBAS & ! cloud base ENDIF ENDDO ENDDO - DO K=KLCLB+1,KBMAX-1 - DO I=ISTS,IENS - spbl(i) = one - gdpm(i,k) * tx1(i) - IF (kb(i) > k .and. spbl(i) > spblmax) THEN - KB(I) = K - ENDIF - ENDDO - ENDDO -! DO K=KBMAX-1,KLCLB+1,-1 -! DO I=ISTS,IENS -! GAMX = FDQS(I,K) / (one+GAM(I,K)) * oneocp -! QSL(i) = GDQS(I,K) + GAMX * (GDH(I,KLCLB)-GDHS(I,K)) -! spbl(i) = one - gdpm(i,k) * tx1(i) -! IF (GDW(I,KLCLB) >= QSL(i) .and. spbl(i) >= spblcrit & -! .and. spbl(i) < spblcrit*6.0) THEN -! .and. spbl(i) < spblcrit*8.0) THEN -! KB(I) = K + KBOFS -! ENDIF -! ENDDO -! if(lprnt) write(0,*)' k=',k,' gdh1=',gdh(ipr,klclb),' gdhs=',gdhs(ipr,k),' kb=',kb(ipr) & -! ,' GDQS=',GDQS(ipr,k),' GDW=',GDW(ipr,KLCLB),' gdpm=',gdpm(ipr,k),' spbl=',spbl(ipr),' qsl=',qsl(ipr) -! ENDDO ENDIF ! do i=ists,iens @@ -1910,8 +2085,8 @@ SUBROUTINE CUMUP & !! in-cloud properties REAL(kind_phys) ACWF (IJSDIM) !< cloud work function REAL(kind_phys) GCLZ (IJSDIM, KMAX) !< cloud liquid water*eta REAL(kind_phys) GCIZ (IJSDIM, KMAX) !< cloud ice*eta - REAL(kind_phys) GPRCIZ(IJSDIM, KMAX) !< rain generation*eta - REAL(kind_phys) GSNWIZ(IJSDIM, KMAX) !< snow generation*eta + REAL(kind_phys) GPRCIZ(IJSDIM, KMAX+1) !< rain generation*eta + REAL(kind_phys) GSNWIZ(IJSDIM, KMAX+1) !< snow generation*eta REAL(kind_phys) GCYT (IJSDIM) !< norm. mass flux @top REAL(kind_phys) GCHT (IJSDIM) !< cloud top MSE*eta REAL(kind_phys) GCQT (IJSDIM) !< cloud top moisture*eta @@ -1924,7 +2099,7 @@ SUBROUTINE CUMUP & !! in-cloud properties REAL(kind_phys) GCwT (IJSDIM) !< cloud top v*eta INTEGER KT (IJSDIM) !< cloud top INTEGER KTMX !< max of cloud top - REAL(kind_phys) WCV (IJSDIM, KMAX) !< updraft velocity (half lev) !DD sigma make output + REAL(kind_phys) WCV (IJSDIM, KMAX+1) !< updraft velocity (half lev) !DD sigma make output ! ! [MODIFIED] REAL(kind_phys) GCYM (IJSDIM, KMAX) !< norm. mass flux @@ -1980,12 +2155,12 @@ SUBROUTINE CUMUP & !! in-cloud properties ! REAL(kind_phys) ELAR (IJSDIM, KMAX) !< entrainment rate REAL(kind_phys) ELAR !< entrainment rate at mid layer ! - REAL(kind_phys) GCHM (IJSDIM, KMAX) !< cloud MSE (half lev) - REAL(kind_phys) GCWM (IJSDIM, KMAX) !< cloud Qt (half lev) !DDsigmadiag - REAL(kind_phys) GCTM (IJSDIM, KMAX) !< cloud T (half lev) !DDsigmadiag make output - REAL(kind_phys) GCQM (IJSDIM, KMAX) !< cloud q (half lev) !DDsigmadiag make output - REAL(kind_phys) GCLM (IJSDIM, KMAX) !< cloud liquid ( half lev) - REAL(kind_phys) GCIM (IJSDIM, KMAX) !< cloud ice (half lev) + REAL(kind_phys) GCHM (IJSDIM, KMAX+1) !< cloud MSE (half lev) + REAL(kind_phys) GCWM (IJSDIM, KMAX+1) !< cloud Qt (half lev) !DDsigmadiag + REAL(kind_phys) GCTM (IJSDIM, KMAX+1) !< cloud T (half lev) !DDsigmadiag make output + REAL(kind_phys) GCQM (IJSDIM, KMAX+1) !< cloud q (half lev) !DDsigmadiag make output + REAL(kind_phys) GCLM (IJSDIM, KMAX+1) !< cloud liquid ( half lev) + REAL(kind_phys) GCIM (IJSDIM, KMAX+1) !< cloud ice (half lev) REAL(kind_phys) GCUM (IJSDIM, KMAX) !< cloud U (half lev) REAL(kind_phys) GCVM (IJSDIM, KMAX) !< cloud V (half lev) REAL(kind_phys) GCtrM (IJSDIM, KMAX,ntrq:ntr) !< cloud tracer (half lev) @@ -2021,8 +2196,9 @@ SUBROUTINE CUMUP & !! in-cloud properties ! REAL(kind_phys) :: WCCRT = zero !m REAL(kind_phys) :: WCCRT = 0.01 REAL(kind_phys) :: WCCRT = 1.0e-6_kind_phys, wvcrt=1.0e-3_kind_phys - REAL(kind_phys) :: TSICE = 268.15_kind_phys ! compatible with macrop_driver - REAL(kind_phys) :: TWICE = 238.15_kind_phys ! compatible with macrop_driver + REAL(kind_phys) :: TSICE = 273.15_kind_phys ! compatible with macrop_driver + REAL(kind_phys) :: TWICE = 233.15_kind_phys ! compatible with macrop_driver + REAL(kind_phys) :: c1t ! REAL(kind_phys) :: wfn_neg = 0.1 REAL(kind_phys) :: wfn_neg = 0.15 @@ -2033,10 +2209,15 @@ SUBROUTINE CUMUP & !! in-cloud properties REAL(kind_phys) :: esat, tem ! REAL(kind_phys) :: esat, tem, rhs_h, rhs_q ! +! [INTERNAL FUNC] + REAL(kind_phys) FPREC ! precipitation ratio in condensate + REAL(kind_phys) FRICE ! ice ratio in cloud water REAL(kind_phys) Z ! altitude REAL(kind_phys) ZH ! scale height REAL(kind_phys) T ! temperature ! + FPREC(Z,ZH) = MIN(MAX(one-EXP(-(Z-PRECZ0)/ZH), zero), one) + FRICE(T) = MIN(MAX((TSICE-T)/(TSICE-TWICE), zero), one) ! ! Note: iteration is not made to diagnose cloud ice for simplicity ! @@ -2052,14 +2233,25 @@ SUBROUTINE CUMUP & !! in-cloud properties GCVT (I) = zero GCwT (I) = zero enddo + do k=1,kmax+1 + do i=ists,iens + GPRCIZ(I,k) = zero + GSNWIZ(I,k) = zero + enddo + enddo + do k=1,kmax + do i=ists,iens + WCV (I,k) = unset_kind_phys + GCLM (I,k) = unset_kind_phys + GCIM (I,k) = unset_kind_phys + enddo + enddo do k=1,kmax do i=ists,iens ACWFK (I,k) = unset_kind_phys ACWFN (I,k) = unset_kind_phys GCLZ (I,k) = zero GCIZ (I,k) = zero - GPRCIZ(I,k) = zero - GSNWIZ(I,k) = zero ! GCHMZ (I,k) = zero GCWMZ (I,k) = zero @@ -2070,15 +2262,12 @@ SUBROUTINE CUMUP & !! in-cloud properties ! BUOY (I,k) = unset_kind_phys BUOYM (I,k) = unset_kind_phys - WCV (I,k) = unset_kind_phys GCY (I,k) = unset_kind_phys ! GCHM (I,k) = unset_kind_phys GCWM (I,k) = unset_kind_phys GCTM (I,k) = unset_kind_phys GCQM (I,k) = unset_kind_phys - GCLM (I,k) = unset_kind_phys - GCIM (I,k) = unset_kind_phys GCUM (I,k) = unset_kind_phys GCVM (I,k) = unset_kind_phys enddo @@ -2199,13 +2388,24 @@ SUBROUTINE CUMUP & !! in-cloud properties FDQSM = GDQSM * tem * (fact1 + fact2*tem) ! calculate d(qs)/dT CPGMI = one / (CP + EL*FDQSM) - PRCZH = PRECZH * MIN(GDZTR(I)*ZTREFI, one) - PRECR = FPREC(GDZM(I,K)-GDZMKB(I), PRCZH ) ! wrk = one / GCYM(I,K) DCTM = (GCHMZ(I,K)*wrk - GDHSM) * CPGMI GCQMZ(i) = min((GDQSM+FDQSM*DCTM)*GCYM(I,K), GCWMZ(I,K)) - GTPRMZ(I,K) = PRECR * (GCWMZ(I,K)-GCQMZ(i)) + if(PRECZH > zero) then + PRCZH = PRECZH * MIN(GDZTR(I)*ZTREFI, one) + PRECR = FPREC(GDZM(I,K)-GDZMKB(I), PRCZH ) + GTPRMZ(I,K) = PRECR * (GCWMZ(I,K)-GCQMZ(i)) + else + DELC=GDZ(I,K)-GDZ(I,KM1) + if(gdtm(i,k)>TSICE) then + c1t=c0t*delc + else + c1t=c0t*exp(d0t*(gdtm(i,k)-TSICE))*delc + end if + c1t=min(c1t, one) + GTPRMZ(I,K) = c1t * (GCWMZ(I,K)-GCQMZ(i)) + end if GTPRMZ(I,K) = MAX(GTPRMZ(I,K), GTPRMZ(I,KM1)) GCCMZ = GCWMZ(I,K) - GCQMZ(i) - GTPRMZ(I,K ) DELC = MIN(GCCMZ, zero) @@ -2274,7 +2474,11 @@ SUBROUTINE CUMUP & !! in-cloud properties wrk = one / GCYM(I,K) DCTM = (GCHMZ(I,K)*wrk - GDHSM) * CPGMI GCQMZ(i) = min((GDQSM+FDQSM*DCTM)*GCYM(I,K), GCWMZ(I,K)) - GTPRMZ(I,K) = PRECR * (GCWMZ(I,K)-GCQMZ(i)) + if(PRECZH > zero) then + GTPRMZ(I,K) = PRECR * (GCWMZ(I,K)-GCQMZ(i)) + else + GTPRMZ(I,K) = c1t * (GCWMZ(I,K)-GCQMZ(i)) + end if GTPRMZ(I,K) = MAX(GTPRMZ(I,K), GTPRMZ(I,KM1)) GCCMZ = GCWMZ(I,K) - GCQMZ(i) - GTPRMZ(I,K) DELC = MIN(GCCMZ, zero) @@ -2399,8 +2603,19 @@ SUBROUTINE CUMUP & !! in-cloud properties wrk = one / gcyt(i) DCT = (GCHT(I)*wrk - GDHS(I,K)) / (CP*(one + GAM(I,K))) GCQT(I) = min((GDQS(I,K) + FDQS(I,K)*DCT) * GCYT(I), GCWT(i)) - PRCZH = PRECZH * MIN(GDZTR(I)*ZTREFI, one) - GTPRT(I) = FPREC(GDZ(I,K)-GDZMKB(I), PRCZH) * (GCWT(i)-GCQT(I)) + if(PRECZH > zero) then + PRCZH = PRECZH * MIN(GDZTR(I)*ZTREFI, one) + GTPRT(I) = FPREC(GDZ(I,K)-GDZMKB(I), PRCZH) * (GCWT(i)-GCQT(I)) + else + DELC=GDZ(I,K)-GDZ(I,K-1) + if(gdtm(i,k)>TSICE) then + c1t=c0t*delc + else + c1t=c0t*exp(d0t*(gdtm(i,k)-TSICE))*delc + end if + c1t=min(c1t, one) + GTPRT(I) = c1t * (GCWT(i)-GCQT(I)) + end if GTPRT(I) = MAX(GTPRT(I), GTPRMZ(I,K)) GCCT = GCWT(i) - GCQT(I) - GTPRT(I) DELC = MIN(GCCT, zero) @@ -2503,24 +2718,6 @@ SUBROUTINE CUMUP & !! in-cloud properties ! ! WRITE( CTNUM, '(I2.2)' ) CTP ! - -contains - - pure function FPREC(Z,ZH) - implicit none - real(kind_phys), intent(in) :: Z - real(kind_phys), intent(in) :: ZH - real(kind_phys) :: FPREC - FPREC = MIN(MAX(one-EXP(-(Z-PRECZ0)/ZH), zero), one) - end function FPREC - - pure function FRICE(T) - implicit none - real(kind_phys), intent(in) :: T - real(kind_phys) :: FRICE - FRICE = MIN(MAX((TSICE-T)/(TSICE-TWICE), zero), one) - end function FRICE - END SUBROUTINE CUMUP !*********************************************************************** !>\ingroup cs_scheme @@ -2562,8 +2759,8 @@ SUBROUTINE CUMBMX & !! cloud base mass flux ! [INTERNAL PARAM] REAL(kind_phys) :: FMAX = 1.5e-2_kind_phys ! maximum flux ! REAL(kind_phys) :: RHMCRT = zero ! critical val. of cloud mean RH -! REAL(kind_phys) :: RHMCRT = 0.25_kind_phys ! critical val. of cloud mean RH - REAL(kind_phys) :: RHMCRT = 0.50_kind_phys ! critical val. of cloud mean RH + REAL(kind_phys) :: RHMCRT = 0.25_kind_phys ! critical val. of cloud mean RH +! REAL(kind_phys) :: RHMCRT = 0.50_kind_phys ! critical val. of cloud mean RH REAL(kind_phys) :: ALP1 = zero REAL(kind_phys) :: TAUD = 1.e3_kind_phys ! REAL(kind_phys) :: TAUD = 6.e2_kind_phys @@ -2624,7 +2821,7 @@ SUBROUTINE CUMFLX & !! cloud mass flux INTEGER, INTENT(IN) :: IJSDIM, KMAX, IM !! DD, for GFS, pass in ! ! [OUTPUT] - REAL(kind_phys) GMFLX (IJSDIM, KMAX) !< mass flux + REAL(kind_phys) GMFLX (IJSDIM, KMAX+1) !< mass flux REAL(kind_phys) CMDET (IJSDIM, KMAX) !< detrainment mass flux REAL(kind_phys) GPRCI (IJSDIM, KMAX) !< rainfall generation REAL(kind_phys) GSNWI (IJSDIM, KMAX) !< snowfall generation @@ -2636,8 +2833,8 @@ SUBROUTINE CUMFLX & !! cloud mass flux REAL(kind_phys) CBMFX (IJSDIM) !< cloud base mass flux REAL(kind_phys) GCYM (IJSDIM, KMAX) !< normalized mass flux REAL(kind_phys) GCYT (IJSDIM) !< detraining mass flux - REAL(kind_phys) GPRCIZ(IJSDIM, KMAX) !< precipitation/M - REAL(kind_phys) GSNWIZ(IJSDIM, KMAX) !< snowfall/M + REAL(kind_phys) GPRCIZ(IJSDIM, KMAX+1) !< precipitation/M + REAL(kind_phys) GSNWIZ(IJSDIM, KMAX+1) !< snowfall/M REAL(kind_phys) GTPRT (IJSDIM) !< rain+snow @top REAL(kind_phys) GCLZ (IJSDIM, KMAX) !< cloud liquid/M REAL(kind_phys) GCIZ (IJSDIM, KMAX) !< cloud ice/M @@ -2773,8 +2970,8 @@ SUBROUTINE CUMSBH & !! adiabat. descent REAL(kind_phys) GDU (IJSDIM, KMAX) REAL(kind_phys) GDV (IJSDIM, KMAX) REAL(kind_phys) DELPI (IJSDIM, KMAX) - REAL(kind_phys) GMFLX (IJSDIM, KMAX) !< mass flux (updraft+downdraft) - REAL(kind_phys) GMFX0 (IJSDIM, KMAX) !< mass flux (updraft only) + REAL(kind_phys) GMFLX (IJSDIM, KMAX+1) !< mass flux (updraft+downdraft) + REAL(kind_phys) GMFX0 (IJSDIM, KMAX+1) !< mass flux (updraft only) INTEGER KB(IJSDIM) !< cloud base index - negative means no convection INTEGER KTMX REAL(kind_phys) CPRES !< pressure factor for cumulus friction @@ -2890,8 +3087,8 @@ SUBROUTINE CUMSBW & !! adiabat. descent REAL(kind_phys) GDU (IJSDIM, KMAX) REAL(kind_phys) GDV (IJSDIM, KMAX) REAL(kind_phys) DELPI (IJSDIM, KMAX) - REAL(kind_phys) GMFLX (IJSDIM, KMAX) !< mass flux (updraft+downdraft) - REAL(kind_phys) GMFX0 (IJSDIM, KMAX) !< mass flux (updraft only) + REAL(kind_phys) GMFLX (IJSDIM, KMAX+1) !< mass flux (updraft+downdraft) + REAL(kind_phys) GMFX0 (IJSDIM, KMAX+1) !< mass flux (updraft only) INTEGER KB(IJSDIM) !< cloud base index - negative means no convection INTEGER KTMX, ISTS, IENS REAL(kind_phys) CPRES !< pressure factor for cumulus friction @@ -2962,7 +3159,7 @@ SUBROUTINE CUMDWN & ! Freeze & Melt & Evaporation ! IMPLICIT NONE - INTEGER, INTENT(IN) :: IM, IJSDIM, KMAX, NTR, ntrq ! DD, for GFS, pass in + INTEGER, INTENT(IN) :: IM, IJSDIM, KMAX, NTR , ntrq, nctp !! DD, for GFS, pass in logical, intent(in) :: do_aw, do_awdd, flx_form ! ! [MODIFY] @@ -2970,13 +3167,13 @@ SUBROUTINE CUMDWN & ! Freeze & Melt & Evaporation REAL(kind_phys) GTQ (IJSDIM, KMAX, NTR) !< Moisture etc tendency REAL(kind_phys) GTU (IJSDIM, KMAX) !< u tendency REAL(kind_phys) GTV (IJSDIM, KMAX) !< v tendency - REAL(kind_phys) GMFLX (IJSDIM, KMAX) !< mass flux + REAL(kind_phys) GMFLX (IJSDIM, KMAX+1) !< mass flux ! ! [OUTPUT] - REAL(kind_phys) GPRCP (IJSDIM, KMAX) !< rainfall flux - REAL(kind_phys) GSNWP (IJSDIM, KMAX) !< snowfall flux + REAL(kind_phys) GPRCP (IJSDIM, KMAX+1) !< rainfall flux + REAL(kind_phys) GSNWP (IJSDIM, KMAX+1) !< snowfall flux REAL(kind_phys) GTEVP (IJSDIM, KMAX) !< evaporation+sublimation - REAL(kind_phys) GMDD (IJSDIM, KMAX) !< downdraft mass flux + REAL(kind_phys) GMDD (IJSDIM, KMAX+1) !< downdraft mass flux !AW microphysical tendencies REAL(kind_phys) gtmelt (IJSDIM, KMAX) !< t tendency ice-liq @@ -2988,8 +3185,6 @@ SUBROUTINE CUMDWN & ! Freeze & Melt & Evaporation REAL(kind_phys) dqldwn (IJSDIM, KMAX) !< ql tendency downdraft detrainment REAL(kind_phys) dqidwn (IJSDIM, KMAX) !< qi tendency downdraft detrainment REAL(kind_phys) dtrdwn (IJSDIM, KMAX, ntrq:ntr) !< tracer tendency downdraft detrainment -! AW downdraft area fraction (assumed zero for now) - REAL(kind_phys) sigmad (IJSDIM,KMAX) !< DDsigma cloud downdraft area fraction ! [INPUT] REAL(kind_phys) GPRCI (IJSDIM, KMAX) !< rainfall generation @@ -3011,6 +3206,8 @@ SUBROUTINE CUMDWN & ! Freeze & Melt & Evaporation REAL(kind_phys) DELPI (IJSDIM, KMAX) INTEGER KB (IJSDIM) INTEGER KTMX, ISTS, IENS + REAL(kind_phys) sigmad (IM,KMAX+1) !< DDsigma cloud downdraft area fraction + ! ! [INTERNAL WORK] ! Note: Some variables have 3-dimensions for the purpose of budget check. @@ -3031,27 +3228,33 @@ SUBROUTINE CUMDWN & ! Freeze & Melt & Evaporation ! profiles of downdraft variables for AW flux tendencies REAL(kind_phys) GCdseD(ISTS:IENS, KMAX) !< downdraft dse REAL(kind_phys) GCqvD (ISTS:IENS, KMAX) !< downdraft qv -! REAL(kind_phys) GCqlD (ISTS:IENS, KMAX) !< downdraft ql -! REAL(kind_phys) GCqiD (ISTS:IENS, KMAX) !< downdraft qi + REAL(kind_phys) GCqlD (ISTS:IENS, KMAX) !< downdraft ql + REAL(kind_phys) GCqiD (ISTS:IENS, KMAX) !< downdraft qi REAL(kind_phys) GCtrD (ISTS:IENS, ntrq:ntr) !< downdraft tracer REAL(kind_phys) GCUD (ISTS:IENS) !< downdraft u REAL(kind_phys) GCVD (ISTS:IENS) !< downdraft v REAL(kind_phys) FSNOW (ISTS:IENS) REAL(kind_phys) GMDDD (ISTS:IENS) - - REAL(kind_phys) GDTW, GCHX, GCTX, GCQSX, GTPRP, EVSU, GTEVE, LVIC, & - DQW, DTW, GDQW, DZ, GCSD, FDET, GDHI, GMDDX, & - GMDDMX, GCHDX, GCWDX, GCUDD, GCVDD, GTHCI, GTQVCI, & - wrk, wrk1, wrk2, wrk3, wrk4, tx1, & - WMX, HMX, DDWMX, DDHMX, dp_above, dp_below, fsigma, & - fmelt, fevp, wrkn, gctrdd(ntrq:ntr) - + INTEGER I, K + REAL(kind_phys) GDTW + REAL(kind_phys) GCHX, GCTX, GCQSX, GTPRP, EVSU, GTEVE, LVIC + REAL(kind_phys) DQW, DTW, GDQW, DZ, GCSD, FDET, GDHI + REAL(kind_phys) GMDDX, GMDDMX + REAL(kind_phys) GCHDX, GCWDX + REAL(kind_phys) GCUDD, GCVDD + REAL(kind_phys) GTHCI, GTQVCI, GTQLCI, GTQICI !M REAL(kind_phys) GTHCI, GTQVCI, GTQLCI, GTQICI, GTUCI, GTVCI + real(kind_phys) wrk, fmelt, fevp, gctrdd(ntrq:ntr) !DD#ifdef OPT_CUMBGT -! Water, energy, downdraft water and downdraft energy budgets -! REAL(kind_phys), dimension(ISTS:IENS) :: WBGT, HBGT, DDWBGT, DDHBGT - integer ij, i, k, kp1, n + REAL(kind_phys) WBGT ( ISTS:IENS ) !! water budget + REAL(kind_phys) HBGT ( ISTS:IENS ) !! energy budget + REAL(kind_phys) DDWBGT( ISTS:IENS ) !! downdraft water budget + REAL(kind_phys) DDHBGT( ISTS:IENS ) !! downdraft energy budget + REAL(kind_phys) WMX, HMX, DDWMX, DDHMX, tx1, wrk1, wrk2, wrk3, wrk4, wrkn + REAL(kind_phys) dp_above, dp_below + real(kind_phys) fsigma + integer ij, n, kp1 !DD#endif ! ! [INTERNAL PARM] @@ -3109,46 +3312,23 @@ SUBROUTINE CUMDWN & ! Freeze & Melt & Evaporation gtsubl(I,k) = zero enddo enddo -! testing on oct 17 2016 - if (flx_form) then - if (.not. do_awdd) then - do k=1,kmax - do i=ists,iens - if (kb(i) > 0) then - dtdwn (i,k) = gtt(i,k) - dqvdwn(i,k) = gtq(i,k,1) - dqldwn(i,k) = gtq(i,k,itl) - dqidwn(i,k) = gtq(i,k,iti) - endif - enddo - enddo - do n=ntrq,ntr - do k=1,kmax - do i=ists,iens - if (kb(i) > 0) then - dtrdwn(i,k,n) = gtq(i,k,n) - endif - enddo - enddo - enddo - else - do k=1,kmax - do i=ists,iens - dtdwn (I,k) = zero - dqvdwn(I,k) = zero - dqldwn(I,k) = zero - dqidwn(I,k) = zero - enddo - enddo - do n=ntrq,ntr - do k=1,kmax - do i=ists,iens - dtrdwn(i,k,n) = zero - enddo - enddo - enddo - endif - endif + +! These are zeroed by the calling routine, cs_cumlus +! do k=1,kmax +! do i=ists,iens +! dtdwn (I,k) = zero +! dqvdwn(I,k) = zero +! dqldwn(I,k) = zero +! dqidwn(I,k) = zero +! enddo +! enddo +! do n=ntrq,ntr +! do k=1,kmax +! do i=ists,iens +! dtrdwn(i,k,n) = zero +! enddo +! enddo +! enddo ! do i=ists,iens GCHD(I) = zero @@ -3178,20 +3358,19 @@ SUBROUTINE CUMDWN & ! Freeze & Melt & Evaporation LVIC = ELocp + EMELTocp*FSNOW(I) GDTW = GDT(I,K) - LVIC*(GDQS(I,K) - GDQ(I,K,1)) & / (one + LVIC*FDQS(I,K)) + + DZ = GDZM(I,KP1) - GDZM(I,K) + FMELT = (one + FTMLT*(GDTW - TWSNOW)) & + * (one - TANH(GMFLX(I,KP1)/GMFLXC)) & + * (one - TANH(VTERMS*MELTAU/DZ)) IF (GDTW < TWSNOW) THEN - GSNWP(I,K) = GSNWP(I,KP1) + GPRCI(I,K) + GSNWI(I,K) - GTTEV(I,K) = EMELToCP * GPRCI(I,K) * DELPI(I,K) - SNMLT(I,K) = -GPRCI(I,K) + SNMLT(I,K) = GPRCP(I,KP1)*min(max(FMELT, one), zero) ELSE - DZ = GDZM(I,KP1) - GDZM(I,K) - FMELT = (one + FTMLT*(GDTW - TWSNOW)) & - * (one - TANH(GMFLX(I,KP1)/GMFLXC)) & - * (one - TANH(VTERMS*MELTAU/DZ)) SNMLT(I,K) = GSNWP(I,KP1)*max(min(FMELT, one), zero) - GSNWP(I,K) = GSNWP(I,KP1)+GSNWI(I,K) - SNMLT(I,K) - GPRCP(I,K) = GPRCP(I,KP1)+GPRCI(I,K) + SNMLT(I,K) - GTTEV(I,K) = -EMELToCP * SNMLT(I,K) * DELPI(I,K) ENDIF + GSNWP(I,K) = GSNWP(I,KP1)+GSNWI(I,K) - SNMLT(I,K) + GPRCP(I,K) = GPRCP(I,KP1)+GPRCI(I,K) + SNMLT(I,K) + GTTEV(I,K) = -EMELToCP * SNMLT(I,K) * DELPI(I,K) !DD heating rate due to precip melting for AW gtmelt(i,k) = gtmelt(i,k) + GTTEV(I,K) endif @@ -3350,8 +3529,15 @@ SUBROUTINE CUMDWN & ! Freeze & Melt & Evaporation GTQ(I,K,1) = GTQ(I,K,1) + GTQEV(I,K) ! GMFLX(I,K) = GMFLX(I,K) - GMDD(I,K) + endif + ENDDO ! end of i loop + ENDDO ! end of k loop ! AW tendencies due to vertical divergence of eddy fluxes + DO K=2,KTMX + kp1 = min(k+1,kmax) + DO I=ISTS,IENS + if (kb(i) > 0) then if (k > 1 .and. flx_form) then fsigma = one - sigmad(i,kp1) dp_below = wrk * (one - sigmad(i,k)) @@ -3381,28 +3567,6 @@ SUBROUTINE CUMDWN & ! Freeze & Melt & Evaporation endif ENDDO ! end of i loop ENDDO ! end of k loop -! - if (.not. do_awdd .and. flx_form) then - do k=1,kmax - do i=ists,iens - if (kb(i) > 0) then - dtdwn(i,k) = gtt(i,k) - dtdwn(i,k) - dqvdwn(i,k) = gtq(i,k,1) - dqvdwn(i,k) - dqldwn(i,k) = gtq(i,k,itl) - dqldwn(i,k) - dqidwn(i,k) = gtq(i,k,iti) - dqidwn(i,k) - endif - enddo - enddo - do n=ntrq,ntr - do k=1,kmax - do i=ists,iens - if (kb(i) > 0) then - dtrdwn(i,k,n) = gtq(i,k,n) - dtrdwn(i,k,n) - endif - enddo - enddo - enddo - endif ! END SUBROUTINE CUMDWN !*********************************************************************** @@ -3428,22 +3592,28 @@ SUBROUTINE CUMCLD & !! cloudiness REAL(kind_phys) FLIQC (IJSDIM, KMAX) !< liquid ratio in cumulus ! ! [INPUT] - REAL(kind_phys) GMFLX (IJSDIM, KMAX) ! cumulus mass flux + REAL(kind_phys) GMFLX (IJSDIM, KMAX+1) ! cumulus mass flux INTEGER KTMX INTEGER ISTS, IENS ! ! [WORK] INTEGER I, K REAL(kind_phys) CUMF, QC, wrk + LOGICAL, SAVE :: OFIRST = .TRUE. ! ! [INTERNAL PARAM] - REAL(kind_phys), parameter :: CMFMIN = 2.e-3_kind_phys, &!< Mc->cloudiness - CMFMAX = 3.e-1_kind_phys, &!< Mc->cloudiness - CLMIN = 1.e-3_kind_phys, &!< cloudiness Min. - CLMAX = 0.1_kind_phys, &!< cloudiness Max. - FACLW = 0.1_kind_phys, &!< Mc->CLW - FACLF = (CLMAX-CLMIN)/LOG(CMFMAX/CMFMIN) -! + REAL(kind_phys) :: FACLW = 0.1_kind_phys !> Mc->CLW + REAL(kind_phys) :: CMFMIN = 2.e-3_kind_phys !> Mc->cloudiness + REAL(kind_phys) :: CMFMAX = 3.e-1_kind_phys !> Mc->cloudiness + REAL(kind_phys) :: CLMIN = 1.e-3_kind_phys !> cloudiness Min. + REAL(kind_phys) :: CLMAX = 0.1_kind_phys !> cloudiness Max. + REAL(kind_phys), SAVE :: FACLF +! + IF ( OFIRST ) THEN + FACLF = (CLMAX-CLMIN)/LOG(CMFMAX/CMFMIN) + OFIRST = .FALSE. + END IF + CUMFRC(ISTS:IENS) = zero DO K=1,KTMX DO I=ISTS,IENS @@ -3668,26 +3838,28 @@ END SUBROUTINE CUMDNR !*********************************************************************** !>\ingroup cs_scheme SUBROUTINE CUMSBR & !! Tracer Subsidence - ( IM , IJSDIM, KMAX , NTR , & !DD dimensions + ( IM , IJSDIM, KMAX, NTR, NCTP, & !DD dimensions GTR , & ! modified - GDR , DELPI , & ! input + GDR , DELP , & ! input GMFLX , KTMX , OTSPT , & ! input + sigmai , sigma , & !DDsigma input ISTS, IENS ) ! input ! IMPLICIT NONE - INTEGER, INTENT(IN) :: IM, IJSDIM, KMAX, NTR !! DD, for GFS, pass in + INTEGER, INTENT(IN) :: IM, IJSDIM, KMAX, NTR, nctp !! DD, for GFS, pass in ! ! [MODIFY] REAL(kind_phys) GTR (IJSDIM, KMAX, NTR) !! tracer tendency ! ! [INPUT] REAL(kind_phys) GDR (IJSDIM, KMAX, NTR) !! tracer - REAL(kind_phys) DELPI (IJSDIM, KMAX) - REAL(kind_phys) GMFLX (IJSDIM, KMAX) !! mass flux + REAL(kind_phys) DELP (IJSDIM, KMAX) + REAL(kind_phys) GMFLX (IJSDIM, KMAX+1) !! mass flux INTEGER KTMX LOGICAL OTSPT (NTR) !! tracer transport on/off INTEGER ISTS, IENS + REAL(kind_phys) sigmai (IM,KMAX+1,NCTP), sigma(IM,KMAX+1) !!DDsigma cloud updraft fraction ! ! [INTERNAL WORK] INTEGER I, K, KM, KP, LT @@ -3703,14 +3875,14 @@ SUBROUTINE CUMSBR & !! Tracer Subsidence KM = MAX(K-1, 1) KP = MIN(K+1, KMAX) DO I=ISTS,IENS - SBR0 = GMFLX(I,KP) * (GDR(I,KP,LT) - GDR(I,K,LT)) - SBR1 = GMFLX(I,K) * (GDR(I,K,LT) - GDR(I,KM,LT)) - IF (GMFLX(I,K) > GMFLX(I,KP)) THEN + SBR0 = GMFLX(I,K+1) * (GDR(I,KP,LT) - GDR(I,K,LT)) + SBR1 = GMFLX(I,K) * (GDR(I,K,LT) - GDR(I,KM,LT)) + IF (GMFLX(I,K) > GMFLX(I,K+1)) THEN FX1 = half ELSE FX1 = zero END IF - GTR(I,K,LT) = GTR(I,K,LT) + DELPI(I,K) & + GTR(I,K,LT) = GTR(I,K,LT) + GRAV/DELP(I,K) & * ((one-FX(I))*SBR0 + FX1*SBR1) FX(I) = FX1 ENDDO @@ -3733,7 +3905,7 @@ SUBROUTINE CUMFXR & ! Tracer mass fixe INTEGER, INTENT(IN) :: IM, IJSDIM, KMAX, NTR !! DD, for GFS, pass in ! ! [MODIFY] - REAL(kind_phys) GTR (IJSDIM, KMAX, NTR) ! tracer tendency + REAL(kind_phys) GTR (IJSDIM, KMAX) ! tracer tendency ! ! [INPUT] REAL(kind_phys) GDR (IJSDIM, KMAX, NTR) ! tracer @@ -3815,14 +3987,14 @@ END SUBROUTINE CUMFXR !********************************************************************* !>\ingroup cs_scheme SUBROUTINE CUMFXR1 & ! Tracer mass fixer - ( IM , IJSDIM, KMAX , & !DD dimensions + ( IM , IJSDIM, KMAX ,nctp, & !DD dimensions GTR , & ! modified GDR , DELP , DELTA , KTMX , IMFXR , & ! input ISTS , IENS ) ! input ! IMPLICIT NONE - INTEGER, INTENT(IN) :: IM, IJSDIM, KMAX ! DD, for GFS, pass in + INTEGER, INTENT(IN) :: IM, IJSDIM, KMAX, nctp !! DD, for GFS, pass in ! ! [MODIFY] REAL(kind_phys) GTR (IJSDIM, KMAX) ! tracer tendency From 633f8d1d5bc8f58b1e17a7196f9d0e698133aba4 Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Thu, 7 Mar 2024 13:57:28 -0500 Subject: [PATCH 13/21] fix array dimension issue --- physics/CONV/Chikira_Sugiyama/cs_conv.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/physics/CONV/Chikira_Sugiyama/cs_conv.F90 b/physics/CONV/Chikira_Sugiyama/cs_conv.F90 index 94dadba87..4ca545bc4 100644 --- a/physics/CONV/Chikira_Sugiyama/cs_conv.F90 +++ b/physics/CONV/Chikira_Sugiyama/cs_conv.F90 @@ -3139,7 +3139,7 @@ END SUBROUTINE CUMSBW !>\ingroup cs_scheme !! This subroution calculates freeze, melt and evaporation in cumulus downdraft. SUBROUTINE CUMDWN & ! Freeze & Melt & Evaporation - ( IM , IJSDIM, KMAX , NTR , ntrq, & !DD dimensions + ( IM , IJSDIM, KMAX , NTR,ntrq,nctp, & !DD dimensions GTT , GTQ , GTU , GTV , & ! modified GMFLX , & ! modified GPRCP , GSNWP , GTEVP , GMDD , & ! output @@ -3905,7 +3905,7 @@ SUBROUTINE CUMFXR & ! Tracer mass fixe INTEGER, INTENT(IN) :: IM, IJSDIM, KMAX, NTR !! DD, for GFS, pass in ! ! [MODIFY] - REAL(kind_phys) GTR (IJSDIM, KMAX) ! tracer tendency + REAL(kind_phys) GTR (IJSDIM, KMAX, NTR) ! tracer tendency ! ! [INPUT] REAL(kind_phys) GDR (IJSDIM, KMAX, NTR) ! tracer From 12b5882e675c70c10468ce6b7ad935bdd51747b5 Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Thu, 7 Mar 2024 14:32:50 -0500 Subject: [PATCH 14/21] remove some unnecessary comment edits and whitespace changes --- physics/CONV/Chikira_Sugiyama/cs_conv.F90 | 193 +++++++++++----------- 1 file changed, 96 insertions(+), 97 deletions(-) diff --git a/physics/CONV/Chikira_Sugiyama/cs_conv.F90 b/physics/CONV/Chikira_Sugiyama/cs_conv.F90 index 4ca545bc4..ae59db7ba 100644 --- a/physics/CONV/Chikira_Sugiyama/cs_conv.F90 +++ b/physics/CONV/Chikira_Sugiyama/cs_conv.F90 @@ -663,96 +663,96 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions ! ! [INTERNAL WORK] REAL(kind_phys), allocatable :: GPRCC (:, :) ! rainfall - REAL(kind_phys) GSNWC ( IJSDIM ) !! snowfall - REAL(kind_phys) CUMCLW( IJSDIM, KMAX ) !! cloud water in cumulus - REAL(kind_phys) CUMFRC( IJSDIM ) !! cumulus cloud fraction + REAL(kind_phys) GSNWC ( IJSDIM ) ! snowfall + REAL(kind_phys) CUMCLW( IJSDIM, KMAX ) ! cloud water in cumulus + REAL(kind_phys) CUMFRC( IJSDIM ) ! cumulus cloud fraction !COSP - REAL(kind_phys) QLIQC ( IJSDIM, KMAX ) !! cumulus cloud liquid water [kg/kg] - REAL(kind_phys) QICEC ( IJSDIM, KMAX ) !! cumulus cloud ice [kg/kg] - REAL(kind_phys) GPRCPF( IJSDIM, KMAX ) !! rainfall flux at full level - REAL(kind_phys) GSNWPF( IJSDIM, KMAX ) !! snowfall flux at full level + REAL(kind_phys) QLIQC ( IJSDIM, KMAX ) ! cumulus cloud liquid water [kg/kg] + REAL(kind_phys) QICEC ( IJSDIM, KMAX ) ! cumulus cloud ice [kg/kg] + REAL(kind_phys) GPRCPF( IJSDIM, KMAX ) ! rainfall flux at full level + REAL(kind_phys) GSNWPF( IJSDIM, KMAX ) ! snowfall flux at full level ! - REAL(kind_phys) GTCFRC( IJSDIM, KMAX ) !! change in cloud fraction - REAL(kind_phys) FLIQC ( IJSDIM, KMAX ) !! liquid ratio in cumulus + REAL(kind_phys) GTCFRC( IJSDIM, KMAX ) ! change in cloud fraction + REAL(kind_phys) FLIQC ( IJSDIM, KMAX ) ! liquid ratio in cumulus ! !#ifdef OPT_CHASER -! REAL(kind_phys) RFXC ( IJSDIM, KMAX+1 ) !! precipi. flx [kg/m2/s] -! REAL(kind_phys) SFXC ( IJSDIM, KMAX+1 ) !! ice/snow flx [kg/m2/s] -! INTEGER LEVCUM( IJSDIM, KMAX ) !! flag for cum. cloud top -! REAL(kind_phys) LNFRC ( IJSDIM, KMAX ) !! areal rates of clouds -! REAL(kind_phys) REVC ( IJSDIM, KMAX ) !! evaporation rates +! REAL(kind_phys) RFXC ( IJSDIM, KMAX+1 ) ! precipi. flx [kg/m2/s] +! REAL(kind_phys) SFXC ( IJSDIM, KMAX+1 ) ! ice/snow flx [kg/m2/s] +! INTEGER LEVCUM( IJSDIM, KMAX ) ! flag for cum. cloud top +! REAL(kind_phys) LNFRC ( IJSDIM, KMAX ) ! areal rates of clouds +! REAL(kind_phys) REVC ( IJSDIM, KMAX ) ! evaporation rates !#endif ! - REAL(kind_phys) GDCFRC( IJSDIM, KMAX ) !! cloud fraction + REAL(kind_phys) GDCFRC( IJSDIM, KMAX ) ! cloud fraction ! -! REAL(kind_phys) GTQL ( IJSDIM, KMAX ) !! tendency of cloud liquid +! REAL(kind_phys) GTQL ( IJSDIM, KMAX ) ! tendency of cloud liquid ! - REAL(kind_phys) GDW ( IJSDIM, KMAX ) !! total water - REAL(kind_phys) GDQS ( IJSDIM, KMAX ) !! saturate moisture + REAL(kind_phys) GDW ( IJSDIM, KMAX ) ! total water + REAL(kind_phys) GDQS ( IJSDIM, KMAX ) ! saturate moisture REAL(kind_phys) FDQS ( IJSDIM, KMAX ) REAL(kind_phys) GAM ( IJSDIM, KMAX ) - REAL(kind_phys) GDS ( IJSDIM, KMAX ) !! dry static energy - REAL(kind_phys) GDH ( IJSDIM, KMAX ) !! moist static energy - REAL(kind_phys) GDHS ( IJSDIM, KMAX ) !! saturate MSE -! - REAL(kind_phys) GCYM ( IJSDIM, KMAX, NCTP ) !! norm. mass flux (half lev) - REAL(kind_phys) GCHB ( IJSDIM ) !! cloud base MSE-Li*Qi - REAL(kind_phys) GCWB ( IJSDIM ) !! cloud base total water - REAL(kind_phys) GCtrB ( IJSDIM, ntrq:ntr ) !! cloud base water vapor tracer - REAL(kind_phys) GCUB ( IJSDIM ) !! cloud base U - REAL(kind_phys) GCVB ( IJSDIM ) !! cloud base V - REAL(kind_phys) GCIB ( IJSDIM ) !! cloud base ice - REAL(kind_phys) ELAM ( IJSDIM, KMAX, NCTP ) !! entrainment (rate*massflux) - REAL(kind_phys) GCYT ( IJSDIM, NCTP ) !! norm. mass flux @top - REAL(kind_phys) GCHT ( IJSDIM, NCTP ) !! cloud top MSE - REAL(kind_phys) GCQT ( IJSDIM, NCTP ) !! cloud top q - REAL(kind_phys) GCwT ( IJSDIM ) !! cloud top total water - REAL(kind_phys) GCUT ( IJSDIM, NCTP ) !! cloud top U - REAL(kind_phys) GCVT ( IJSDIM, NCTP ) !! cloud top V - REAL(kind_phys) GCLT ( IJSDIM, NCTP ) !! cloud top cloud water - REAL(kind_phys) GCIT ( IJSDIM, NCTP ) !! cloud top cloud ice + REAL(kind_phys) GDS ( IJSDIM, KMAX ) ! dry static energy + REAL(kind_phys) GDH ( IJSDIM, KMAX ) ! moist static energy + REAL(kind_phys) GDHS ( IJSDIM, KMAX ) ! saturate MSE +! + REAL(kind_phys) GCYM ( IJSDIM, KMAX, NCTP ) ! norm. mass flux (half lev) + REAL(kind_phys) GCHB ( IJSDIM ) ! cloud base MSE-Li*Qi + REAL(kind_phys) GCWB ( IJSDIM ) ! cloud base total water + REAL(kind_phys) GCtrB ( IJSDIM, ntrq:ntr ) ! cloud base water vapor tracer + REAL(kind_phys) GCUB ( IJSDIM ) ! cloud base U + REAL(kind_phys) GCVB ( IJSDIM ) ! cloud base V + REAL(kind_phys) GCIB ( IJSDIM ) ! cloud base ice + REAL(kind_phys) ELAM ( IJSDIM, KMAX, NCTP ) ! entrainment (rate*massflux) + REAL(kind_phys) GCYT ( IJSDIM, NCTP ) ! norm. mass flux @top + REAL(kind_phys) GCHT ( IJSDIM, NCTP ) ! cloud top MSE + REAL(kind_phys) GCQT ( IJSDIM, NCTP ) ! cloud top q + REAL(kind_phys) GCwT ( IJSDIM ) ! cloud top total water + REAL(kind_phys) GCUT ( IJSDIM, NCTP ) ! cloud top U + REAL(kind_phys) GCVT ( IJSDIM, NCTP ) ! cloud top V + REAL(kind_phys) GCLT ( IJSDIM, NCTP ) ! cloud top cloud water + REAL(kind_phys) GCIT ( IJSDIM, NCTP ) ! cloud top cloud ice REAL(kind_phys) GCtrT (IJSDIM, ntrq:ntr, NCTP) ! cloud top tracer - REAL(kind_phys) GTPRT ( IJSDIM, NCTP ) !! precipitation/M - REAL(kind_phys) GCLZ ( IJSDIM, KMAX ) !! cloud liquid for each CTP - REAL(kind_phys) GCIZ ( IJSDIM, KMAX ) !! cloud ice for each CTP + REAL(kind_phys) GTPRT ( IJSDIM, NCTP ) ! precipitation/M + REAL(kind_phys) GCLZ ( IJSDIM, KMAX ) ! cloud liquid for each CTP + REAL(kind_phys) GCIZ ( IJSDIM, KMAX ) ! cloud ice for each CTP - REAL(kind_phys) ACWF ( IJSDIM ) !! cloud work function - REAL(kind_phys) GPRCIZ( IJSDIM, KMAX+1, NCTP ) !! precipitation - REAL(kind_phys) GSNWIZ( IJSDIM, KMAX+1, NCTP ) !! snowfall - REAL(kind_phys) GTPRC0( IJSDIM ) !! precip. before evap. + REAL(kind_phys) ACWF ( IJSDIM ) ! cloud work function + REAL(kind_phys) GPRCIZ( IJSDIM, KMAX+1, NCTP ) ! precipitation + REAL(kind_phys) GSNWIZ( IJSDIM, KMAX+1, NCTP ) ! snowfall + REAL(kind_phys) GTPRC0( IJSDIM ) ! precip. before evap. - REAL(kind_phys) GMFLX ( IJSDIM, KMAX+1 ) !! mass flux (updraft+downdraft) - REAL(kind_phys) QLIQ ( IJSDIM, KMAX ) !! total cloud liquid - REAL(kind_phys) QICE ( IJSDIM, KMAX ) !! total cloud ice - REAL(kind_phys) GPRCI ( IJSDIM, KMAX ) !! rainfall generation - REAL(kind_phys) GSNWI ( IJSDIM, KMAX ) !! snowfall generation + REAL(kind_phys) GMFLX ( IJSDIM, KMAX+1 ) ! mass flux (updraft+downdraft) + REAL(kind_phys) QLIQ ( IJSDIM, KMAX ) ! total cloud liquid + REAL(kind_phys) QICE ( IJSDIM, KMAX ) ! total cloud ice + REAL(kind_phys) GPRCI ( IJSDIM, KMAX ) ! rainfall generation + REAL(kind_phys) GSNWI ( IJSDIM, KMAX ) ! snowfall generation - REAL(kind_phys) GPRCP ( IJSDIM, KMAX+1 ) !! rainfall flux + REAL(kind_phys) GPRCP ( IJSDIM, KMAX+1 ) ! rainfall flux ! - REAL(kind_phys) GTEVP ( IJSDIM, KMAX ) !! evaporation+sublimation - REAL(kind_phys) GMDD ( IJSDIM, KMAX+1 ) !! downdraft mass flux + REAL(kind_phys) GTEVP ( IJSDIM, KMAX ) ! evaporation+sublimation + REAL(kind_phys) GMDD ( IJSDIM, KMAX+1 ) ! downdraft mass flux - REAL(kind_phys) CUMHGT( IJSDIM, NCTP ) !! cloud top height - REAL(kind_phys) CTOPP ( IJSDIM ) !! cloud top pressure + REAL(kind_phys) CUMHGT( IJSDIM, NCTP ) ! cloud top height + REAL(kind_phys) CTOPP ( IJSDIM ) ! cloud top pressure - REAL(kind_phys) GDZTR ( IJSDIM ) !! tropopause height - REAL(kind_phys) FLIQOU( IJSDIM, KMAX ) !! liquid ratio in cumulus + REAL(kind_phys) GDZTR ( IJSDIM ) ! tropopause height + REAL(kind_phys) FLIQOU( IJSDIM, KMAX ) ! liquid ratio in cumulus !#ifdef OPT_CHASER ! REAL(kind_phys) TOPFLX( IJSDIM, NCTP ) !! flux at each cloud top !#endif INTEGER KB ( IJSDIM ) - INTEGER KSTRT ( IJSDIM ) !! tropopause level + INTEGER KSTRT ( IJSDIM ) ! tropopause level REAL(kind_phys) GAMX REAL(kind_phys) CIN ( IJSDIM ) INTEGER JBUOY ( IJSDIM ) REAL(kind_phys) DELZ, BUOY, DELWC, DELER -!M REAL(kind_phys) WCB ( NCTP ) !! updraft velocity**2 @base +!M REAL(kind_phys) WCB ( NCTP ) ! updraft velocity**2 @base !M SAVE WCB REAL(kind_phys) WCBX (IJSDIM) -! REAL(kind_phys) ERMR ( NCTP ) !! entrainment rate (ASMODE) +! REAL(kind_phys) ERMR ( NCTP ) ! entrainment rate (ASMODE) ! SAVE ERMR - INTEGER KTMX ( NCTP ) !! max of cloud top - INTEGER KTMXT !! max of cloud top + INTEGER KTMX ( NCTP ) ! max of cloud top + INTEGER KTMXT ! max of cloud top REAL(kind_phys) TIMED REAL(kind_phys) GDCLDX, GDMU2X, GDMU3X ! @@ -760,26 +760,26 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions INTEGER KBMX, I, K, CTP, ierr, n, kp1, l, l1, kk, kbi, kmi, km1 real(kind_phys) tem1, tem2, tem3, cbmfl, mflx_e, teme, tems - REAL(kind_phys) HBGT ( IJSDIM ) !! imbalance in column heat - REAL(kind_phys) WBGT ( IJSDIM ) !! imbalance in column water + REAL(kind_phys) HBGT ( IJSDIM ) ! imbalance in column heat + REAL(kind_phys) WBGT ( IJSDIM ) ! imbalance in column water !DDsigma begin local work variables - all on model interfaces (sfc=1) - REAL(kind_phys) lamdai( IJSDIM, KMAX+1, nctp ) !! lamda for cloud type ctp - REAL(kind_phys) lamdaprod( IJSDIM, KMAX+1 ) !! product of (1+lamda) through cloud type ctp - REAL(kind_phys) gdrhom !! density - REAL(kind_phys) gdtvm !! virtual temperature - REAL(kind_phys) gdqm, gdwm,gdlm, gdim !! water vaper + REAL(kind_phys) lamdai( IJSDIM, KMAX+1, nctp ) ! lamda for cloud type ctp + REAL(kind_phys) lamdaprod( IJSDIM, KMAX+1 ) ! product of (1+lamda) through cloud type ctp + REAL(kind_phys) gdrhom ! density + REAL(kind_phys) gdtvm ! virtual temperature + REAL(kind_phys) gdqm, gdwm,gdlm, gdim ! water vaper REAL(kind_phys) gdtrm(ntrq:ntr) ! tracer character(len=4) :: cproc !DDsigmadiag ! the following are new arguments to cumup to get them out - REAL(kind_phys) wcv( IJSDIM, KMAX+1, nctp) !! in-cloud vertical velocity - REAL(kind_phys) GCTM ( IJSDIM, KMAX+1 ) !! cloud T (half lev) !DDsigmadiag make output - REAL(kind_phys) GCQM ( IJSDIM, KMAX+1, nctp ) !! cloud q (half lev) !DDsigmadiag make output - REAL(kind_phys) GCwM ( IJSDIM, KMAX+1, nctp ) !! cloud q (half lev) !DDsigmadiag make output - REAL(kind_phys) GCiM ( IJSDIM, KMAX+1 ) !! cloud q (half lev) !DDsigmadiag make output - REAL(kind_phys) GClM ( IJSDIM, KMAX+1 ) !! cloud q (half lev) !DDsigmadiag make output - REAL(kind_phys) GChM ( IJSDIM, KMAX+1, nctp ) !! cloud q (half lev) !DDsigmadiag make output + REAL(kind_phys) wcv( IJSDIM, KMAX+1, nctp) ! in-cloud vertical velocity + REAL(kind_phys) GCTM ( IJSDIM, KMAX+1 ) ! cloud T (half lev) !DDsigmadiag make output + REAL(kind_phys) GCQM ( IJSDIM, KMAX+1, nctp ) ! cloud q (half lev) !DDsigmadiag make output + REAL(kind_phys) GCwM ( IJSDIM, KMAX+1, nctp ) ! cloud q (half lev) !DDsigmadiag make output + REAL(kind_phys) GCiM ( IJSDIM, KMAX+1 ) ! cloud q (half lev) !DDsigmadiag make output + REAL(kind_phys) GClM ( IJSDIM, KMAX+1 ) ! cloud q (half lev) !DDsigmadiag make output + REAL(kind_phys) GChM ( IJSDIM, KMAX+1, nctp ) ! cloud q (half lev) !DDsigmadiag make output REAL(kind_phys) GCtrM (IJSDIM, KMAX, ntrq:ntr) ! cloud tracer (half lev) !DDsigmadiag make output ! these are the fluxes at the interfaces - AW will operate on them @@ -800,30 +800,30 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions !DDsigma end local work variables ! ! [INTERNAL PARM] - REAL(kind_phys) :: WCBMIN = 0._kind_phys !! min. of updraft velocity at cloud base -!M REAL(kind_phys) :: WCBMAX = 1.4_kind_phys !! max. of updraft velocity at cloud base + REAL(kind_phys) :: WCBMIN = 0._kind_phys ! min. of updraft velocity at cloud base +!M REAL(kind_phys) :: WCBMAX = 1.4_kind_phys ! max. of updraft velocity at cloud base !M wcbas commented by Moorthi since it is not used -!M REAL(kind_phys) :: WCBAS = 2._kind_phys !! updraft velocity**2 at cloud base (ASMODE) -!M REAL(kind_phys) :: ERAMIN = 1.e-5_kind_phys !! min. of entrainment rate - !! used only in OPT_ASMODE -!M REAL(kind_phys) :: ERAMAX = 2.e-3_kind_phys !! max. of entrainment rate - !! used only in OPT_ASMODE +!M REAL(kind_phys) :: WCBAS = 2._kind_phys ! updraft velocity**2 at cloud base (ASMODE) +!M REAL(kind_phys) :: ERAMIN = 1.e-5_kind_phys ! min. of entrainment rate + ! used only in OPT_ASMODE +!M REAL(kind_phys) :: ERAMAX = 2.e-3_kind_phys ! max. of entrainment rate + ! used only in OPT_ASMODE ! downdraft mass flux terms now slot nctp+1 in the *fluxterm arrays - REAL(kind_phys) dtdwn ( IJSDIM, KMAX ) !! t tendency downdraft detrainment - REAL(kind_phys) dqvdwn ( IJSDIM, KMAX ) !! qv tendency downdraft detrainment - REAL(kind_phys) dqldwn ( IJSDIM, KMAX ) !! ql tendency downdraft detrainment - REAL(kind_phys) dqidwn ( IJSDIM, KMAX ) !! qi tendency downdraft detrainment + REAL(kind_phys) dtdwn ( IJSDIM, KMAX ) ! t tendency downdraft detrainment + REAL(kind_phys) dqvdwn ( IJSDIM, KMAX ) ! qv tendency downdraft detrainment + REAL(kind_phys) dqldwn ( IJSDIM, KMAX ) ! ql tendency downdraft detrainment + REAL(kind_phys) dqidwn ( IJSDIM, KMAX ) ! qi tendency downdraft detrainment REAL(kind_phys), dimension(ijsdim,kmax,ntrq:ntr) :: dtrdwn ! tracer tendency downdraft detrainment - LOGICAL :: OINICB = .false. !! set 0.d0 to CBMFX + LOGICAL :: OINICB = .false. ! set 0.d0 to CBMFX - REAL(kind_phys) :: VARMIN = 1.e-13_kind_phys !! minimum of PDF variance - REAL(kind_phys) :: VARMAX = 5.e-7_kind_phys !! maximum of PDF variance - REAL(kind_phys) :: SKWMAX = 0.566_kind_phys !! maximum of PDF skewness + REAL(kind_phys) :: VARMIN = 1.e-13_kind_phys ! minimum of PDF variance + REAL(kind_phys) :: VARMAX = 5.e-7_kind_phys ! maximum of PDF variance + REAL(kind_phys) :: SKWMAX = 0.566_kind_phys ! maximum of PDF skewness - REAL(kind_phys) :: PSTRMX = 400.e2_kind_phys !! max P of tropopause - REAL(kind_phys) :: PSTRMN = 50.e2_kind_phys !! min P of tropopause - REAL(kind_phys) :: GCRSTR = 1.e-4_kind_phys !! crit. dT/dz tropopause + REAL(kind_phys) :: PSTRMX = 400.e2_kind_phys ! max P of tropopause + REAL(kind_phys) :: PSTRMN = 50.e2_kind_phys ! min P of tropopause + REAL(kind_phys) :: GCRSTR = 1.e-4_kind_phys ! crit. dT/dz tropopause ! 0: mass fixer is not applied ! tracers which may become negative values @@ -835,10 +835,9 @@ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions real(kind=kind_phys), parameter :: zero=0.0, one=1.0 real(kind=kind_phys) :: tem, esat ! - LOGICAL, SAVE :: OFIRST = .TRUE. !! called first time? + LOGICAL, SAVE :: OFIRST = .TRUE. ! called first time? ! - IF ( OFIRST ) THEN - + IF (OFIRST) THEN OFIRST = .FALSE. IF (OINICB) THEN CBMFX = zero From a0801b69045ad6f6bf656dc451fcb6d85ccee298 Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Fri, 8 Mar 2024 15:20:43 +0000 Subject: [PATCH 15/21] uncomment calculation of cf_upi in order to allow calculation of w_upi --- physics/CONV/Chikira_Sugiyama/cs_conv.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/physics/CONV/Chikira_Sugiyama/cs_conv.F90 b/physics/CONV/Chikira_Sugiyama/cs_conv.F90 index ae59db7ba..4e7030dd5 100644 --- a/physics/CONV/Chikira_Sugiyama/cs_conv.F90 +++ b/physics/CONV/Chikira_Sugiyama/cs_conv.F90 @@ -466,7 +466,7 @@ subroutine cs_conv_run( IJSDIM , KMAX , ntracp1 , NN, & ! CNV_PRC3(i,k) = 0.0 CNV_NDROP(i,k) = 0.0 CNV_NICE(i,k) = 0.0 -! cf_upi(i,k) = max(0.0,min(0.01*log(1.0+500*ud_mf(i,k)),0.1)) + cf_upi(i,k) = max(0.0,min(0.01*log(1.0+500*ud_mf(i,k)),0.1)) ! CLCN(i,k) = cf_upi(i,k) !downdraft is below updraft !! clcn(i,k) = max(0.0,min(0.01*log(1.0+500*ud_mf(i,k)/delta),0.25)) @@ -500,7 +500,7 @@ subroutine cs_conv_run( IJSDIM , KMAX , ntracp1 , NN, & ! CNV_PRC3(i,k) = 0.0 CNV_NDROP(i,k) = 0.0 CNV_NICE(i,k) = 0.0 -! cf_upi(i,k) = max(0.0,min(0.01*log(1.0+500*ud_mf(i,k)),0.1)) + cf_upi(i,k) = max(0.0,min(0.01*log(1.0+500*ud_mf(i,k)),0.1)) ! & 500*ud_mf(i,k)),0.60)) ! CLCN(i,k) = cf_upi(i,k) !downdraft is below updraft From e1db7f20047f292fa7fe23f47b25fe2c535d1679 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Fri, 8 Mar 2024 09:49:44 -0700 Subject: [PATCH 16/21] In physics/Interstitials/UFS_SCM_NEPTUNE/maximum_hourly_diagnostics.meta: change units 'flashes 5 min-1' to 'flashes min-1' and update long name to make clear this is per 5 minutes --- .../UFS_SCM_NEPTUNE/maximum_hourly_diagnostics.meta | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/maximum_hourly_diagnostics.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/maximum_hourly_diagnostics.meta index 0c2d1bcbe..5e9b97333 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/maximum_hourly_diagnostics.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/maximum_hourly_diagnostics.meta @@ -295,24 +295,24 @@ intent = in [ltg1_max] standard_name = lightning_threat_index_1 - long_name = lightning threat index 1 - units = flashes 5 min-1 + long_name = lightning threat index 1 in flashes per 5 minutes + units = flashes min-1 dimensions = (horizontal_loop_extent) type = real kind = kind_phys intent = inout [ltg2_max] standard_name = lightning_threat_index_2 - long_name = lightning threat index 2 - units = flashes 5 min-1 + long_name = lightning threat index 2 in flashes per 5 minutes + units = flashes min-1 dimensions = (horizontal_loop_extent) type = real kind = kind_phys intent = inout [ltg3_max] standard_name = lightning_threat_index_3 - long_name = lightning threat index 3 - units = flashes 5 min-1 + long_name = lightning threat index 3 in flashes per 5 minutes + units = flashes min-1 dimensions = (horizontal_loop_extent) type = real kind = kind_phys From 26f95141ce3ef4cbc6d78dc9c98eda8dca48466b Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Fri, 8 Mar 2024 12:40:51 -0700 Subject: [PATCH 17/21] In physics/Interstitials/UFS_SCM_NEPTUNE/maximum_hourly_diagnostics.F90, scale lightning threat from flashes per 5 minutes to flashes per minute to match units in metadata --- .../UFS_SCM_NEPTUNE/maximum_hourly_diagnostics.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/maximum_hourly_diagnostics.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/maximum_hourly_diagnostics.F90 index cd1016053..f66aa77e9 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/maximum_hourly_diagnostics.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/maximum_hourly_diagnostics.F90 @@ -80,6 +80,12 @@ subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics, !Lightning threat indices if (lightning_threat) then call lightning_threat_indices + ! Lightning threat indices are calculated as flashes per 5 minutes. + ! In order to scale that to flashes per minute (standard units), + ! we must divide the indices by a factor of 5 / multiply by 0.2 + ltg1_max = 0.2_kind_phys * ltg1_max + ltg2_max = 0.2_kind_phys * ltg2_max + ltg3_max = 0.2_kind_phys * ltg3_max endif !Calculate hourly max 1-km agl and -10C reflectivity From 743dc85aa132a41502c678671a2de93517da8cf8 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Mon, 11 Mar 2024 16:25:16 +0000 Subject: [PATCH 18/21] correctly convert from flashes per five minutes to flashes per minute --- .../maximum_hourly_diagnostics.F90 | 21 ++++++++++++------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/maximum_hourly_diagnostics.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/maximum_hourly_diagnostics.F90 index f66aa77e9..5ac28afe8 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/maximum_hourly_diagnostics.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/maximum_hourly_diagnostics.F90 @@ -15,6 +15,9 @@ module maximum_hourly_diagnostics real(kind=kind_phys), parameter ::PQ0=379.90516E0, A2A=17.2693882, A3=273.16, A4=35.86, RHmin=1.0E-6 ! *DH + ! Conversion from flashes per five minutes to flashes per minute. + real(kind=kind_phys), parameter :: scaling_factor = 0.2 + contains #if 0 @@ -80,12 +83,6 @@ subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics, !Lightning threat indices if (lightning_threat) then call lightning_threat_indices - ! Lightning threat indices are calculated as flashes per 5 minutes. - ! In order to scale that to flashes per minute (standard units), - ! we must divide the indices by a factor of 5 / multiply by 0.2 - ltg1_max = 0.2_kind_phys * ltg1_max - ltg2_max = 0.2_kind_phys * ltg2_max - ltg3_max = 0.2_kind_phys * ltg3_max endif !Calculate hourly max 1-km agl and -10C reflectivity @@ -201,7 +198,10 @@ subroutine lightning_threat_indices endif IF ( ltg1 .LT. clim1 ) ltg1 = 0. - + + ! Scale to flashes per minue + ltg1 = ltg1 * scaling_factor + IF ( ltg1 .GT. ltg1_max(i) ) THEN ltg1_max(i) = ltg1 ENDIF @@ -214,14 +214,19 @@ subroutine lightning_threat_indices ltg2 = coef2 * totice_colint(i) IF ( ltg2 .LT. clim2 ) ltg2 = 0. + + ! Scale to flashes per minute + ltg2 = ltg2 * scaling_factor IF ( ltg2 .GT. ltg2_max(i) ) THEN ltg2_max(i) = ltg2 ENDIF + ! This calculation is already in flashes per minute. ltg3_max(i) = 0.95 * ltg1_max(i) + 0.05 * ltg2_max(i) - IF ( ltg3_max(i) .LT. clim3 ) ltg3_max(i) = 0. + ! Thus, we must scale clim3. The compiler will optimize this away. + IF ( ltg3_max(i) .LT. clim3 * scaling_factor ) ltg3_max(i) = 0. enddo end subroutine lightning_threat_indices From 7e74ada4d3a1fcbdeab272aa81a8af7b2b878300 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Mon, 11 Mar 2024 16:39:18 +0000 Subject: [PATCH 19/21] correct the meta file --- .../UFS_SCM_NEPTUNE/maximum_hourly_diagnostics.meta | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/maximum_hourly_diagnostics.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/maximum_hourly_diagnostics.meta index 5e9b97333..5d18bd9bd 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/maximum_hourly_diagnostics.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/maximum_hourly_diagnostics.meta @@ -295,7 +295,7 @@ intent = in [ltg1_max] standard_name = lightning_threat_index_1 - long_name = lightning threat index 1 in flashes per 5 minutes + long_name = lightning threat index 1 units = flashes min-1 dimensions = (horizontal_loop_extent) type = real @@ -303,7 +303,7 @@ intent = inout [ltg2_max] standard_name = lightning_threat_index_2 - long_name = lightning threat index 2 in flashes per 5 minutes + long_name = lightning threat index 2 units = flashes min-1 dimensions = (horizontal_loop_extent) type = real @@ -311,7 +311,7 @@ intent = inout [ltg3_max] standard_name = lightning_threat_index_3 - long_name = lightning threat index 3 in flashes per 5 minutes + long_name = lightning threat index 3 units = flashes min-1 dimensions = (horizontal_loop_extent) type = real From c77b9e8f4772450b5c7ec89575ccb53aa5bcaa22 Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Wed, 13 Mar 2024 10:11:08 -0400 Subject: [PATCH 20/21] fix metadata error in cs_conv.meta --- physics/CONV/Chikira_Sugiyama/cs_conv.meta | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/CONV/Chikira_Sugiyama/cs_conv.meta b/physics/CONV/Chikira_Sugiyama/cs_conv.meta index 49e460ed6..d75ab1006 100644 --- a/physics/CONV/Chikira_Sugiyama/cs_conv.meta +++ b/physics/CONV/Chikira_Sugiyama/cs_conv.meta @@ -258,7 +258,7 @@ standard_name = convective_updraft_area_fraction_at_model_interfaces long_name = convective updraft area fraction at model interfaces units = frac - dimensions = (horizontal_loop_extent,vertical_layer_dimension) + dimensions = (horizontal_loop_extent,vertical_interface_dimension) type = real kind = kind_phys intent = out From bbdec2ffceeb0cda70aaaa4f44e619a5468e7be3 Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Wed, 13 Mar 2024 10:50:38 -0400 Subject: [PATCH 21/21] fix convective_updraft_area_fraction_at_model_interfaces metadata --- physics/CONV/Chikira_Sugiyama/cs_conv_post.meta | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/CONV/Chikira_Sugiyama/cs_conv_post.meta b/physics/CONV/Chikira_Sugiyama/cs_conv_post.meta index 75de3fca7..5877c051b 100644 --- a/physics/CONV/Chikira_Sugiyama/cs_conv_post.meta +++ b/physics/CONV/Chikira_Sugiyama/cs_conv_post.meta @@ -33,7 +33,7 @@ standard_name = convective_updraft_area_fraction_at_model_interfaces long_name = convective updraft area fraction at model interfaces units = frac - dimensions = (horizontal_loop_extent,vertical_layer_dimension) + dimensions = (horizontal_loop_extent,vertical_interface_dimension) type = real kind = kind_phys intent = in