From e198c5252ac222be5f5d1a3b1fea138f70fe4421 Mon Sep 17 00:00:00 2001 From: "Jun.Wang" Date: Tue, 5 Nov 2019 15:51:23 +0000 Subject: [PATCH] this PR includes: 1)vlab #69814: post update,2)ufs issue #3, Add 3D reflectivity to restart file and restart reproducibility for regional fv3,3)ufs issue #5, Updates to WW3, 4)vlab #69735, update netcdf time units attribute when iau_offset --- .gitmodules | 4 +-- atmos_cubed_sphere | 2 +- gfsphysics/GFS_layer/GFS_physics_driver.F90 | 2 +- gfsphysics/GFS_layer/GFS_restart.F90 | 11 +++++- gfsphysics/GFS_layer/GFS_typedefs.F90 | 12 ++++++- io/FV3GFS_io.F90 | 9 ++--- io/module_write_netcdf.F90 | 39 ++++++++++++++++++--- io/module_wrt_grid_comp.F90 | 4 +-- io/post_gfs.F90 | 32 +++++++++++++---- 9 files changed, 90 insertions(+), 25 deletions(-) diff --git a/.gitmodules b/.gitmodules index 949d298df..04f5ebec0 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,7 +1,7 @@ [submodule "atmos_cubed_sphere"] path = atmos_cubed_sphere - url = https://github.com/NOAA-EMC/GFDL_atmos_cubed_sphere - branch = dev/emc + url = https://github.com/junwang-noaa/GFDL_atmos_cubed_sphere + branch = regfv3_rst [submodule "ccpp/framework"] path = ccpp/framework url = https://github.com/NCAR/ccpp-framework diff --git a/atmos_cubed_sphere b/atmos_cubed_sphere index 786447c83..cc28aa5aa 160000 --- a/atmos_cubed_sphere +++ b/atmos_cubed_sphere @@ -1 +1 @@ -Subproject commit 786447c8391a6806cd7b869bfa9dca69e3c95a48 +Subproject commit cc28aa5aa94817b24aea033cc19c1792e5927459 diff --git a/gfsphysics/GFS_layer/GFS_physics_driver.F90 b/gfsphysics/GFS_layer/GFS_physics_driver.F90 index 401fbbf86..5b67f7faa 100644 --- a/gfsphysics/GFS_layer/GFS_physics_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_physics_driver.F90 @@ -1149,7 +1149,7 @@ subroutine GFS_physics_driver & if (fice(i) < one) then wet(i) = .true. ! Sfcprop%tsfco(i) = tgice - Sfcprop%tsfco(i) = max(Sfcprop%tisfc(i), tgice) + if (.not. Model%cplflx) Sfcprop%tsfco(i) = max(Sfcprop%tisfc(i), tgice) ! Sfcprop%tsfco(i) = max((Sfcprop%tsfc(i) - fice(i)*sfcprop%tisfc(i)) & ! / (one - fice(i)), tgice) endif diff --git a/gfsphysics/GFS_layer/GFS_restart.F90 b/gfsphysics/GFS_layer/GFS_restart.F90 index eafbcb9ba..a24cc0fc6 100644 --- a/gfsphysics/GFS_layer/GFS_restart.F90 +++ b/gfsphysics/GFS_layer/GFS_restart.F90 @@ -117,6 +117,9 @@ subroutine GFS_restart_populate (Restart, Model, Statein, Stateout, Sfcprop, & #endif Restart%num3d = Model%ntot3d + if(Model%lrefres) then + Restart%num3d = Model%ntot3d+1 + endif #ifdef CCPP ! GF if (Model%imfdeepcnv == 3) then @@ -252,7 +255,13 @@ subroutine GFS_restart_populate (Restart, Model, Statein, Stateout, Sfcprop, & Restart%data(nb,num)%var3p => Tbd(nb)%phy_f3d(:,:,num) enddo enddo - + if (Model%lrefres) then + num = Model%ntot3d+1 + restart%name3d(num) = 'ref_f3d' + do nb = 1,nblks + Restart%data(nb,num)%var3p => IntDiag(nb)%refl_10cm(:,:) + enddo + endif #ifdef CCPP !--- RAP/HRRR-specific variables, 3D num = Model%ntot3d diff --git a/gfsphysics/GFS_layer/GFS_typedefs.F90 b/gfsphysics/GFS_layer/GFS_typedefs.F90 index 7f8239a5a..39520b0d4 100644 --- a/gfsphysics/GFS_layer/GFS_typedefs.F90 +++ b/gfsphysics/GFS_layer/GFS_typedefs.F90 @@ -701,6 +701,9 @@ module GFS_typedefs !--- GFDL microphysical paramters logical :: lgfdlmprad !< flag for GFDL mp scheme and radiation consistency + !--- Thompson,GFDL mp parameter + logical :: lrefres !< flag for radar reflectivity in restart file + !--- land/surface model parameters integer :: lsm !< flag for land surface model lsm=1 for noah lsm integer :: lsm_noah=1 !< flag for NOAH land surface model @@ -2740,6 +2743,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & !--- GFDL microphysical parameters logical :: lgfdlmprad = .false. !< flag for GFDLMP radiation interaction + !--- Thompson,GFDL microphysical parameter + logical :: lrefres = .false. !< flag for radar reflectivity in restart file + !--- land/surface model parameters integer :: lsm = 1 !< flag for land surface model to use =0 for osu lsm; =1 for noah lsm; =2 for noah mp lsm; =3 for RUC lsm integer :: lsoil = 4 !< number of soil layers @@ -3023,7 +3029,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & mg_do_graupel, mg_do_hail, mg_nccons, mg_nicons, mg_ngcons, & mg_ncnst, mg_ninst, mg_ngnst, sed_supersat, do_sb_physics, & mg_alf, mg_qcmin, mg_do_ice_gmao, mg_do_liq_liu, & - ltaerosol, lradar, ttendlim, lgfdlmprad, & + ltaerosol, lradar, lrefres, ttendlim, lgfdlmprad, & !--- max hourly avg_max_length, & !--- land/surface model control @@ -3312,6 +3318,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%ttendlim = ttendlim !--- gfdl MP parameters Model%lgfdlmprad = lgfdlmprad +!--- Thompson,GFDL MP parameter + Model%lrefres = lrefres !--- land/surface model parameters Model%lsm = lsm @@ -4310,6 +4318,7 @@ subroutine control_print(Model) print *, ' Thompson microphysical parameters' print *, ' ltaerosol : ', Model%ltaerosol print *, ' lradar : ', Model%lradar + print *, ' lrefres : ', Model%lrefres print *, ' ttendlim : ', Model%ttendlim print *, ' ' endif @@ -4327,6 +4336,7 @@ subroutine control_print(Model) if (Model%imp_physics == Model%imp_physics_gfdl) then print *, ' GFDL microphysical parameters' print *, ' GFDL MP radiation inter: ', Model%lgfdlmprad + print *, ' lrefres : ', Model%lrefres print *, ' ' endif diff --git a/io/FV3GFS_io.F90 b/io/FV3GFS_io.F90 index 3e7b7d2e7..e7f7dbd57 100644 --- a/io/FV3GFS_io.F90 +++ b/io/FV3GFS_io.F90 @@ -1119,15 +1119,10 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain) Sfcprop(nb)%zorll(ix) = Sfcprop(nb)%zorlo(ix) Sfcprop(nb)%zorl(ix) = Sfcprop(nb)%zorlo(ix) Sfcprop(nb)%tsfc(ix) = Sfcprop(nb)%tsfco(ix) - if (Sfcprop(nb)%slmsk(ix) < 0.1 .or. Sfcprop(nb)%slmsk(ix) > 1.9) then + if (Sfcprop(nb)%slmsk(ix) > 1.9) then Sfcprop(nb)%landfrac(ix) = 0.0 - if (Sfcprop(nb)%oro_uf(ix) > 0.01) then - Sfcprop(nb)%lakefrac(ix) = 1.0 ! lake - else - Sfcprop(nb)%lakefrac(ix) = 0.0 ! ocean - endif else - Sfcprop(nb)%landfrac(ix) = 1.0 ! land + Sfcprop(nb)%landfrac(ix) = Sfcprop(nb)%slmsk(ix) endif enddo enddo diff --git a/io/module_write_netcdf.F90 b/io/module_write_netcdf.F90 index 1fce3d8b9..5de1362af 100644 --- a/io/module_write_netcdf.F90 +++ b/io/module_write_netcdf.F90 @@ -22,7 +22,7 @@ module module_write_netcdf contains !---------------------------------------------------------------------------------------- - subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc) + subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, idate, rc) ! type(ESMF_FieldBundle), intent(in) :: fieldbundle type(ESMF_FieldBundle), intent(in) :: wrtfb @@ -30,6 +30,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc integer, intent(in) :: mpi_comm integer, intent(in) :: mype integer, intent(in) :: im, jm + integer, intent(in) :: idate(7) integer, optional,intent(out) :: rc ! !** local vars @@ -146,7 +147,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc call add_dim(ncid, "phalf", phalf_dimid, wrtgrid, rc) end if - call add_dim(ncid, "time", time_dimid, wrtgrid, rc) + call add_dim(ncid, "time", time_dimid, wrtgrid, rc, idate=idate) call get_global_attr(wrtfb, ncid, rc) @@ -503,11 +504,12 @@ subroutine get_grid_attr(grid, prefix, ncid, varid, rc) end subroutine get_grid_attr - subroutine add_dim(ncid, dim_name, dimid, grid, rc) + subroutine add_dim(ncid, dim_name, dimid, grid, idate, rc) integer, intent(in) :: ncid character(len=*), intent(in) :: dim_name - integer, intent(inout) :: dimid + integer, intent(inout) :: dimid type(ESMF_Grid), intent(in) :: grid + integer, intent(in), optional :: idate(7) integer, intent(out) :: rc ! local variable @@ -556,9 +558,38 @@ subroutine add_dim(ncid, dim_name, dimid, grid, rc) call get_grid_attr(grid, dim_name, ncid, dim_varid, rc) + ! if write grid comp changes time units + if ( present (idate) ) then + if ( trim(dim_name) == "time") then + ncerr = nf90_get_att(ncid, dim_varid, 'units', time_units); NC_ERR_STOP(ncerr) + time_units = get_time_units_from_idate(idate) + ncerr = nf90_put_att(ncid, dim_varid, 'units', trim(time_units)); NC_ERR_STOP(ncerr) + endif + endif + end subroutine add_dim ! !---------------------------------------------------------------------------------------- + function get_time_units_from_idate(idate, time_measure) result(time_units) + ! create time units attribute of form 'hours since YYYY-MM-DD HH:MM:SS' + ! from integer array with year,month,day,hour,minute,second + ! optional argument 'time_measure' can be used to change 'hours' to + ! 'days', 'minutes', 'seconds' etc. + character(len=*), intent(in), optional :: time_measure + integer, intent(in) :: idate(7) + character(len=12) :: timechar + character(len=nf90_max_name) :: time_units + if (present(time_measure)) then + timechar = trim(time_measure) + else + timechar = 'hours' + endif + write(time_units,101) idate(1:6) +101 format(' since ',i4.4,'-',i2.2,'-',i2.2,' ',& + i2.2,':',i2.2,':',i2.2) + time_units = trim(adjustl(timechar))//time_units + end function get_time_units_from_idate +! subroutine nccheck(status) use netcdf implicit none diff --git a/io/module_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90 index d0846be53..f8546614b 100644 --- a/io/module_wrt_grid_comp.F90 +++ b/io/module_wrt_grid_comp.F90 @@ -1386,7 +1386,7 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) wbeg = MPI_Wtime() call write_netcdf(file_bundle,wrt_int_state%wrtFB(nbdl),trim(filename), & - wrt_mpi_comm,wrt_int_state%mype,imo,jmo,rc) + wrt_mpi_comm,wrt_int_state%mype,imo,jmo,idate,rc) wend = MPI_Wtime() if (lprnt) then write(*,'(A,F10.5,A,I4.2,A,I2.2)')' netcdf Write Time is ',wend-wbeg & @@ -1431,7 +1431,7 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) wbeg = MPI_Wtime() call write_netcdf(file_bundle,wrt_int_state%wrtFB(nbdl),trim(filename), & - wrt_mpi_comm,wrt_int_state%mype,imo,jmo,rc) + wrt_mpi_comm,wrt_int_state%mype,imo,jmo,idate,rc) wend = MPI_Wtime() if (mype == lead_write_task) then write(*,'(A,F10.5,A,I4.2,A,I2.2)')' netcdf Write Time is ',wend-wbeg & diff --git a/io/post_gfs.F90 b/io/post_gfs.F90 index 7b5a87026..ce6bb281e 100644 --- a/io/post_gfs.F90 +++ b/io/post_gfs.F90 @@ -451,8 +451,7 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & end do ! ! GFS does not output PD -! pt = 10000. ! this is for 100 hPa added by Moorthi - pt = 0. + pt = ak5(1) ! GFS may not have model derived radar ref. ! TKE @@ -1296,7 +1295,7 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & endif ! inst incoming sfc longwave - if(trim(fieldname)=='dlwsf') then + if(trim(fieldname)=='dlwrf') then !$omp parallel do private(i,j) do j=jsta,jend do i=ista, iend @@ -1845,6 +1844,16 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & ! time averaged surface clear sky outgoing LW if(trim(fieldname)=='csulf') then + !$omp parallel do private(i,j) + do j=jsta,jend + do i=ista, iend + alwoutc(i,j) = arrayr42d(i,j) + enddo + enddo + endif + + ! time averaged TOA clear sky outgoing LW + if(trim(fieldname)=='csulftoa') then !$omp parallel do private(i,j) do j=jsta,jend do i=ista, iend @@ -1864,7 +1873,7 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & endif ! time averaged TOA clear sky outgoing SW - if(trim(fieldname)=='csusf') then + if(trim(fieldname)=='csusftoa') then !$omp parallel do private(i,j) do j=jsta,jend do i=ista, iend @@ -2271,7 +2280,6 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & enddo end do -!??? reset pint(lev=1) !$omp parallel do private(i,j) do j=jsta,jend do i=1,im @@ -2282,7 +2290,7 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & ! print *,'in setvar, pt=',pt,'ak5(lp1)=', ak5(lp1),'ak5(1)=',ak5(1) ! compute alpint - do l=lp1,2,-1 + do l=lp1,1,-1 !$omp parallel do private(i,j) do j=jsta,jend do i=1,im @@ -2321,6 +2329,18 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & enddo enddo +! compute cwm for gfdlmp + if( imp_physics == 11 ) then + do l=1,lm +!$omp parallel do private(i,j) + do j=jsta,jend + do i=ista,iend + cwm(i,j,l)=qqg(i,j,l)+qqs(i,j,l)+qqr(i,j,l)+qqi(i,j,l)+qqw(i,j,l) + enddo + enddo + enddo + endif + ! estimate 2m pres and convert t2m to theta !$omp parallel do private(i,j) do j=jsta,jend