diff --git a/sorc/chgres_cube.fd/atmosphere.F90 b/sorc/chgres_cube.fd/atmosphere.F90 index 12d117fd4..f3feba3a2 100644 --- a/sorc/chgres_cube.fd/atmosphere.F90 +++ b/sorc/chgres_cube.fd/atmosphere.F90 @@ -319,8 +319,10 @@ subroutine atmosphere_driver(localpet) call vintg if( wam_cold_start ) then + if (lev_input.lt.149) then call vintg_wam (cycle_year,cycle_mon,cycle_day,cycle_hour) endif + endif !----------------------------------------------------------------------------------- ! Compute height. @@ -807,7 +809,7 @@ subroutine newpr1(localpet) real(esmf_kind_r8), allocatable :: pi(:,:,:) print*,"COMPUTE 3-D PRESSURE FROM ADJUSTED SURFACE PRESSURE." - + ! for WAM idvc = 3, but vcoord(:,3)=0.0 idvc = 2 ! hard wire for now. idsl = 2 ! hard wire for now. @@ -2147,7 +2149,7 @@ subroutine compute_zh enddo enddo - + print *,'zhptr',minval(zhptr),maxval(zhptr) deallocate(pe0, pn0) end subroutine compute_zh diff --git a/sorc/chgres_cube.fd/input_data.F90 b/sorc/chgres_cube.fd/input_data.F90 index 193d295de..c79352952 100644 --- a/sorc/chgres_cube.fd/input_data.F90 +++ b/sorc/chgres_cube.fd/input_data.F90 @@ -103,6 +103,15 @@ module input_data !! for 7 layers of soil for the RUC LSM character(len=50), private, allocatable :: slevs(:) !< The atmospheric levels in the GRIB2 input file. + !---WAM + real(esmf_kind_r8), allocatable, public :: ri(:),cpi(:) !R & cp for multi gases + type(esmf_field) :: psx_input_grid ! ps gradient + type(esmf_field) :: psy_input_grid + type(esmf_field) :: div_input_grid ! divergence + type(esmf_field) :: pm_input_grid ! for diag omega + type(esmf_field) :: pd_input_grid ! for diag omega + type(esmf_field) :: dpm_input_grid ! for diag omega + type(esmf_field) :: dpd_input_grid ! for diag omega ! Fields associated with the nst model. @@ -138,6 +147,7 @@ module input_data public :: convert_winds public :: init_sfc_esmf_fields public :: dint2p + public :: get_omega contains @@ -531,7 +541,69 @@ subroutine init_atm_esmf_fields call error_handler("IN FieldCreate", rc) end subroutine init_atm_esmf_fields +!--diag omega field init + subroutine init_atm_omega_esmf_fields + implicit none + integer :: rc + print*,"- CALL FieldCreate FOR INPUT GRID SURFACE PRESSURE GRADIENT." + psx_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + psy_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + print*,"- CALL FieldCreate FOR INPUT GRID DIVERGENCE." + div_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + ungriddedLBound=(/1/), & + ungriddedUBound=(/lev_input/), rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + print*,"- CALL FieldCreate FOR INPUT PM." + pm_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + ungriddedLBound=(/1/), & + ungriddedUBound=(/lev_input/), rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + print*,"- CALL FieldCreate FOR INPUT PD." + pd_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + ungriddedLBound=(/1/), & + ungriddedUBound=(/lev_input/), rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + print*,"- CALL FieldCreate FOR INPUT PM." + dpm_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + ungriddedLBound=(/1/), & + ungriddedUBound=(/lev_input/), rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + print*,"- CALL FieldCreate FOR INPUT PM." + dpd_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + ungriddedLBound=(/1/), & + ungriddedUBound=(/lev_input/), rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + end subroutine init_atm_omega_esmf_fields +!--- !> Create surface input grid esmf fields !! !! @author George Gayno NCEP/EMC @@ -746,7 +818,7 @@ end subroutine init_sfc_esmf_fields subroutine read_input_atm_gfs_sigio_file(localpet) use sigio_module - + use mpi implicit none integer, intent(in) :: localpet @@ -767,6 +839,16 @@ subroutine read_input_atm_gfs_sigio_file(localpet) type(sigio_head) :: sighead type(sigio_dbta) :: sigdata + !WAM + real(esmf_kind_r8) :: cp0= 1039.6 !?? N2 + real(esmf_kind_r8) :: psmean + integer :: thermodyn_id + integer :: nvcoord + real, allocatable :: vcoord(:,:) + real(esmf_kind_r8), allocatable :: dummy2d2(:,:) + real(esmf_kind_r8), allocatable :: dummy(:,:,:),dummy3d_dpd(:,:,:),dummy3d_dpm(:,:,:) + real(esmf_kind_r8), pointer :: tptr(:,:,:), qptr(:,:,:), wptr(:,:,:) + integer:: im,idvc,idsl,idvm the_file = trim(data_dir_input_grid) // "/" // trim(atm_files_input_grid(1)) print*,"- ATMOSPHERIC DATA IN SIGIO FORMAT." @@ -786,6 +868,25 @@ subroutine read_input_atm_gfs_sigio_file(localpet) lev_input = sighead%levs levp1_input = lev_input + 1 + nvcoord = sighead%nvcoord + allocate(vcoord(levp1_input,nvcoord)) + vcoord = sighead%vcoord + idvc = sighead%idvc + idsl = sighead%idsl + idvm = sighead%idvm + im = i_input*j_input + !---WAM + thermodyn_id = mod(idvm/10,10) + allocate(ri(0:sighead%ntrac)) + allocate(cpi(0:sighead%ntrac)) + + if (thermodyn_id == 3) then + ri = sighead%ri + cpi = sighead%cpi + else + ri = 0.0 + cpi = 0.0 + endif if (num_tracers_input /= sighead%ntrac) then call error_handler("WRONG NUMBER OF TRACERS EXPECTED.", 99) endif @@ -796,6 +897,14 @@ subroutine read_input_atm_gfs_sigio_file(localpet) trim(tracers_input(3)) /= 'clwmr') then call error_handler("TRACERS SELECTED DO NOT MATCH FILE CONTENTS.", 99) endif + elseif(sighead%idvt == 200) then ! WAM 'spfh,,o3mr,clwmr,o,o2' + if (trim(tracers_input(1)) /= 'spfh' .or. & + trim(tracers_input(2)) /= 'o3mr' .or. & + trim(tracers_input(3)) /= 'clwmr' .or. & + trim(tracers_input(4)) /= 'omr' .or. & + trim(tracers_input(5)) /= 'o2mr') then + call error_handler("TRACERS SELECTED DO NOT MATCH FILE CONTENTS.", 99) + end if else print*,'- UNRECOGNIZED IDVT: ', sighead%idvt call error_handler("UNRECOGNIZED IDVT", 99) @@ -807,14 +916,24 @@ subroutine read_input_atm_gfs_sigio_file(localpet) call init_atm_esmf_fields + call init_atm_omega_esmf_fields + if (localpet == 0) then allocate(dummy2d(i_input,j_input)) + allocate(dummy2d2(i_input,j_input)) + allocate(dummy(i_input,j_input,lev_input)) allocate(dummy3d(i_input,j_input,lev_input)) allocate(dummy3d2(i_input,j_input,lev_input)) + allocate(dummy3d_dpd(i_input,j_input,lev_input)) + allocate(dummy3d_dpm(i_input,j_input,lev_input)) else allocate(dummy2d(0,0)) + allocate(dummy2d2(0,0)) + allocate(dummy(0,0,0)) allocate(dummy3d(0,0,0)) allocate(dummy3d2(0,0,0)) + allocate(dummy3d_dpd(0,0,0)) + allocate(dummy3d_dpm(0,0,0)) endif if (localpet == 0) then @@ -829,11 +948,19 @@ subroutine read_input_atm_gfs_sigio_file(localpet) call error_handler("READING SIGDATA.", rc) endif call sptez(0,sighead%jcap,4,i_input, j_input, sigdata%ps, dummy2d, 1) + select case(mod(idvm,10)) + case(0,1) dummy2d = exp(dummy2d) * 1000.0 + case(2) + dummy2d = dummy2d* 1000.0 + case default + print *,' default selected: psi is p in pascal ' + end select + ! dummy2d = exp(dummy2d) * 1000.0 print*,'surface pres ',maxval(dummy2d),minval(dummy2d) - endif + endif ! localpet == 0 - print*,"- CALL FieldScatter FOR SURFACE PRESSURE." + if (localpet == 0) print*,"- CALL FieldScatter FOR SURFACE PRESSURE." call ESMF_FieldScatter(ps_input_grid, dummy2d, rootpet=0, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldScatter", rc) @@ -843,19 +970,34 @@ subroutine read_input_atm_gfs_sigio_file(localpet) print*,'terrain ',maxval(dummy2d),minval(dummy2d) endif - print*,"- CALL FieldScatter FOR TERRAIN." + if (localpet == 0) print*,"- CALL FieldScatter FOR TERRAIN." call ESMF_FieldScatter(terrain_input_grid, dummy2d, rootpet=0, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldScatter", rc) - + if (localpet == 0) then + if (thermodyn_id == 3) then + dummy3d2 = 0.0 !sumq + dummy = 0.0 !xcp + endif + endif do k = 1, num_tracers_input if (localpet == 0) then call sptezm(0,sighead%jcap,4,i_input, j_input, lev_input, sigdata%q(:,:,k), dummy3d, 1) print*,trim(tracers_input(k)),maxval(dummy3d),minval(dummy3d) + if (thermodyn_id == 3) then + if( cpi(k) .ne. 0.0 .and. ri(k) .ne. 0.0) then + dummy = dummy + cpi(k)*dummy3d !xcp + dummy3d2 = dummy3d2 + dummy3d !sumq + endif + else + if (k == 1) then + dummy = (1.+(461.50/287.05-1)*dummy3d) !virt + endif + endif endif - print*,"- CALL FieldScatter FOR INPUT ", trim(tracers_input(k)) + if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT ", trim(tracers_input(k)) call ESMF_FieldScatter(tracers_input_grid(k), dummy3d, rootpet=0, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldScatter", rc) @@ -864,16 +1006,30 @@ subroutine read_input_atm_gfs_sigio_file(localpet) if (localpet == 0) then call sptezm(0,sighead%jcap,4,i_input, j_input, lev_input, sigdata%t, dummy3d, 1) + print *,"-: WAM ENTHAPHY==>TEMPERATURE" + select case( thermodyn_id ) + case(0,1) + dummy3d = dummy3d/dummy + case(2) + case(3) + if (cpi(0) == 0) then + dummy3d = dummy3d/cp0 + else + dummy = (1.-dummy3d2)*cpi(0)+dummy + dummy3d = dummy3d/dummy + endif + case default + end select print*,'temp ',maxval(dummy3d),minval(dummy3d) endif - print*,"- CALL FieldScatter FOR INPUT GRID TEMPERATURE." + if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT GRID TEMPERATURE." call ESMF_FieldScatter(temp_input_grid, dummy3d, rootpet=0, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldScatter", rc) !--------------------------------------------------------------------------- -! The spectral gfs files have omega, not vertical velocity. Set to +! The spectral gfs files have omega?? not vertical velocity. Set to ! zero for now. Convert from omega to vv in the future? !--------------------------------------------------------------------------- @@ -882,7 +1038,7 @@ subroutine read_input_atm_gfs_sigio_file(localpet) dummy3d = 0.0 endif - print*,"- CALL FieldScatter FOR INPUT DZDT." + if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT DZDT." call ESMF_FieldScatter(dzdt_input_grid, dummy3d, rootpet=0, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldScatter", rc) @@ -893,22 +1049,102 @@ subroutine read_input_atm_gfs_sigio_file(localpet) print*,'v ',maxval(dummy3d2),minval(dummy3d2) endif - print*,"- CALL FieldScatter FOR INPUT U-WIND." + if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT U-WIND." call ESMF_FieldScatter(u_input_grid, dummy3d, rootpet=0, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldScatter", rc) - print*,"- CALL FieldScatter FOR INPUT V-WIND." + if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT V-WIND." call ESMF_FieldScatter(v_input_grid, dummy3d2, rootpet=0, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldScatter", rc) +!--------------------------------------------------------------------------- +! DIV-divergence here. idrt = 4 !GAUSSIAN GRID +!--------------------------------------------------------------------------- + if (localpet == 0) then + call sptezm(0,sighead%jcap,4,i_input, j_input, lev_input, sigdata%d, dummy3d, 1) + print*,'div ',maxval(dummy3d),minval(dummy3d) + endif - deallocate(dummy2d, dummy3d, dummy3d2) + if (localpet == 0) print*,"- CALL FieldScatter FOR INPUT DIV." + call ESMF_FieldScatter(div_input_grid, dummy3d, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) +!--------------------------------------------------------------------------- +! ps gradient here. idrt = 4 !GAUSSIAN GRID +!--------------------------------------------------------------------------- + if (localpet == 0) then + call sptezd(0,sighead%jcap,4,i_input,j_input,sigdata%ps,psmean,dummy2d,dummy2d2,1) + print *,'psx,psy',maxval(dummy2d),minval(dummy2d),maxval(dummy2d2),minval(dummy2d2) + endif + + if (localpet == 0) print*,"- CALL FieldScatter FOR SURFACE PRESSURE GRADIENT-X." + call ESMF_FieldScatter(psx_input_grid, dummy2d, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + if (localpet == 0) print*,"- CALL FieldScatter FOR SURFACE PRESSURE GRADIENT-Y." + call ESMF_FieldScatter(psy_input_grid, dummy2d2, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + deallocate(dummy2d2) if (localpet == 0) call sigio_axdbta(sigdata, iret) call sigio_sclose(21, iret) +!--------------------------------------------------------------------------- +! OMEGA RELATED +!--------------------------------------------------------------------------- + if (localpet == 0) print*,"- CALL FieldGather TEMPERATURE." + call ESMF_FieldGather(temp_input_grid,dummy,rootPet=0, tile=1, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGether", rc) + + if (localpet == 0) print*,"- CALL FieldGather SURFACE PRESSURE." + call ESMF_FieldGather(ps_input_grid,dummy2d,rootPet=0, tile=1, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGether", rc) + + if (localpet==0) then + call sigio_modprd(im,im,lev_input,nvcoord,idvc,idsl,vcoord,iret, & + ps=dummy2d,t=dummy,pm=dummy3d,pd=dummy3d2,dpmdps=dummy3d_dpm,dpddps=dummy3d_dpd) + print *,'modprd',minval(dummy3d),maxval(dummy3d) + endif + if (localpet == 0) print*,"- CALL FieldScatter MODPRD-PD." + call ESMF_FieldScatter(pd_input_grid, dummy3d2, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + + if (localpet == 0) print*,"- CALL FieldScatter MODPRD-PM." + call ESMF_FieldScatter(pm_input_grid, dummy3d, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + + if (localpet == 0) print*,"- CALL FieldScatter MODPRD-DPD." + call ESMF_FieldScatter(dpd_input_grid, dummy3d_dpd, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + + if (localpet == 0) print*,"- CALL FieldScatter MODPRD-DPM." + call ESMF_FieldScatter(dpm_input_grid, dummy3d_dpm, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + + deallocate(dummy2d, dummy3d, dummy3d2) + deallocate(dummy,dummy3d_dpd, dummy3d_dpm) + deallocate(vcoord) + +!--- WAM GET OMEGA +if (localpet == 0) print *,'GET OMEGA HERE' +nullify(wptr) +if (localpet == 0) print*,"- CALL FieldGet DZDT." +call ESMF_FieldGet(dzdt_input_grid, & + farrayPtr=wptr, rc=rc) +if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) +call get_omega(wptr,idvm,localpet) +print *,'OMEGA',minval(wptr),maxval(wptr) +call cleanup_input_atm_omega_data !--------------------------------------------------------------------------- ! Convert from 2-d to 3-d component winds. !--------------------------------------------------------------------------- @@ -975,6 +1211,19 @@ subroutine read_input_atm_gfs_sigio_file(localpet) print*,'pres ',psptr(clb(1),clb(2)),pptr(clb(1),clb(2),:) endif + nullify(tptr) + if (localpet == 0) print*,"- CALL FieldGet TEMPERATURE." + call ESMF_FieldGet(temp_input_grid, & + farrayPtr=tptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + if (localpet == 0) print*,"- CONVERT FROM OMEGA TO DZDT." + nullify(qptr) + call ESMF_FieldGet(tracers_input_grid(1), & + farrayPtr=qptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) end subroutine read_input_atm_gfs_sigio_file !> Read input atmospheric data from spectral gfs (global gaussian in @@ -7330,6 +7579,23 @@ subroutine cleanup_input_atm_data deallocate(tracers_input_grid) end subroutine cleanup_input_atm_data +!> WAM-OMEGA related fields +!! + subroutine cleanup_input_atm_omega_data + implicit none + + integer :: rc + print*,'- DESTROY OMEGA RELATED INPUT DATA.' + + call ESMF_FieldDestroy(psx_input_grid, rc=rc) + call ESMF_FieldDestroy(psy_input_grid, rc=rc) + call ESMF_FieldDestroy(div_input_grid, rc=rc) + call ESMF_FieldDestroy(pm_input_grid, rc=rc) + call ESMF_FieldDestroy(pd_input_grid, rc=rc) + call ESMF_FieldDestroy(dpd_input_grid, rc=rc) + call ESMF_FieldDestroy(dpm_input_grid, rc=rc) + end subroutine cleanup_input_atm_omega_data + !> Free up memory associated with nst data. !! @@ -7726,5 +7992,151 @@ SUBROUTINE DINT2P(PPIN,XXIN,NPIN,PPOUT,XXOUT,NPOUT & RETURN END SUBROUTINE DINT2P + subroutine get_omega(wptr,idvm,localpet) + use mpi + implicit none + + integer,intent(in) :: idvm,localpet + real(esmf_kind_r8), pointer :: wptr(:,:,:) + real(esmf_kind_r8), pointer :: psptr(:,:), psxptr(:,:), psyptr(:,:), & + uptr(:,:,:), vptr(:,:,:), div(:,:,:) + real(esmf_kind_r8), pointer :: pd(:,:,:),pm(:,:,:), & + dpmdps(:,:,:),dpddps(:,:,:) + integer :: clb(3), cub(3) + real,allocatable :: psx(:,:), psy(:,:) + real,allocatable :: pi(:,:,:), dpidps(:,:,:) + real,allocatable :: w(:,:,:) + + real :: vgradp,os + integer :: i, j , k, rc, iret + + if (localpet == 0) print*,"- CALL FieldGet FOR U-WIND" + nullify(uptr) + call ESMF_FieldGet(u_input_grid, & + computationalLBound=clb, & + computationalUBound=cub, & + farrayPtr=uptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + if (localpet == 0) print*,"- CALL FieldGet FOR V-WIND" + nullify(vptr) + call ESMF_FieldGet(v_input_grid, & + farrayPtr=vptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + if (localpet == 0) print*,"- CALL FieldGet FOR SURFACE PRESSURE." + nullify(psptr) + call ESMF_FieldGet(ps_input_grid, & + farrayPtr=psptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + if (localpet == 0) print*,"- CALL FieldGet FOR PS GRADIENT-X" + nullify(psxptr) + call ESMF_FieldGet(psx_input_grid, & + farrayPtr=psxptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + if (localpet == 0) print*,"- CALL FieldGet FOR PS GRADIENT-Y" + nullify(psyptr) + call ESMF_FieldGet(psy_input_grid, & + farrayPtr=psyptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + if (localpet == 0) print*,"- CALL FieldGet FOR DIV" + nullify(div) + call ESMF_FieldGet(div_input_grid, & + farrayPtr=div, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + if (localpet == 0) print*,"- CALL FieldGet FOR MODPRD-PM" + nullify(pm) + call ESMF_FieldGet(pm_input_grid, & + farrayPtr=pm, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + if (localpet == 0) print*,"- CALL FieldGet FOR MODPRD-PD" + nullify(pd) + call ESMF_FieldGet(pd_input_grid, & + farrayPtr=pd, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + if (localpet == 0) print*,"- CALL FieldGet FOR MODPRD-DPD" + nullify(dpddps) + call ESMF_FieldGet(dpd_input_grid, & + farrayPtr=dpddps, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + if (localpet == 0) print*,"- CALL FieldGet FOR MODPRD-DPM" + nullify(dpmdps) + call ESMF_FieldGet(dpm_input_grid, & + farrayPtr=dpmdps, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldGet", rc) + + allocate(psx(clb(1):cub(1),clb(2):cub(2))) + allocate(psy(clb(1):cub(1),clb(2):cub(2))) + allocate(pi(clb(1):cub(1),clb(2):cub(2),1:levp1_input)) + allocate(dpidps(clb(1):cub(1),clb(2):cub(2),1:levp1_input)) + allocate(w(clb(1):cub(1),clb(2):cub(2),1:lev_input)) + + select case(mod(idvm,10)) + case(0,1) + continue + case(2) +!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j) + do i = clb(1), cub(1) + do j = clb(2), cub(2) + psx(i,j) = psxptr(i,j)/(psptr(i,j)*1.0e-3) + psy(i,j) = psyptr(i,j)/(psptr(i,j)*1.0e-3) + enddo + enddo +!$OMP END PARALLEL DO + case default +!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j) + do i = clb(1), cub(1) + do j = clb(2), cub(2) + psx(i,j) = psxptr(i,j)/psptr(i,j) + psy(i,j) = psxptr(i,j)/psptr(i,j) + enddo + enddo +!$OMP END PARALLEL DO + end select + +!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j) + do i = clb(1), cub(1) + do j = clb(2), cub(2) + pi(i,j,1) = psptr(i,j) + dpidps(i,j,1) = 1. + do k = 1,lev_input + pi(i,j,k+1) = pi(i,j,k)-pd(i,j,k) + dpidps(i,j,k+1) = dpidps(i,j,k)-dpddps(i,j,k) + enddo + os = 0.0 + do k = lev_input, 1, -1 + vgradp = uptr(i,j,k)*psx(i,j)+vptr(i,j,k)*psy(i,j) + os = os-vgradp*psptr(i,j)*(dpmdps(i,j,k)-dpidps(i,j,k+1))- & + div(i,j,k)*(pm(i,j,k)-pi(i,j,k+1)) + w(i,j,k) = vgradp*psptr(i,j)*dpmdps(i,j,k)+os + os = os-vgradp*psptr(i,j)*(dpidps(i,j,k)-dpmdps(i,j,k))- & + div(i,j,k)*(pi(i,j,k)-pm(i,j,k)) + wptr(i,j,k) = w(i,j,k) + enddo + enddo + enddo +!$OMP END PARALLEL DO + deallocate(psx,psy) + deallocate(pi,dpidps) + deallocate(w) + + end subroutine get_omega end module input_data diff --git a/util/gdas_init/config_wam b/util/gdas_init/config_wam new file mode 100644 index 000000000..d10a83482 --- /dev/null +++ b/util/gdas_init/config_wam @@ -0,0 +1,115 @@ +#----------------------------------------------------------- +# +# 1) Compile the chgres_cube program. Invoke +# ./sorc/build_chgres_cube.sh +# +# 2) Ensure links to the 'fixed' directories are +# set. See the ./sorc/link_fixdirs.sh script prolog +# for details. +# +# 3) Set all config variables. See definitions +# below. +# +# 4) Invoke the driver script for your machine (with no +# arguments). +# +# Variable definitions: +# -------------------- +# EXTRACT_DIR - Directory where data extracted from HPSS +# is stored. +# EXTRACT_DATA - Set to 'yes' to extract data from HPSS. +# If data has been extracted and is located +# in EXTRACT_DIR, set to 'no'. +# RUN_CHGRES - To run chgres, set to 'yes'. To extract +# data only, set to 'no'. +# yy/mm/dd/hh - The year/month/day/hour of your desired +# experiment. Currently, does not support +# pre-ENKF GFS data, prior to +# 2012 May 21 00z. Use two digits. +# LEVS - Number of hybrid levels plus 1. To +# run with 64 levels, set LEVS to 65. +# CRES_HIRES - Resolution of the hires component of +# your experiment. +# CRES_ENKF - Resolution of the enkf component of the +# your experiment. +# UFS_DIR - Location of your checked out UFS_UTILS +# repo. +# OUTDIR - Directory where the coldstart data output +# from chgres is stored. +# CDUMP - When 'gdas', will process gdas and enkf +# members. When 'gfs', will process gfs +# member for running free forecast only. +# use_v16retro - When 'yes', use v16 retro parallel data. +# The retro parallel tarballs can be missing +# or incomplete. So this option may not +# always work. Contact george.gayno@noaa.gov +# if you encounter problems. +# +#----------------------------------------------------------- + +# EXTRACT_DIR=/lfs/h2/emc/stmp/$USER/gdas.init/input +# EXTRACT_DATA=no +EXTRACT_DIR=/scratch1/NCEPDEV/stmp2/$USER/chgres/input +EXTRACT_DATA=no + +RUN_CHGRES=yes + +yy=2018 +mm=02 +dd=01 +hh=00 + +use_v16retro=no + +LEVS=150 + +CDUMP=wdas + +CRES_HIRES=C96 +CRES_ENKF=C96 + +UFS_DIR=$PWD/../.. + +OUTDIR=/scratch1/NCEPDEV/stmp2/$USER/chgres/output +if [[ "$CDUMP" = "wfs" ]] || [[ "$CDUMP" = "wdas" ]]; then + gfs_ver=wam +else +#--------------------------------------------------------- +# Dont touch anything below here. +#--------------------------------------------------------- + +if [ "$use_v16retro" = "yes" ]; then + + gfs_ver=v16retro + +else + + gfs_ver=v16 + +# No ENKF data prior to 2012/05/21/00z + if [ $yy$mm$dd$hh -lt 2012052100 ]; then + set +x + echo FATAL ERROR: SCRIPTS DO NOT SUPPORT OLD GFS DATA + exit 2 + elif [ $yy$mm$dd$hh -lt 2016051000 ]; then + gfs_ver=v12 + elif [ $yy$mm$dd$hh -lt 2017072000 ]; then + gfs_ver=v13 + elif [ $yy$mm$dd$hh -lt 2019061200 ]; then + gfs_ver=v14 + elif [ $yy$mm$dd$hh -lt 2021032100 ]; then + gfs_ver=v15 +# The way the v16 switch over was done, there is no complete +# set of v16 or v15 data for 2021032100. And although +# v16 was officially implemented 2021032212, the v16 prod +# tarballs were archived starting 2021032106. + elif [ $yy$mm$dd$hh -lt 2021032106 ]; then + set +x + echo FATAL ERROR: NO V15 OR V16 DATA FOR 2021032100 + exit 1 + fi + +fi +fi +export EXTRACT_DIR yy mm dd hh UFS_DIR OUTDIR CRES_HIRES CRES_ENKF +export LEVS gfs_ver diff --git a/util/gdas_init/driver.hera.sh b/util/gdas_init/driver.hera.sh index 2694b7b18..bdf3fa3df 100755 --- a/util/gdas_init/driver.hera.sh +++ b/util/gdas_init/driver.hera.sh @@ -1,4 +1,4 @@ -#!/bin/bash + #!/bin/bash #--------------------------------------------------------------------- # Driver script for running on Hera. @@ -7,7 +7,7 @@ #--------------------------------------------------------------------- set -x - +ulimit -s unlimited compiler=${compiler:-"intel"} source ../../sorc/machine-setup.sh > /dev/null 2>&1 module use ../../modulefiles @@ -19,9 +19,8 @@ module use -a /scratch2/NCEPDEV/nwprod/NCEPLIBS/modulefiles module load prod_util/1.1.0 PROJECT_CODE=fv3-cpu -QUEUE=batch - -source config +QUEUE=debug +source config_wam if [ $EXTRACT_DATA == yes ]; then @@ -108,6 +107,9 @@ if [ $EXTRACT_DATA == yes ]; then DEPEND="-d afterok:$DATAH:$DATA1:$DATA2:$DATA3:$DATA4:$DATA5:$DATA6:$DATA7:$DATA8" fi ;; + wam) + echo "TTBD" + ;; esac else # do not extract data. @@ -162,6 +164,10 @@ if [ $RUN_CHGRES == yes ]; then sbatch --parsable --ntasks-per-node=6 --nodes=${NODES} -t $WALLT -A $PROJECT_CODE -q $QUEUE -J chgres_${CDUMP} \ -o log.${CDUMP} -e log.${CDUMP} ${DEPEND} run_v16.chgres.sh ${CDUMP} ;; + wam) + sbatch --parsable --ntasks-per-node=6 --nodes=${NODES} -t $WALLT -A $PROJECT_CODE -q $QUEUE -J chgres_${CDUMP} \ + -o log.${CDUMP}.$$.txt -e log.${CDUMP}.$$.txt ${DEPEND} run_wam.chgres.sh ${CDUMP} + ;; esac if [ "$CDUMP" = "gdas" ]; then diff --git a/util/gdas_init/run_wam.chgres.sh b/util/gdas_init/run_wam.chgres.sh new file mode 100755 index 000000000..ed378008a --- /dev/null +++ b/util/gdas_init/run_wam.chgres.sh @@ -0,0 +1,112 @@ +#!/bin/bash + +copy_data() +{ + +mkdir -p $SAVEDIR +cp gfs_ctrl.nc $SAVEDIR + +for tile in 'tile1' 'tile2' 'tile3' 'tile4' 'tile5' 'tile6' +do + cp out.atm.${tile}.nc ${SAVEDIR}/gfs_data.${tile}.nc + cp out.sfc.${tile}.nc ${SAVEDIR}/sfc_data.${tile}.nc +done +} + +#--------------------------------------------------------------------------- +# Run chgres using v16 netcdf history data as input. These history +# files are part of the OPS v16 gfs/gdas/enkf tarballs, and the +# v16 retro parallel gfs tarballs. To run using the v16 retro +# gdas tarballs (which contain warm restart files), the +# run_v16retro.chgres.sh is used. +#--------------------------------------------------------------------------- + +set -x + +MEMBER=$1 + +FIX_FV3=$UFS_DIR/fix +FIX_ORO=${FIX_FV3}/orog +FIX_AM=${FIX_FV3}/am + +WORKDIR=${WORKDIR:-$OUTDIR/work.${MEMBER}} + +if [ ${MEMBER} == 'wdas' ] || [ ${MEMBER} == 'wfs' ] ; then + CTAR=${CRES_HIRES} +#--------------------------------------------------------------------------- +# Some gfs tarballs from the v16 retro parallels dont have 'atmos' +# in their path. Account for this. +#--------------------------------------------------------------------------- + INPUT_DATA_DIR="${EXTRACT_DIR}/${MEMBER}.${yy}${mm}${dd}/${hh}/atmos" + if [ ! -d ${INPUT_DATA_DIR} ]; then + INPUT_DATA_DIR="${EXTRACT_DIR}/${MEMBER}.${yy}${mm}${dd}/${hh}" + fi + ATMFILE="${MEMBER}.t${hh}z.atmanl" + SFCFILE="${MEMBER}.t${hh}z.sfcanl" +else + date10=`$NDATE -6 $yy$mm$dd$hh` + yy_d=$(echo $date10 | cut -c1-4) + mm_d=$(echo $date10 | cut -c5-6) + dd_d=$(echo $date10 | cut -c7-8) + hh_d=$(echo $date10 | cut -c9-10) + CTAR=${CRES_ENKF} + INPUT_DATA_DIR="${EXTRACT_DIR}/enkfgdas.${yy_d}${mm_d}${dd_d}/${hh_d}/atmos/mem${MEMBER}" + ATMFILE="gdas.t${hh_d}z.atmf006.nc" + SFCFILE="gdas.t${hh_d}z.sfcf006.nc" +fi + +rm -fr $WORKDIR +mkdir -p $WORKDIR +cd $WORKDIR + +cat << EOF > fort.41 + +&config + fix_dir_target_grid="${FIX_ORO}/${CTAR}/fix_sfc" + mosaic_file_target_grid="${FIX_ORO}/${CTAR}/${CTAR}_mosaic.nc" + orog_dir_target_grid="${FIX_ORO}/${CTAR}" + orog_files_target_grid="${CTAR}_oro_data.tile1.nc","${CTAR}_oro_data.tile2.nc","${CTAR}_oro_data.tile3.nc","${CTAR}_oro_data.tile4.nc","${CTAR}_oro_data.tile5.nc","${CTAR}_oro_data.tile6.nc" + data_dir_input_grid="${INPUT_DATA_DIR}" + atm_files_input_grid="${ATMFILE}" + sfc_files_input_grid="${SFCFILE}" + vcoord_file_target_grid="${FIX_AM}/global_hyblev.l${LEVS}.txt" + cycle_mon=$mm + cycle_day=$dd + cycle_hour=$hh + convert_atm=.true. + convert_sfc=.true. + convert_nst=.false. + wam_cold_start=.true. + input_type="gfs_sigio" + tracers="sphum","spo3","liq_wat","spo","spo2" + tracers_input="spfh","o3mr","clwmr","omr","o2mr" +/ +EOF + +$APRUN $UFS_DIR/exec/chgres_cube +rc=$? + +if [ $rc != 0 ]; then + exit $rc +fi + +if [ ${MEMBER} == 'gdas' ] || [ ${MEMBER} == 'gfs' ]; then + SAVEDIR=$OUTDIR/${MEMBER}.${yy}${mm}${dd}/${hh}/atmos/INPUT + copy_data + touch $SAVEDIR/../${MEMBER}.t${hh}z.loginc.txt + if [ ${MEMBER} == 'gdas' ]; then + cp ${INPUT_DATA_DIR}/*abias* $SAVEDIR/.. + cp ${INPUT_DATA_DIR}/*radstat $SAVEDIR/.. + fi +else + SAVEDIR=$OUTDIR/${yy}${mm}${dd}/${hh}/${MEMBER}/INPUT + copy_data + touch $SAVEDIR/../t${hh}z.loginc.txt +fi + +rm -fr $WORKDIR + +set +x +echo CHGRES COMPLETED FOR MEMBER $MEMBER + +exit 0