[WAM-IPE] WAM-IPE r94274: Commit Houjun and Phil IAU code
Samuel.Trahan at noaa.gov
Samuel.Trahan at noaa.gov
Fri Jun 16 16:56:09 UTC 2017
Friendly WAM-IPE developers,
This is an automated email about a WAM-IPE commit.
Project: WAM-IPE
URL: https://svnemc.ncep.noaa.gov/projects/gsm/branches/WAM-IPE/quasitrunk
Revision: 94274
Author: weiyu.yang at noaa.gov
Date: 2017-06-12T18:58:09.842413Z
Message:
Commit Houjun and Phil IAU code
See attached file for full differences.
First 4000 bytes of differences:
Index: checkout/phys/do_physics_one_step.f
===================================================================
--- checkout/phys/do_physics_one_step.f (revision 93160)
+++ checkout/phys/do_physics_one_step.f (revision 94274)
@@ -53,6 +53,7 @@
!! Feb 25 2016 S Moorthi - add kdt_dif and associated chnages
!! March 2016 Hang Lei Add DDT mdl_parm for physics driver
!! Mar 24 2016 Xingren Wu Connect Ice/Snow thickness (volume)
+!! Aug 17 2016 P Pegion add logic to know if in iauinterval
!!
USE machine, ONLY: KIND_GRID, KIND_RAD, kind_phys
use resol_def
@@ -89,6 +90,7 @@
TYPE(G3D_Var_Data) :: g3d_fld
TYPE(G2D_Var_Data) :: g2d_fld
type(aoi_var_data) :: aoi_fld
+ logical :: iniauinterval
! type(gfs_phy_tracer_type) :: gfs_phy_tracer
!* REAL(KIND=KIND_GRID) GRID_GR(lonr*lats_node_r_max,lotgr)
CHARACTER(16) :: CFHOUR1
@@ -160,6 +162,7 @@
! &,' ldfi=',ldfi,' me=',me,' kdt=',kdt,' kdt_start=',kdt_start
lssav = .true.
+ iniauinterval = .false.
if(ndfi > 0 .and. kdt_dif > ndfi/2 .and.
& kdt_dif <= ndfi .and. ldfi ) then
lssav = .false.
@@ -560,6 +563,7 @@
!
!!
!* call gloopb ( grid_gr,
+ if (iau) call checkiauforcing(phour*3600,iniauinterval)
call gloopb (grid_fld, g3d_fld, sfc_fld,
& flx_fld, aoi_fld, nst_fld,
& lats_nodes_r, global_lats_r, lonsperlar,
@@ -570,7 +574,7 @@
& ozplin, jindx1, jindx2, ddy,
& phy_f3d, phy_f2d, phy_fctd, nctp,
& xlat, nblck, kdt, restart_step,
- & mdl_parm)
+ & mdl_parm, iniauinterval)
!
!!
endif ! if (comp_task) then
Index: checkout/phys/compns_physics.f
===================================================================
--- checkout/phys/compns_physics.f (revision 93160)
+++ checkout/phys/compns_physics.f (revision 94274)
@@ -169,7 +169,7 @@
& sfcio_out,climate,a2oi_out,cplflx,ngrid_a2oi,griboro,
& sppt,skeb,shum,vc,vcamp,
! iau parameters
- & iau,iaufiles_fg,iaufiles_anl,iaufhrs,iau_delthrs
+ & iau,iaufiles_fg,iaufiles_anl,iaufhrs,iau_delthrs,iaulnp
!
shuff_lats_r = .false.
! reshuff_lats_r = .false.
Index: checkout/phys/makefile
===================================================================
--- checkout/phys/makefile (revision 93160)
+++ checkout/phys/makefile (revision 94274)
@@ -84,6 +84,7 @@
wrt3d.o \
aoicpl_prep.o \
para_fixio_w.o \
+ checkiauforcing.o \
para_nst_w.o
#
@@ -175,6 +176,9 @@
para_nstio_w.o: para_nstio_w.f90
$(FC) $(FFLAGS1) -c para_nstio_w.f90
+checkiauforcing.o: checkiauforcing.f90
+ $(FC) $(FFLAG90) -c checkiauforcing.f90
+
wrtout_physics.o: wrtout_physics.f90
$(FC) $(FFLAGS1) -c wrtout_physics.f90
Index: checkout/phys/namelist_physics_def.f
===================================================================
--- checkout/phys/namelist_physics_def.f (revision 93160)
+++ checkout/phys/namelist_physics_def.f (revision 94274)
@@ -39,6 +39,7 @@
! iau parameters
logical :: iau = .false. ! iau forcing included
+ integer :: iaulnp=2 ! 0: no pressure update, 1: log pressure, 2: pressure
integer :: iau_delthrs = 6 ! iau time interval (to scale increments)
character(len=120), dimension(7) :: iaufiles_fg,iaufiles_anl
real(kind=kind_evod), dimension(7) :: iaufhrs
Index: checkout/phys/wrtout_physics.f
===================================================================
--- checkout/phys/wrtout_physics.f (revision 93160)
+++ checkout/phys/wrtout_physics.f (revision 94274)
@@ -86,7 +86,7 @@
& nodes_comp
use namelist_physics_def, ONLY: gen
... see attachment for the rest ...
-------------- next part --------------
Index: checkout/phys/do_physics_one_step.f
===================================================================
--- checkout/phys/do_physics_one_step.f (revision 93160)
+++ checkout/phys/do_physics_one_step.f (revision 94274)
@@ -53,6 +53,7 @@
!! Feb 25 2016 S Moorthi - add kdt_dif and associated chnages
!! March 2016 Hang Lei Add DDT mdl_parm for physics driver
!! Mar 24 2016 Xingren Wu Connect Ice/Snow thickness (volume)
+!! Aug 17 2016 P Pegion add logic to know if in iauinterval
!!
USE machine, ONLY: KIND_GRID, KIND_RAD, kind_phys
use resol_def
@@ -89,6 +90,7 @@
TYPE(G3D_Var_Data) :: g3d_fld
TYPE(G2D_Var_Data) :: g2d_fld
type(aoi_var_data) :: aoi_fld
+ logical :: iniauinterval
! type(gfs_phy_tracer_type) :: gfs_phy_tracer
!* REAL(KIND=KIND_GRID) GRID_GR(lonr*lats_node_r_max,lotgr)
CHARACTER(16) :: CFHOUR1
@@ -160,6 +162,7 @@
! &,' ldfi=',ldfi,' me=',me,' kdt=',kdt,' kdt_start=',kdt_start
lssav = .true.
+ iniauinterval = .false.
if(ndfi > 0 .and. kdt_dif > ndfi/2 .and.
& kdt_dif <= ndfi .and. ldfi ) then
lssav = .false.
@@ -560,6 +563,7 @@
!
!!
!* call gloopb ( grid_gr,
+ if (iau) call checkiauforcing(phour*3600,iniauinterval)
call gloopb (grid_fld, g3d_fld, sfc_fld,
& flx_fld, aoi_fld, nst_fld,
& lats_nodes_r, global_lats_r, lonsperlar,
@@ -570,7 +574,7 @@
& ozplin, jindx1, jindx2, ddy,
& phy_f3d, phy_f2d, phy_fctd, nctp,
& xlat, nblck, kdt, restart_step,
- & mdl_parm)
+ & mdl_parm, iniauinterval)
!
!!
endif ! if (comp_task) then
Index: checkout/phys/compns_physics.f
===================================================================
--- checkout/phys/compns_physics.f (revision 93160)
+++ checkout/phys/compns_physics.f (revision 94274)
@@ -169,7 +169,7 @@
& sfcio_out,climate,a2oi_out,cplflx,ngrid_a2oi,griboro,
& sppt,skeb,shum,vc,vcamp,
! iau parameters
- & iau,iaufiles_fg,iaufiles_anl,iaufhrs,iau_delthrs
+ & iau,iaufiles_fg,iaufiles_anl,iaufhrs,iau_delthrs,iaulnp
!
shuff_lats_r = .false.
! reshuff_lats_r = .false.
Index: checkout/phys/makefile
===================================================================
--- checkout/phys/makefile (revision 93160)
+++ checkout/phys/makefile (revision 94274)
@@ -84,6 +84,7 @@
wrt3d.o \
aoicpl_prep.o \
para_fixio_w.o \
+ checkiauforcing.o \
para_nst_w.o
#
@@ -175,6 +176,9 @@
para_nstio_w.o: para_nstio_w.f90
$(FC) $(FFLAGS1) -c para_nstio_w.f90
+checkiauforcing.o: checkiauforcing.f90
+ $(FC) $(FFLAG90) -c checkiauforcing.f90
+
wrtout_physics.o: wrtout_physics.f90
$(FC) $(FFLAGS1) -c wrtout_physics.f90
Index: checkout/phys/namelist_physics_def.f
===================================================================
--- checkout/phys/namelist_physics_def.f (revision 93160)
+++ checkout/phys/namelist_physics_def.f (revision 94274)
@@ -39,6 +39,7 @@
! iau parameters
logical :: iau = .false. ! iau forcing included
+ integer :: iaulnp=2 ! 0: no pressure update, 1: log pressure, 2: pressure
integer :: iau_delthrs = 6 ! iau time interval (to scale increments)
character(len=120), dimension(7) :: iaufiles_fg,iaufiles_anl
real(kind=kind_evod), dimension(7) :: iaufhrs
Index: checkout/phys/wrtout_physics.f
===================================================================
--- checkout/phys/wrtout_physics.f (revision 93160)
+++ checkout/phys/wrtout_physics.f (revision 94274)
@@ -86,7 +86,7 @@
& nodes_comp
use namelist_physics_def, ONLY: gen_coord_hybrid, ldiag3d,
& hybrid, fhlwr, fhswr, ens_nam,
- & nstf_name, sfcio_out
+ & nstf_name, sfcio_out,iau
use mpi_def, ONLY: liope, info, mpi_comm_all,
& mc_comp, mpi_comm_null
use gfs_physics_sfc_flx_mod, ONLY: Sfc_Var_Data, Flx_Var_Data
@@ -125,6 +125,10 @@
!!
character CFHOUR*40,CFORM*40
integer jdate(4),ndigyr,ndig,kh,IOPROC
+! --- locals for iau:
+ integer :: jdate_iau(4),idat(8),jdat(8)
+ real (kind=kind_evod) :: zhour_iau,fhour_iau
+ real :: rinc(5)
!!
real(kind=kind_io8) slmsk(lonr,lats_node_r)
!!
@@ -154,11 +158,32 @@
IF(NDIGYR == 2) THEN
JDATE(4) = MOD(IDATE(4)-1,100) + 1
ENDIF
+ if ( iau .and. fhour >= 6.0 ) then
+ idat = 0
+ idat(1) = idate(4)
+ idat(2) = idate(2)
+ idat(3) = idate(3)
+ idat(5) = idate(1)
+ rinc = 0.
+ rinc(2) = 6.
+ call w3movdat(rinc,idat,jdat)
+ jdate_iau(4) = jdat(1)
+ jdate_iau(2) = jdat(2)
+ jdate_iau(3) = jdat(3)
+ jdate_iau(1) = jdat(5)
+ fhour_iau = fhour - 6.0
+ zhour_iau = zhour - 6.0
+ if ( zhour_iau < 0.0 ) zhour_iau = 0.0
+ else
+ jdate_iau = jdate
+ fhour_iau = fhour
+ zhour_iau = zhour
+ endif
!sela set lfnhr to false for writing one step output etc.
lfnhr = .true. ! no output
! lfnhr = 3600*abs(fhour-nint(fhour)).le.1.or.phour.eq.0
- lfnhr = 3600*abs(fhour-nint(fhour)) <= 1
+ lfnhr = 3600*abs(fhour_iau-nint(fhour_iau)) <= 1
IF (LFNHR) THEN
KH = NINT(FHOUR)
NDIG = MAX(LOG10(KH+0.5)+1.,2.)
@@ -165,7 +190,7 @@
WRITE(CFORM,'("(I",I1,".",I1,")")') NDIG,NDIG
WRITE(CFHOUR,CFORM) KH
ELSE
- KS = NINT(FHOUR*3600)
+ KS = NINT(fhour_iau*3600)
KH = KS/3600
KM = (KS-KH*3600)/60
KS = KS-KH*3600-KM*60
@@ -195,7 +220,7 @@
t4 = rtc()
tbd = t4-t3
t3 = rtc()
- SECPHY = (FHOUR-ZHOUR)*3600.
+ SECPHY = (fhour_iau-zhour_iau)*3600.
SECSWR = MAX(SECPHY,FHSWR)
SECLWR = MAX(SECPHY,FHLWR)
@@ -223,13 +248,14 @@
! predating the I/O task updates.
!
call wrtflx_a
- & (IOPROC,noflx,ZHOUR,FHOUR,IDATE,colat1,SECSWR,SECLWR,
+ & (IOPROC,noflx,zhour_iau,fhour_iau,jdate_iau,colat1,
+ & SECSWR,SECLWR,
& sfc_fld, flx_fld, fluxr, global_lats_r,lonsperlar,
& slmsk)
if ( ngrids_aer > 0) then
call wrtaer
- & (IOPROC,noaer,ZHOUR,FHOUR,IDATE,
+ & (IOPROC,noaer,ZHOUR_iau,FHOUR_iau,JDATE_iau,
& sfc_fld, g2d_fld, global_lats_r, lonsperlar)
endif
@@ -246,12 +272,12 @@
if (me == 0) print *,'---- start sfc.f section -----'
call sfc_only_move(ioproc)
cosfc = sfcf//CFHOUR
- call sfc_wrt(ioproc,nosfc,cosfc,fhour,jdate
+ call sfc_wrt(ioproc,nosfc,cosfc,fhour_iau,jdate_iau
&, global_lats_r,lonsperlar)
if ( nstf_name(1) > 0 ) then
call nst_only_move(ioproc)
const = nstf//CFHOUR
- call nst_wrt(ioproc,nonst,const,fhour,jdate
+ call nst_wrt(ioproc,nonst,const,fhour_iau,jdate_iau
&, global_lats_r,lonsperlar)
endif
call flx_only_move(ioproc)
@@ -265,8 +291,8 @@
call baopenwt(noflx,cosfc,iostat)
! print *,'after open flx file,',trim(cosfc)
endif
- call wrtflx_w(IOPROC,noflx,ZHOUR,FHOUR,IDATE,colat1,SECSWR,
- & SECLWR,slmsk, global_lats_r,lonsperlar)
+ call wrtflx_w(IOPROC,noflx,zhour_iau,fhour_iau,jdate_iau,
+ & colat1,SECSWR,SECLWR,slmsk, global_lats_r,lonsperlar)
endif ! sfcio _out
!
t4 = rtc()
@@ -279,11 +305,11 @@
if(me == IOPROC)
& call BAOPENWT(NO3D,d3df//CFHOUR,iostat)
if (hybrid .or. gen_coord_hybrid) then
- call WRT3D_hyb(IOPROC,no3d,ZHOUR,FHOUR,IDATE,colat1,
- & global_lats_r,lonsperlar,pl_coeff,
+ call WRT3D_hyb(IOPROC,no3d,zhour_iau,fhour_iau,jdate_iau,
+ & colat1,global_lats_r,lonsperlar,pl_coeff,
& SECSWR,SECLWR,sfc_fld%slmsk,flx_fld%psurf)
else
- call WRT3D(IOPROC,no3d,ZHOUR,FHOUR,IDATE,colat1,
+ call WRT3D(IOPROC,no3d,zhour_iau,fhour_iau,jdate_iau,colat1,
& global_lats_r,lonsperlar,pl_coeff,
& SECSWR,SECLWR,sl,si,sfc_fld%slmsk,flx_fld%psurf)
endif
Index: checkout/phys/gloopb.f
===================================================================
--- checkout/phys/gloopb.f (revision 93160)
+++ checkout/phys/gloopb.f (revision 94274)
@@ -8,7 +8,7 @@
& ozplin, jindx1, jindx2, ddy,
& phy_f3d, phy_f2d, phy_fctd, nctp,
& xlat, nblck, kdt, restart_step,
- & mdl_parm)
+ & mdl_parm, iniauinterval)
!!
!! Code Revision:
!! Sep 2009 Shrinivas Moorthi added nst_fld
@@ -70,6 +70,7 @@
! SPW drivers
! new idea_init
! new call of idea_phys, no Unified GW physics
+! Aug 17 2016 P Pegion - add logic for not doing mass tracer fix if in iau interval
!======================================================================
use resol_def
use layout1
@@ -124,7 +125,7 @@
!
integer nblck, kdt, nbdsw, nbdlw, nctp
- logical restart_step
+ logical restart_step,iniauinterval
!!
real(kind=kind_phys) tstep, phour,slag, sdec, cdec
real(kind=kind_phys) plyr(levs)
@@ -879,7 +880,7 @@
! For tracer fixer
do n=1,ntrac
- if (fixtrc(n)) then
+ if (fixtrc(n).and..not.iniauinterval) then
do k=1,levs
do i=1,njeff
trcp(lon+i-1,n,1) = trcp(lon+i-1,n,1)
@@ -1322,7 +1323,7 @@
enddo
else
do n=1,ntrac
- if (fixtrc(n)) then
+ if (fixtrc(n).and..not.iniauinterval) then
do k=1,levs
do i=1,njeff
trcp(lon+i-1,n,2) = trcp(lon+i-1,n,2)
@@ -1632,7 +1633,7 @@
tem = 0.5 / lonsperlar(lat)
!$omp parallel do private(j,n)
do n=1,ntrac
- if (fixtrc(n)) then
+ if (fixtrc(n).and..not.iniauinterval) then
trcj(lan,n,1) = 0.0
trcj(lan,n,2) = 0.0
do j=1,lons_lat
@@ -1652,7 +1653,7 @@
!$omp parallel do private(n,lat,tem)
do n=1,ntrac
- if (fixtrc(n)) then
+ if (fixtrc(n).and..not.iniauinterval) then
sumtrc(n,1) = 0.0
sumtrc(n,2) = 0.0
do lat=1,latr
Index: checkout/phys/checkiauforcing.f
===================================================================
--- checkout/phys/checkiauforcing.f (nonexistent)
+++ checkout/phys/checkiauforcing.f (revision 94274)
@@ -0,0 +1,44 @@
+ subroutine checkiauforcing(t,INIAUINTERVAL)
+
+! check to see if model is in IAU window
+ use namelist_physics_def,only : iaufhrs,iau_delthrs,iaufiles_anl
+ use machine, only:kind_evod
+ implicit none
+ real(kind_evod), intent(in) :: t
+ LOGICAL,intent(out) :: INIAUINTERVAL
+ real(kind_evod) dt
+ integer n,t1,t2,nfiles,nfilesall
+
+ dt = iau_delthrs*3600.
+ INIAUINTERVAL=.false.
+
+ nfilesall = size(iaufiles_anl)
+ nfiles = 0
+ do n=1,nfilesall
+ if (trim(iaufiles_anl(n)) .eq. '' .or. iaufhrs(n) .lt. 0) exit
+ nfiles = nfiles + 1
+ enddo
+ ! set forcing to zero and return if outside iau window.
+ if ( nfiles > 1) then ! IAU forcing files bookend interval
+ if (t <= iaufhrs(1)*3600. .or. t > iaufhrs(nfiles)*3600.) then
+ return
+ endif
+ else ! single file at middle of window
+ t1=iaufhrs(1)*3600 - dt*0.5
+ t2=iaufhrs(1)*3600 + dt*0.5
+ if ( t <= t1 .or. t > t2 ) then
+ return
+ endif
+ endif
+ if (nfiles > 1) then
+ if (t .eq. 3600.*iaufhrs(nfiles)) then
+ return
+ else if (t .eq. 3600.*iaufhrs(1)) then
+ return
+ endif
+ do n=1,nfiles
+ if (iaufhrs(n)*3600. > t) exit
+ enddo
+ endif
+ INIAUINTERVAL=.true.
+ end subroutine checkiauforcing
Index: checkout/dyn/namelist_dynamics_def.f
===================================================================
--- checkout/dyn/namelist_dynamics_def.f (revision 93160)
+++ checkout/dyn/namelist_dynamics_def.f (revision 94274)
@@ -34,7 +34,11 @@
logical ndslfv
! hmhj idea add
- logical lsidea
+ logical lsidea,iau
+ integer ::iaulnp=2 ! 0: no pressure update, 1: log pressure, 2: pressure
+ integer :: iau_delthrs = 6 ! iau time interval (to scale increments)
+ character(len=120), dimension(7) :: iaufiles_fg,iaufiles_anl
+ real(kind=kind_evod), dimension(7) :: iaufhrs
! WAM IPE coupling flags.
!------------------------
logical :: wam_ipe_coupling, height_dependent_g
Index: checkout/dyn/gfs_dynamics_grid_comp_mod.f
===================================================================
--- checkout/dyn/gfs_dynamics_grid_comp_mod.f (revision 93160)
+++ checkout/dyn/gfs_dynamics_grid_comp_mod.f (revision 94274)
@@ -28,6 +28,7 @@
! Dec 2013 Jun Wang, add restart in gfs_dynamics_start_time_get_mod argument list
! Feb 26 2016 S Moorthi add kdt_start and kdt_dif for grid-point digital filter
! logic. Filter works for initial state or for restart
+! Aug 2016 Add calls for iau forcing
!
! !interface:
!
@@ -48,6 +49,7 @@
use gfs_dyn_tracer_config, only: gfs_dyn_tracer
use namelist_dynamics_def, only: nemsio_in, ldfi_spect
USE gfs_dynamics_initialize_slg_mod, ONLY: gfs_dynamics_initialize_slg
+ use gfs_dyn_iau_module, only : getiauforcing,applyiauforcing,gq_iau
implicit none
@@ -516,10 +518,13 @@
logical,save :: first_reset=.true.
logical,save :: first_dfiend=.true.
TYPE(ESMF_TimeInterval) :: HALFDFIINTVAL
+ integer :: timestep_sec
+ real(kind_evod) dtiau
!! debug print for tracking import and export state (Sarah Lu)
TYPE(ESMF_Field) :: ESMFField !chlu_debug
TYPE(ESMF_FieldBundle) :: ESMFBundle !chlu_debug
+ real(kind_grid), dimension(:,:,:),allocatable :: grid_gr_iau
REAL(ESMF_KIND_R8), DIMENSION(:,:,:), POINTER :: fArr3D !chlu_debug
integer :: rc1, rcfinal, DFIHR &
, localPE,ii1,ii2,ii3 & !chlu_debug
@@ -584,6 +589,24 @@
! ' start_step=',int_state%start_step,' restart_step=',int_state%restart_step
!
IF(.NOT. int_state%restart_step .AND. .NOT. int_state%start_step ) THEN
+ IF (iau) THEN
+ allocate(grid_gr_iau(lonf,lats_node_a_max,gq_iau))
+ call esmf_clockget(clock, timestep = timestep, rc = rc1)
+ call esmf_timeintervalget(timestep, s = timestep_sec, rc = rc1)
+ IF (semilag) THEN
+ dtiau = timestep_sec
+ ELSE
+ dtiau = 2.*timestep_sec
+ ENDIF
+ int_state%iniauinterval=.false.
+ call getiauforcing(grid_gr_iau,int_state%phour*3600.,int_state,rc)
+ IF (rc == 0) THEN
+ IF (me == 0) print*,'applying IAU forcing'
+ int_state%iniauinterval=.true.
+ call applyiauforcing(grid_gr_iau,imp_gfs_dyn,int_state,dtiau,rc)
+ ENDIF
+ deallocate(grid_gr_iau)
+ ENDIF
IF(.NOT. int_state%reset_step) THEN
CALL gfs_dynamics_import2internal(imp_gfs_dyn, int_state, rc1)
ELSE
Index: checkout/dyn/gfs_dynamics_run_mod.f
===================================================================
--- checkout/dyn/gfs_dynamics_run_mod.f (revision 93160)
+++ checkout/dyn/gfs_dynamics_run_mod.f (revision 94274)
@@ -25,6 +25,7 @@
! dec 2014 Weiyu Yang, Add NEMSIO output.
! feb 2016 S Moorthi - grid-point digital filter fix to filter at initial
! time and during restart
+! aug 2016 Phil Pegion - add logic for iau
!
!
@@ -192,7 +193,8 @@
gis_dyn%dfiend_step, gis_dyn%pdryini, &
gis_dyn%nblck, gis_dyn%zhour, &
gis_dyn%pwat, gis_dyn%ptot, &
- gis_dyn%ptrc, gis_dyn%nfcstdate7)
+ gis_dyn%ptrc, gis_dyn%nfcstdate7, &
+ gis_dyn%iniauinterval)
gis_dyn%phour = fhour
@@ -293,7 +295,7 @@
gis_dyn%start_step ,gis_dyn% restart_step , &
gis_dyn%reset_step ,gis_dyn% end_step , &
gis_dyn% dfiend_step , & ! jw
- gis_dyn% nfcstdate7)
+ gis_dyn% nfcstdate7 ,gis_dyn%iniauinterval)
else
print *,' number of spectral loop is wrong. it is ', &
Index: checkout/dyn/grid_collect.f
===================================================================
--- checkout/dyn/grid_collect.f (revision 93160)
+++ checkout/dyn/grid_collect.f (revision 94274)
@@ -16,7 +16,7 @@
use gfs_dyn_write_state, only: buff_mult_pieceg,ngrid,ngridg
use gfs_dyn_layout1
use gfs_dyn_mpi_def
- use namelist_dynamics_def, only: gen_coord_hybrid
+ use namelist_dynamics_def, only: gen_coord_hybrid,iau
use gfs_dyn_coordinate_def, only: ak5, bk5
use gfs_dyn_physcons, rk => con_rocp
@@ -369,6 +369,11 @@
REAL(KIND=KIND_io4) ,allocatable :: buff_multg(:,:)
character * 16 :: recname1, reclevtyp1
integer reclev1
+! --- locals for iau:
+ integer :: jdate_iau(4),idat(8),jdat(8)
+ real (kind=kind_evod) :: fhour_iau
+ real :: rinc(5)
+
real(kind=kind_io4),allocatable :: tmp(:)
type(nemsio_gfile) :: gfileout
!
@@ -609,16 +614,36 @@
!end first
endif
!
+ print*,'in grid collect',iau,xhour
+ if ( iau .and. xhour >= 6.0 ) then
+ idat = 0
+ idat(1) = idate(4)
+ idat(2) = idate(2)
+ idat(3) = idate(3)
+ idat(5) = idate(1)
+ rinc = 0.
+ rinc(2) = 6.
+ call w3movdat(rinc,idat,jdat)
+ jdate_iau(4) = jdat(1)
+ jdate_iau(2) = jdat(2)
+ jdate_iau(3) = jdat(3)
+ jdate_iau(1) = jdat(5)
+ fhour_iau = xhour - 6.0
+ else
+ jdate_iau = idate
+ fhour_iau = xhour
+ endif
+
pdryini4 = pdryini
iorder_l = 2
irealf = 2
- yhour = xhour
+ yhour = fhour_iau
idvt = (ntoz-1) + 10 * (ntcw-1)
idate7=0
- idate7(1)=idate(4)
- idate7(2)=idate(2)
- idate7(3)=idate(3)
- idate7(4)=idate(1)
+ idate7(1)=jdate_iau(4)
+ idate7(2)=jdate_iau(2)
+ idate7(3)=jdate_iau(3)
+ idate7(4)=jdate_iau(1)
idate7(7)=100
!
nfhour=int(yhour)
Index: checkout/dyn/gfs_dyn_iau_mod.f
===================================================================
--- checkout/dyn/gfs_dyn_iau_mod.f (nonexistent)
+++ checkout/dyn/gfs_dyn_iau_mod.f (revision 94274)
@@ -0,0 +1,553 @@
+module gfs_dyn_iau_module
+
+! module for iau data, IO and time interpolation.
+ use ESMF
+ use gfs_dyn_layout1, only: me, len_trie_ls, len_trio_ls, ls_dim, nodes, &
+ lats_node_a, me_l_0,lats_node_a_max
+ use gfs_dyn_machine, only: kind_evod, kind_grid,kind_io4,kind_io8,kind_phys
+ use namelist_dynamics_def, only: iaufiles_fg,iaufiles_anl,iaufhrs,iau,iau_delthrs,ens_nam,nemsio_in,semilag
+ use gfs_dyn_resol_def
+ use gfs_dyn_mpi_def, only: mc_comp,mpi_sum,mpi_real4,mpi_complex,mpi_r_io_r,mpi_real8
+ use gfs_dynamics_states_mod, only: gfs_dynamics_internal_state
+ USE gfs_dynamics_namelist_mod, only: GFS_Dyn_State_Namelist
+ USE gfs_dynamics_grid_create_mod
+ use gfs_dyn_tracer_config, only: gfs_dyn_tracer
+ USE gfs_dynamics_add_get_state_mod
+ USE gfs_dynamics_err_msg_mod
+ USE gfs_dyn_coordinate_def, ONLY: ak5, bk5
+
+
+! iaufiles_fg: filenames for first-guess fields.
+! iaufiles_anl: filenames for analysis fields.
+! iaufhrs: forecast hours for input files.
+! iau_delthrs: length of IAU window (in hours).
+
+ implicit none
+ private
+
+ public :: init_iau, destroy_iau, getiauforcing,applyiauforcing
+
+ real(kind_evod), dimension(:,:,:,:),allocatable :: grid_gr_iauall
+ integer, public :: nfiles,gq_iau,gu_iau,gv_iau,gt_iau,grq_iau
+ logical, public :: iau_initialized = .false.,iauinc=.false.
+
+
+
+ contains
+
+ subroutine init_iau(int_state)
+
+! read in first-guess and analysis files, compute and store increments in
+! spectral space.
+
+ type(gfs_dynamics_internal_state),intent(in) :: int_state
+ integer, allocatable, dimension(:) :: idt
+ integer n,nfilesall,idate(4),ierr
+
+ integer lan,k,lat,i,lon_dim,lons_lat
+
+ gu_iau=1 ! u wind
+ gv_iau=gu_iau+levs ! v wind
+ gt_iau=gv_iau+levs ! virtual temperature
+ grq_iau=gt_iau+levs ! tracers
+ gq_iau=grq_iau+levh ! surface pressure
+ if (me.EQ.me_l_0) print*,'gu_iau,gv_iau,gt_iau=',gu_iau,gv_iau,gt_iau
+ if (me.EQ.me_l_0) print*,'grq_iau,gq_iau=',grq_iau,gq_iau
+
+ iau_initialized = .true.
+ iauinc=.false.
+ if (trim(iaufiles_fg(1)) .eq. '' .and. trim(iaufiles_anl(1)) .ne. '') then
+ iauinc=.true.
+ endif
+
+ nfilesall = size(iaufiles_anl)
+ nfiles = 0
+ do n=1,nfilesall
+ if (trim(iaufiles_anl(n)) .eq. '' .or. iaufhrs(n) .lt. 0) exit
+ if (me .eq. me_l_0) then
+ print *,n,trim(adjustl(iaufiles_anl(n)))
+ if (.not. iauinc) print *,n,trim(adjustl(iaufiles_fg(n)))
+ endif
+ nfiles = nfiles + 1
+ enddo
+ if (me .eq. me_l_0) print *,'nfiles = ',nfiles
+ call mpi_barrier(mc_comp,ierr)
+ if (nfiles < 1) then
+ print *,'must be at least one file in iaufiles_fg and iaufiles_anal'
+ call mpi_quit(9999)
+ endif
+ allocate(idt(nfiles-1))
+ idt = iaufhrs(2:nfiles)-iaufhrs(1:nfiles-1)
+ do n=1,nfiles-1
+ if (idt(n) .ne. iaufhrs(2)-iaufhrs(1)) then
+ print *,'forecast intervals in iaufhrs must be constant'
+ call mpi_quit(9999)
+ endif
+ enddo
+ if (me .eq. me_l_0) print *,'iau interval = ',iau_delthrs,' hours'
+ deallocate(idt)
+ allocate(grid_gr_iauall(lonf,lats_node_a_max,gq_iau,nfiles))
+ if (nemsio_in) then
+ call read_iau_nemsio(int_state)
+ else
+ call read_iau_sigio(int_state)
+ endif
+ end subroutine init_iau
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine read_iau_nemsio(int_state)
+! read in first-guess and analysis files, compute and store increments in
+! grid space.
+ implicit none
+ type(gfs_dynamics_internal_state),intent(in) :: int_state
+ real(kind_grid), dimension(:,:,:),allocatable :: grid_gr_fg
+ character(len=120) filename
+ integer n,iprint,idate(4),ierr
+
+ REAL(KIND=KIND_GRID) ptrc (lonf,lats_node_a,ntrac) !glbsum
+
+ iprint = 1
+ if (me.eq.me_l_0) print*,'reading iau_nemsio',iauinc,me
+ if (.not. iauinc) allocate(grid_gr_fg(lonf,lats_node_a_max,gq_iau))
+ if (me.eq.me_l_0) print*,'allocated',lonf,lats_node_a_max,gq_iau
+ do n=1,nfiles
+ if (.not. iauinc) then ! skip reading first guess and only read in increment file
+ filename = iaufiles_fg(n)
+ if (int_state%ens) then
+ filename = trim(filename) // ens_nam
+ endif
+ if (me .eq. me_l_0) print *,'reading fg ',trim(filename)
+ CALL TREADEO_nemsio_iau(filename,IDATE, &
+ grid_gr_fg(:,1:lats_node_a,gq_iau), &
+ grid_gr_fg(:,1:lats_node_a,gu_iau:gu_iau+levs-1), &
+ grid_gr_fg(:,1:lats_node_a,gv_iau:gv_iau+levs-1), &
+ grid_gr_fg(:,1:lats_node_a,gt_iau:gt_iau+levs-1), &
+ grid_gr_fg(:,1:lats_node_a,grq_iau:grq_iau+levh-1), &
+ int_state%LS_NODE,int_state%LS_NODES,int_state%MAX_LS_NODES,IPRINT, &
+ int_state%global_lats_a,int_state%lats_nodes_a,int_state%lonsperlat)
+ if (me .eq. me_l_0)print*,'fg grid_gr ps',me, grid_gr_fg(1,1,gq_iau),grid_gr_fg(2,lats_node_a,gq_iau)
+ if(iaulnp.EQ.1) grid_gr_fg(:,:,gq_iau)=alog(grid_gr_fg(:,:,gq_iau))
+ endif ! not iauinc
+ filename = iaufiles_anl(n) ! analysis or increment file
+ if (int_state%ens) then
+ filename = trim(filename) // ens_nam
+ endif
+ if (me .eq. me_l_0) print *,'reading anl ',trim(filename)
+ CALL TREADEO_nemsio_iau(filename,IDATE, &
+ grid_gr_iauall(:,1:lats_node_a,gq_iau,n), &
+ grid_gr_iauall(:,1:lats_node_a,gu_iau:gu_iau+levs-1,n), &
+ grid_gr_iauall(:,1:lats_node_a,gv_iau:gv_iau+levs-1,n), &
+ grid_gr_iauall(:,1:lats_node_a,gt_iau:gt_iau+levs-1,n), &
+ grid_gr_iauall(:,1:lats_node_a,grq_iau:grq_iau+levh-1,n), &
+ int_state%LS_NODE,int_state%LS_NODES,int_state%MAX_LS_NODES,IPRINT, &
+ int_state%global_lats_a,int_state%lats_nodes_a,int_state%lonsperlat)
+ if(iaulnp.EQ.1) grid_gr_iauall(:,:,gq_iau,n)=alog(grid_gr_iauall(:,:,gq_iau,n))
+ if (me.EQ.me_l_0) print*,'anl grid_gr ps',me, grid_gr_iauall(1,1,gq_iau,n),grid_gr_iauall(2,lats_node_a,gq_iau,n)
+!!$omp workshare
+ if (.not. iauinc)grid_gr_iauall(:,:,:,n)=grid_gr_iauall(:,:,:,n)-grid_gr_fg(:,:,:)
+!!$omp end workshare
+ enddo
+ if (me.EQ.me_l_0) print*,'iau grid_gr ps',me, grid_gr_iauall(1,1,gq_iau,1),grid_gr_iauall(2,lats_node_a,gq_iau,1)
+ call mpi_barrier(mc_comp,ierr)
+ if (.not.iauinc) deallocate(grid_gr_fg)
+ end subroutine read_iau_nemsio
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine read_iau_sigio(int_state)
+
+! read in first-guess and analysis files, compute and store increments in
+! spectral space.
+
+ implicit none
+ type(gfs_dynamics_internal_state),intent(in) :: int_state
+ real(kind_grid), dimension(:,:,:),allocatable :: grid_gr_fg
+ real(kind_grid), dimension(:,:,:),allocatable :: grid_gr_anl
+ real(kind_grid), dimension(:,:) ,allocatable :: grid_gr_ps
+ real(kind_grid), dimension(:,:,:),allocatable :: grid_gr_tr
+ real(kind_grid), dimension(:,:,:),allocatable :: grid_gr_dummy
+ character(len=120) filename
+ integer n,iprint,idate(4),ierr
+
+ integer N1,lan,k,lat,i,lon_dim,lons_lat,j
+
+ iprint=1
+ if (.not. iauinc) allocate(grid_gr_fg(lonf,lats_node_a_max,gq_iau))
+ allocate(grid_gr_anl(lonf,lats_node_a_max,gq_iau))
+ allocate(grid_gr_tr(lonf,lats_node_a_max,levh))
+ allocate(grid_gr_ps(lonf,lats_node_a_max))
+ allocate(grid_gr_dummy(lonf,lats_node_a_max,levs))
+ do n=1,nfiles
+ if (.not. iauinc) then ! skip reading first guess and only read in increment file
+ filename = iaufiles_fg(n)
+ if (int_state%ens) then
+ filename = trim(filename) // ens_nam
+ endif
+ if (me .eq. me_l_0) print *,'reading ',trim(filename)
+ CALL TREADEO_IAU(IDATE,grid_gr_fg, &
+ int_state%LS_NODE,int_state%LS_NODES,int_state%MAX_LS_NODES,IPRINT, &
+ int_state%global_lats_a,int_state%lats_nodes_a,int_state%lonsperlat, filename, &
+ int_state%epse, int_state%epso, int_state%epsedn, &
+ int_state%epsodn, int_state%plnew_a, int_state%plnow_a, &
+ int_state%plnev_a, int_state%plnod_a,int_state%snnp1ev, int_state%snnp1od)
+ if (me .eq. me_l_0) print *,'done reading ',trim(filename)
+ if (me .eq. me_l_0) print *,'calling model_to_common_vars'
+ call model_to_common_vars (grid_gr_fg(1,1,gq_iau), &
+ grid_gr_fg(1,1,gt_iau), &
+ grid_gr_fg(1,1,grq_iau), &
+ grid_gr_fg(1,1,gu_iau), &
+ grid_gr_fg(1,1,gv_iau), &
+ grid_gr_dummy(1,1,1), &
+ grid_gr_dummy(1,1,1), &
+ grid_gr_dummy(1,1,1), &
+ int_state%global_lats_a, &
+ int_state%lonsperlat,1)
+ if(iaulnp.EQ.1) grid_gr_fg(:,:,gq_iau)=alog(grid_gr_fg(:,:,gq_iau))
+ if (me .eq. me_l_0) print *,'back from model_to_common_vars'
+ call mpi_barrier(mc_comp,i)
+ endif
+ filename = iaufiles_anl(n) ! analysis or increment file
+ if (int_state%ens) then
+ filename = trim(filename) // ens_nam
+ endif
+ if (me.eq.me_l_0) print*,'calling treadeo',trim(filename)
+ CALL TREADEO_IAU(IDATE,grid_gr_anl, &
+ int_state%LS_NODE,int_state%LS_NODES,int_state%MAX_LS_NODES,IPRINT, &
+ int_state%global_lats_a,int_state%lats_nodes_a,int_state%lonsperlat, filename, &
+ int_state%epse, int_state%epso, int_state%epsedn, &
+ int_state%epsodn, int_state%plnew_a, int_state%plnow_a, &
+ int_state%plnev_a, int_state%plnod_a,int_state%snnp1ev, int_state%snnp1od)
+ if (me.eq.me_l_0) print*,'back from treadeo',trim(filename)
+ if (iauinc) then
+ grid_gr_tr=int_state%grid_gr(:,:,g_rq:g_rq+3*levs-1)
+ else
+ grid_gr_tr=grid_gr_anl(:,:,grq_iau:grq_iau+3*levs-1)
+ endif
+ grid_gr_ps=grid_gr_anl(:,:,gq_iau)
+ call model_to_common_vars (grid_gr_ps(1,1), &
+ grid_gr_anl(1,1,gt_iau), &
+ grid_gr_tr(1,1,1), &
+ grid_gr_anl(1,1,gu_iau), &
+ grid_gr_anl(1,1,gv_iau), &
+ grid_gr_dummy(1,1,1), &
+ grid_gr_dummy(1,1,1), &
+ grid_gr_dummy(1,1,1), &
+ int_state%global_lats_a, &
+ int_state%lonsperlat,1)
+ ! copy back converted pressure if not increments
+ if (.not. iauinc) grid_gr_anl(:,:,gq_iau)=grid_gr_ps
+ if(iaulnp.EQ.1) grid_gr_anl(:,:,gq_iau)=alog(grid_gr_anl(:,:,gq_iau))
+ print*,'b4 grid_gr ps',me, grid_gr_anl(1,1,gq_iau),grid_gr_fg(1,1,gq_iau)
+!!$omp workshare
+ if (iauinc) then
+ grid_gr_iauall(:,:,gq_iau,n)=grid_gr_anl(:,:,gq_iau)
+ grid_gr_iauall(:,:,gu_iau:gu_iau+levs-1,n)=grid_gr_anl(:,:,gu_iau:gu_iau+levs-1)
+ grid_gr_iauall(:,:,gv_iau:gv_iau+levs-1,n)=grid_gr_anl(:,:,gv_iau:gv_iau+levs-1)
+ grid_gr_iauall(:,:,gt_iau:gt_iau+levs-1,n)=grid_gr_anl(:,:,gt_iau:gt_iau+levs-1)
+ grid_gr_iauall(:,:,grq_iau:grq_iau+levh-1,n)=grid_gr_anl(:,:,grq_iau:grq_iau+levh-1)
+ else
+ grid_gr_iauall(:,:,gq_iau,n)=grid_gr_anl(:,:,gq_iau)-grid_gr_fg(:,:,gq_iau)
+ grid_gr_iauall(:,:,gu_iau:gu_iau+levs-1,n)=grid_gr_anl(:,:,gu_iau:gu_iau+levs-1)-grid_gr_fg(:,:,gu_iau:gu_iau+levs-1)
+ grid_gr_iauall(:,:,gv_iau:gv_iau+levs-1,n)=grid_gr_anl(:,:,gv_iau:gv_iau+levs-1)-grid_gr_fg(:,:,gv_iau:gv_iau+levs-1)
+ grid_gr_iauall(:,:,gt_iau:gt_iau+levs-1,n)=grid_gr_anl(:,:,gt_iau:gt_iau+levs-1)-grid_gr_fg(:,:,gt_iau:gt_iau+levs-1)
+ grid_gr_iauall(:,:,grq_iau:grq_iau+levh-1,n)=grid_gr_anl(:,:,grq_iau:grq_iau+levh-1)-grid_gr_fg(:,:,grq_iau:grq_iau+levh-1)
+ endif
+!!$omp end workshare
+ enddo
+ if (me.EQ.me_l_0) print*,'iau grid_gr ps',me, grid_gr_iauall(1,1,gq_iau,1),grid_gr_iauall(2,lats_node_a,gq_iau,1)
+ if (me.EQ.me_l_0) print*,'iau grid_gr t',me, grid_gr_iauall(1,1,gt_iau+30,1),grid_gr_iauall(2,lats_node_a,gt_iau+30,1)
+ if (me.EQ.me_l_0) print*,'iau grid_gr u',me, grid_gr_iauall(1,1,gu_iau+30,1),grid_gr_iauall(2,lats_node_a,gu_iau+30,1)
+ if (me.EQ.me_l_0) print*,'iau grid_gr q',me, grid_gr_iauall(1,1,grq_iau+30,1),grid_gr_iauall(2,lats_node_a,grq_iau+30,1)
+ if (.not.iauinc) deallocate(grid_gr_fg)
+ deallocate(grid_gr_anl)
+ deallocate(grid_gr_ps,grid_gr_tr)
+ end subroutine read_iau_sigio
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ subroutine getiauforcing(grid_gr_iau,t,int_state,rc)
+
+! compute an IAU forcing by interpolating increments to model time set
+! and dividing by length of IAU window (in seconds).
+
+ implicit none
+ REAL(KIND=KIND_GRID),intent(out) :: grid_gr_iau(lonf,lats_node_a_max,gq_iau)
+ type(gfs_dynamics_internal_state),intent(in) :: int_state
+ real(kind_evod), intent(in) :: t
+ INTEGER,intent(out) :: rc
+ real(kind_evod) delt, dt
+ integer n,t1,t2
+ integer i,j,k,lat
+ !real(kind_io4), allocatable, dimension(:,:) :: workg,workg_out
+ !real(kind_phys), allocatable, dimension(:,:) :: grid_gr_tmp
+ !real (kind=kind_io8) glolal(lonf,lats_node_a)
+ !integer kmsk0(lonf,lats_node_a)
+ !kmsk0 = 0
+
+ grid_gr_iau = 0.
+ dt = iau_delthrs*3600.
+ if (me .eq. me_l_0) print *,'in getiauforcing1',nfiles,iaufhrs(1:nfiles)
+ if (me .eq. me_l_0) print *,'in getiauforcing2',t, iaufhrs(1)*3600,iaufhrs(nfiles)*3600.
+ ! set forcing to zero and return if outside iau window.
+ if ( nfiles > 1) then ! IAU forcing files bookend interval
+ if (t <= iaufhrs(1)*3600. .or. t > iaufhrs(nfiles)*3600.) then
+ if (me .eq. me_l_0) print *,'no iau forcing'
+ rc=1
+ return
+ endif
+ else ! single file at middle of window
+ t1=iaufhrs(1)*3600 - dt*0.5
+ t2=iaufhrs(1)*3600 + dt*0.5
+ if ( t <= t1 .or. t > t2 ) then
+ if (me .eq. me_l_0) print *,'no iau forcing'
+ rc=1
+ return
+ endif
+ endif
+ if (nfiles > 1) then
+ if (t .eq. 3600.*iaufhrs(nfiles)) then
+!!$omp workshare
+ grid_gr_iau = grid_gr_iauall(:,:,:,nfiles)/dt
+!!$omp end workshare
+ return
+ else if (t .eq. 3600.*iaufhrs(1)) then
+!!$omp workshare
+ grid_gr_iau = grid_gr_iauall(:,:,:,1)/dt
+!!$omp end workshare
+ return
+ endif
+ do n=1,nfiles
+ if (iaufhrs(n)*3600. > t) exit
+ enddo
+ if (me .eq. me_l_0) print *,'n,t,to',n,t/3600.,iaufhrs(n)
+ delt = (iaufhrs(n)-(t/3600.))/(iaufhrs(n)-iaufhrs(n-1))
+!!$omp workshare
+ grid_gr_iau = ((1.-delt)*grid_gr_iauall(:,:,:,n) + delt*grid_gr_iauall(:,:,:,n-1))/dt
+!!$omp end workshare
+ if (me .eq. me_l_0) print *,'getiauforcing:',t/3600.,1.-delt,n,iaufhrs(n),delt,n-1,iaufhrs(n-1)
+ else
+ grid_gr_iau = grid_gr_iauall(:,:,:,1)/dt
+ endif
+ if (me .eq. me_l_0) print *,'have iau forcing',grid_gr_iau(1,1,gq_iau)
+ rc=0
+ !allocate(workg(lonf,latg))
+ !allocate(workg_out(lonf,latg))
+ !allocate(grid_gr_tmp(lonf,lats_node_a_max))
+! write out data
+ !if (me.eq.me_l_0) open(77,file='../iau_forcing.dat',form='unformatted')
+ !do k=1,gq_iau
+ ! workg= 0.
+ ! grid_gr_tmp(:,:)=grid_gr_iauall(:,:,k,1)
+ ! CALL uninterpred(2,kmsk0,glolal,grid_gr_tmp,&
+ ! int_state%global_lats_a,int_state%lonsperlat)
+ ! do j=1,lats_node_a
+ ! lat=int_state%global_lats_a(ipt_lats_node_a-1+j)
+ ! do i=1,lonf
+ ! workg(i,lat) = glolal(i,j)
+ ! enddo
+ ! enddo
+ ! call mpi_reduce(workg,workg_out,lonf*latg,&
+ ! mpi_real4,mpi_sum,me_l_0,mc_comp,i)
+ ! if (me .eq. me_l_0) then
+ ! write(77) workg_out
+ ! print*,'array a k',grid_gr_iauall(1,1,k,1),glolal(1,1)
+ ! print*,'writing k',workg(1,1),workg_out(1,1)
+ ! endif
+ !enddo
+ !close(77)
+ !deallocate(workg,workg_out)
+ !call mpi_barrier(mc_comp,i)
+ !call mpi_quit(9999)
+! end write
+ end subroutine getiauforcing
+
+SUBROUTINE applyiauforcing(grid_gr_iau,imp_gfs_dyn,int_state,deltim,rc)
+
+! this subroutine added the iau forcing term to the gfs import state
+!-----------------------------------------------------------
+
+ implicit none
+ type(esmf_state), intent(inout) :: imp_gfs_dyn
+ type(gfs_dynamics_internal_state),intent(in) :: int_state
+ REAL(KIND=KIND_GRID),intent(in) :: grid_gr_iau(lonf,lats_node_a_max,gq_iau)
+ REAL(kind_evod),intent(in) :: deltim
+ INTEGER, OPTIONAL, INTENT(out) :: rc
+
+ TYPE(ESMF_Field) :: Field
+ TYPE(ESMF_FieldBundle) :: Bundle
+ TYPE(GFS_Dyn_State_Namelist) :: cf
+
+ INTEGER :: rc1, rcfinal,irec
+ INTEGER :: k, kstr, kend,i,j
+ INTEGER :: mstr, mend,lat,ierr
+
+ REAL, DIMENSION(:, :), POINTER :: FArr2D
+ REAL, DIMENSION(:, :, :), POINTER :: FArr3D
+ ! real(kind_io4), allocatable, dimension(:,:) :: workg,workg_out
+ ! real(kind_phys), allocatable, dimension(:,:) :: grid_gr_tmp
+ ! real (kind=kind_io8) glolal(lonf,lats_node_a)
+ ! integer kmsk0(lonf,lats_node_a)
+ ! character*3,fstr
+ ! kmsk0 = 0
+
+! initialize the error signal variables.
+!---------------------------------------
+ rc1 = esmf_success
+ rcfinal = esmf_success
+
+ CALL esmf_logwrite( &
+ " update import state with IAU forcing", &
+ ESMF_LOGMSG_INFO, rc = rc1)
+
+ cf = int_state%esmf_sta_list
+
+! get the surface pressure array from the esmf import state.
+! for the detailed comments for every computational steps
+! please refer to the surface orography array code.
+!-----------------------------------------------------------
+ !write(fstr,FMT='(I3.3)') int_state%kdt
+ !if (me.eq.me_l_0) open(78,file='../psiau_'//fstr//'.dat',form='unformatted')
+ !allocate(workg(lonf,latg))
+ !allocate(workg_out(lonf,latg))
+ !allocate(grid_gr_tmp(lonf,lats_node_a_max))
+! write out data
+ !workg= 0.
+ !grid_gr_tmp(:,:)=grid_gr_iau(:,:,gq_iau)
+ !CALL uninterpred(2,kmsk0,glolal,grid_gr_tmp,&
+ ! int_state%global_lats_a,int_state%lonsperlat)
+ ! do j=1,lats_node_a
+ ! lat=int_state%global_lats_a(ipt_lats_node_a-1+j)
+ ! do i=1,lonf
+ ! workg(i,lat) = glolal(i,j)
+ ! enddo
+ ! enddo
+ ! call mpi_reduce(workg,workg_out,lonf*latg,&
+ ! mpi_real4,mpi_sum,me_l_0,mc_comp,i)
+ ! if (me .eq. me_l_0) then
+ ! write(78) workg_out
+ ! endif
+ IF(cf%ps_import == 1) THEN
+ CALL getf90arrayfromstate(imp_gfs_dyn, 'ps', FArr2D, 0, rc = rc1)
+ CALL gfs_dynamics_err_msg(rc1,"get esmf state - ps",rcfinal)
+! write out data
+ ! workg= 0.
+ ! grid_gr_tmp(:,:)=FArr2D
+ ! CALL uninterpred(2,kmsk0,glolal,grid_gr_tmp,&
+ ! int_state%global_lats_a,int_state%lonsperlat)
+ ! do j=1,lats_node_a
+ ! lat=int_state%global_lats_a(ipt_lats_node_a-1+j)
+ ! do i=1,lonf
+ ! workg(i,lat) = glolal(i,j)
+ ! enddo
+ ! enddo
+ !call mpi_reduce(workg,workg_out,lonf*latg,&
+ ! mpi_real4,mpi_sum,me_l_0,mc_comp,i)
+ !if (me .eq. me_l_0) then
+ ! write(78) workg_out
+ !endif
+!! add iau forcing
+ if (me .eq. me_l_0) print *,'before forcing ps:',farr2d(4,4)
+ If (iaulnp.GT.0) THEN
+ If (iaulnp.EQ.1) THEN
+ FArr2D = exp(alog(FArr2D) + grid_gr_iau(:,:,gq_iau) * deltim)
+ ELSE
+ IF (iauinc .and. .not.semilag) THEN
+ FArr2D = FArr2D * (1.0 + grid_gr_iau(:,:,gq_iau) * deltim)
+!
+ ELSE
+ FArr2D = FArr2D + grid_gr_iau(:,:,gq_iau) * deltim
+ ENDIF
+ ENDIF
+ ENDIF
+ if (me .eq. me_l_0) print *,'applyiauforcing ps:',grid_gr_iau(4,4,gq_iau),farr2d(4,4)
+ END IF
+! write out data
+ !workg= 0.
+ !grid_gr_tmp(:,:)=FArr2D
+ !CALL uninterpred(2,kmsk0,glolal,grid_gr_tmp,&
+ ! int_state%global_lats_a,int_state%lonsperlat)
+ !do j=1,lats_node_a
+ ! lat=int_state%global_lats_a(ipt_lats_node_a-1+j)
+ ! do i=1,lonf
+ ! workg(i,lat) = glolal(i,j)
+ ! enddo
+ !enddo
+ !call mpi_reduce(workg,workg_out,lonf*latg,&
+ ! mpi_real4,mpi_sum,me_l_0,mc_comp,i)
+ !if (me .eq. me_l_0) then
+ ! write(78) workg_out
+ ! close(78)
+ !endif
+ !deallocate(workg,workg_out)
+
+! get the temperature array from the esmf import state.
+!------------------------------------------------------
+ IF(cf%temp_import == 1) THEN
+ CALL getf90arrayfromstate(imp_gfs_dyn, 't', FArr3D, 0, rc = rc1)
+ if (me .eq. me_l_0) print *,'before forcing t:',farr3d(4,4,4)
+ CALL gfs_dynamics_err_msg(rc1,"get esmf state - t",rcfinal)
+ FArr3D = FArr3D + grid_gr_iau(:,:,gt_iau:gt_iau+levs-1) * deltim
+ if (me .eq. me_l_0) print *,'applyiauforcing t:',grid_gr_iau(4,4,gt_iau+3),farr3d(4,4,4)
+ END IF
+
+! get the zonal-wind array from the esmf import state.
+! for detailed line by line comments please refer to
+! the temperature code.
+!-----------------------------------------------------
+ IF(cf%u_import == 1) THEN
+ CALL getf90arrayfromstate(imp_gfs_dyn, 'u', FArr3D, 0, rc = rc1)
+ if (me .eq. me_l_0) print *,'before forcing u:',farr3d(4,4,4)
+ CALL gfs_dynamics_err_msg(rc1,"get esmf state - u",rcfinal)
+ FArr3D = FArr3D + grid_gr_iau(:,:,gu_iau:gu_iau+levs-1) * deltim
+ if (me .eq. me_l_0) print *,'applyiauforcing u:',grid_gr_iau(4,4,gu_iau+3),farr3d(4,4,4)
+ END IF
+
+! get the meridian-wind array from the esmf import state.
+!-----------------------------------------------------
+ IF(cf%v_import == 1) THEN
+ CALL getf90arrayfromstate(imp_gfs_dyn, 'v', FArr3D, 0, rc = rc1)
+ if (me .eq. me_l_0) print *,'before forcing v:',farr3d(4,4,4)
+ CALL gfs_dynamics_err_msg(rc1,"get esmf state - v",rcfinal)
+ FArr3D = FArr3D + grid_gr_iau(:,:,gv_iau:gv_iau+levs-1) * deltim
+ if (me .eq. me_l_0) print *,'applyiauforcing v:',grid_gr_iau(4,4,gv_iau+3),farr3d(4,4,4)
+ END IF
+
+! get the tracer array from the esmf import state.
+!-------------------------------------------------
+ IF(cf%tracer_import == 1) THEN
+ CALL ESMF_StateGet(imp_gfs_dyn, 'tracers', Bundle, rc = rc1 )
+ CALL gfs_dynamics_err_msg(rc1, 'retrieve Ebundle from state', rcfinal)
+ DO k = 1, int_state%ntrac
+ IF(ASSOCIATED(FArr3D)) NULLIFY(FArr3D)
+ CALL ESMF_FieldBundleGet(Bundle, &
+ trim(gfs_dyn_tracer%vname(k, 1)), &
+ field = Field, &
+ rc = rc1)
+ CALL gfs_dynamics_err_msg(rc1, 'retrieve Efield from bundle', rcfinal)
+
+ CALL ESMF_FieldGet(Field, farrayPtr = FArr3D, localDE = 0, rc = rc1)
+
+ CALL gfs_dynamics_err_msg(rc1, 'retrieve Farray from field', rcfinal)
+ kstr = grq_iau + (k-1)*levs
+ kend = kstr + levs - 1
+ if (me .eq. me_l_0) print *,'before forcing tr:',k,farr3d(4,4,4)
+ FArr3D = FArr3D + grid_gr_iau(:,:,kstr:kend) * deltim
+ if (me .eq. me_l_0) print *,'applyiauforcing tr:',grid_gr_iau(4,4,kstr+3),farr3d(4,4,4)
+ END DO
+ END IF
+
+!
+!
+! PRINT out the final error signal message and put it to rc.
+!-----------------------------------------------------------
+ IF(PRESENT(rc)) CALL gfs_dynamics_err_msg_final(rcfinal, &
+ "applyiauforcing",rc)
+
+ END SUBROUTINE applyiauforcing
+
+
+ subroutine destroy_iau()
+
+ implicit none
+! deallocate array.
+ deallocate(grid_gr_iauall)
+
+ end subroutine destroy_iau
+
+end module gfs_dyn_iau_module
Index: checkout/dyn/gfs_dynamics_initialize_slg_mod.f
===================================================================
--- checkout/dyn/gfs_dynamics_initialize_slg_mod.f (revision 93160)
+++ checkout/dyn/gfs_dynamics_initialize_slg_mod.f (revision 94274)
@@ -45,6 +45,7 @@
use gfs_dynamics_getcf_mod
use gfs_dyn_machine, only : kind_io4
use nemsio_module , only : nemsio_init
+ use gfs_dyn_iau_module, only: init_iau
!
use gfs_dyn_write_state, only : buff_mult_pieceg
use gfs_dyn_layout1, only : ipt_lats_node_a, lats_node_a_max,lon_dim_a
@@ -867,6 +868,9 @@
!!
!!
gis_dyn%zhour = fhour
+ if (iau) then
+ call init_iau(gis_dyn)
+ endif
rc = 0
dyn_ini_time = dyn_ini_time + (timef() - btime)
Index: checkout/dyn/treadeo_nemsio_iau.f
===================================================================
--- checkout/dyn/treadeo_nemsio_iau.f (nonexistent)
+++ checkout/dyn/treadeo_nemsio_iau.f (revision 94274)
@@ -0,0 +1,288 @@
+ subroutine treadeo_nemsio_iau(cfile,IDATE
+ &, psg,uug,vvg,ttg,rqg
+ &, LS_NODE,LS_NODES,MAX_LS_NODES
+ &, IPRINT
+ &, global_lats_a,lats_nodes_a
+ &, lonsperlat)
+
+!!
+!! Revision history:
+! 2008 Henry Juang, original code
+! Nov 23 2009 Sarah Lu, tracer read-in is generalized (loop through ntrac, with
+! tracer name specified in gfs_dyn_tracer_config)
+! Apr 09 2010 Sarah Lu, set rqg initial value to 1.e-15
+! Aug 17 2010 Sarah Lu, clean debug print
+! Aug 25 2010 Sarah Lu, modified to compute tracer global sum
+! Sep 08 2010 Jun Wang, change to nemsio format file
+! Dec 16 2010 Jun Wang, change to nemsio library
+! Feb 20 2011 Henry Juang, change to have mass dp and ndslfv
+! Feb 26 2011 Sarah Lu, modify to read both cold-start (from chgres) and
+! warm-start (from replay) ICs
+! Jun 2011 Jun Wang reading fields from grib file with either w3_d
+! or w3_4 lib
+! Nov 11 2011 Sarah Lu, change floor value for tracer initial values;
+! remove nvcoord read-in and check
+! Sep 20 2012 Jun Wang, set n time step trie/o to be consistnet with sigio input
+! Feb 18 2015 S Moorthi - adapted to two time level semi_Lagrangian model
+!
+
+ use gfs_dyn_resol_def
+ use gfs_dyn_layout1
+ use gfs_dyn_coordinate_def
+ use namelist_dynamics_def
+ use gfs_dyn_vert_def
+ use gfs_dyn_mpi_def
+ use gfs_dyn_physcons, rkap => con_rocp
+ &, cpd => con_cp
+ use nemsio_module
+ use nemsio_def
+ use gfs_dyn_tracer_config, ONLY : gfs_dyn_tracer ! generalized tracer
+ use gfs_dyn_io_header
+
+!
+ IMPLICIT NONE
+ character*(*) cfile
+ INTEGER IDATE(4), idate7(7)
+ &, levsi, jcapi, latgi, lonfi
+ logical fromiau
+ integer ls_node(ls_dim,3)
+!
+ INTEGER, dimension(nodes) :: MAX_LS_NODES, lats_nodes_a
+ INTEGER LS_NODES(LS_DIM,NODES)
+ INTEGER J,K,L,LOCL,N,lv,kk,w3rlkind,w3ikind,iprint
+ &, i,lan,lat,iblk,lons_lat,il,lon,njeff,nn
+ &, indev,indod,indev1,indev2,indod1,indod2
+ &, nfhour,nfminute,nfsecondn,nfsecondd
+!
+
+! for generalized tracers
+ integer nreci
+ character*8 :: vname
+ character*8, allocatable :: recnamei(:), reclevtypi(:)
+ integer, allocatable :: reclevi(:)
+!
+ integer iret, num_dta, ijm, tlmeta
+ real(kind=kind_evod) ga2, psurfff, pressk, tem
+ real(kind=kind_evod), parameter :: rkapi=1.0/rkap,
+ & rkapp1=1.0+rkap
+!
+ integer kmsk(lonf,latg), global_lats_a(latg), lonsperlat(latg)
+ real(kind=kind_io8), dimension(lonf,lats_node_a) :: buffo, buff2
+!
+ real (kind=kind_io4), allocatable :: nemsio_data(:)
+!!
+ real(kind=kind_grid), dimension(lonf,lats_node_a) :: psg
+ real(kind=kind_grid), dimension(lonf,lats_node_a,levs) :: uug
+ &, vvg, ttg
+ real(kind=kind_grid) rqg(lonf,lats_node_a,levh)
+!
+! REAL(KIND=KIND_GRID) ptrc (lonf,lats_node_a,ntrac) !glbsum
+!
+ real(8) timef,stime,etime
+!
+ INTEGER INDLSEV,JBASEV,INDLSOD,JBASOD
+ INCLUDE 'function2'
+!------------------------------------------------------------------------
+! Input file is in grid-point space - use gfs_io package
+!
+ if (me == 0) write(0,*)' before nemsio_open cfile=',cfile
+
+ stime = timef()
+ call nemsio_open(gfile_in,trim(cfile),'read',iret)
+ etime = timef()
+
+ if (me == 0) write(0,*)'in read nemsio file, open time='
+ &, timef()-stime
+!
+ call nemsio_getfilehead(gfile_in,iret=iret,
+ & version=ivsupa,idate=idate7,
+ & nfhour=nfhour,nfminute=nfminute,
+ & nfsecondn=nfsecondn,nfsecondd=nfsecondd,
+ & dimy=latgi,dimx=lonfi,dimz=levsi,
+ & jcap=jcapi,idvc=idvc,
+ & ncldt=ncldt,tlmeta=tlmeta)
+
+! write(0,*)' after nemsio_getfilehead levsi=',levsi,latgi,lonfi
+! &,' me=',me
+ idate(1) = idate7(4)
+ idate(2:3) = idate7(2:3)
+ idate(4) = idate7(1)
+!
+ call nemsio_getheadvar(gfile_in,'iorder',iorder,iret=iret)
+ call nemsio_getheadvar(gfile_in,'irealf',irealf,iret=iret)
+ call nemsio_getheadvar(gfile_in,'igen',igen,iret=iret)
+ call nemsio_getheadvar(gfile_in,'dimx',lonb,iret=iret)
+ call nemsio_getheadvar(gfile_in,'dimy',latb,iret=iret)
+ call nemsio_getheadvar(gfile_in,'iens',iens,iret=iret)
+ call nemsio_getheadvar(gfile_in,'idpp',idpp,iret=iret)
+ call nemsio_getheadvar(gfile_in,'idrun',idrun,iret=iret)
+ call nemsio_getheadvar(gfile_in,'idusr',idusr,iret=iret)
+
+ if (me == 0) then
+ write(0,*)'iret=',iret,
+ & ' levsi=',levsi,' idate=',idate,
+ & 'lonf=',lonf,'lonfi=',lonfi,'latg=',latg,'latgi=',latgi,
+ & 'jcap=',jcap,'jcapi=',jcapi,'levs=',levs,'levsi=',levsi,
+ & 'idvc=',idvc,'tlmeta=',tlmeta,
+ & 'gen_coord_hybrid=',gen_coord_hybrid
+ if(lonf .ne. lonfi .or. latg .ne. latgi .or.
+ & jcap .ne. jcapi .or. levs .ne. levsi) then
+ print *,' Input resolution and the model resolutions are'
+ &, ' different- run aborted'
+ call mpi_quit(555)
+ endif
+ if ( gen_coord_hybrid ) then
+ if(me==0) print *, ' Use sigma-theta-p hybrid coordinate'
+ if (idvc == 3 ) then
+ if(me==0)
+ & print *, ' Cold_start input is consistent, run continues'
+ else
+ if(me==0)
+ & print *, ' Cold_start input is different, run aborted'
+ call mpi_quit(556)
+ endif
+ endif
+ if ( hybrid ) then
+ if(me==0)print *, ' Use sigma-p hybrid coordinate'
+ if (idvc == 2 ) then
+ if(me==0)
+ & print *, ' Cold_start input is consistent, run continues'
+ else
+ if(me==0)
+ & print *, ' Cold_start input is different, run aborted'
+ call mpi_quit(557)
+ endif
+ endif
+ endif
+!
+! retrieve nreci, recnamei, reclevtypi, and reclevi
+ call nemsio_getfilehead(gfile_in,iret=iret,nrec=nreci)
+ if (me == 0) then
+ print *, 'LU_TRC: nreci =', nreci, iret
+ endif
+
+ allocate (recnamei(nreci))
+ allocate (reclevtypi(nreci))
+ allocate (reclevi(nreci))
+ call nemsio_getfilehead(gfile_in,iret=iret,recname=recnamei,
+ & reclevtyp=reclevtypi,reclev=reclevi)
+ stime=timef()
+ igen = igen
+ ienst = iens(1)
+ iensi = iens(2)
+!Weiyu, pdryini4 did not initiaqlized, temp give value zero.
+!-----------------------------------------------------------
+!
+!
+ IF (me.eq.0) THEN
+ write(0,*)'cfile,in treadeo idate=',cfile,idate
+ &, ' idvc=',idvc,' jcap=',jcap
+ ENDIF
+!
+ allocate (nemsio_data(lonf*latg))
+ stime=timef()
+
+ call w3kind(w3rlkind,w3ikind)
+!
+! Read ps
+ if(w3rlkind==8) then
+ call nemsio_readrecv(gfile_in,'pres','sfc',1,nemsio_data,
+ & iret=iret)
+ elseif(w3rlkind==4) then
+ call nemsio_readrecvw34(gfile_in,'pres','sfc',1,nemsio_data,
+ & iret=iret)
+ endif
+ call split2d(nemsio_data,buffo,global_lats_a)
+ CALL interpred(1,kmsk,buffo,psg,global_lats_a,lonsperlat)
+!
+! Read u
+ do k=1,levs
+ if(w3rlkind==8) then
+ call nemsio_readrecv(gfile_in,'ugrd','mid layer',k,nemsio_data,
+ & iret=iret)
+ elseif(w3rlkind==4) then
+ call nemsio_readrecvw34(gfile_in,'ugrd','mid layer',k,
+ & nemsio_data,iret=iret)
+ endif
+! print *,'in treadeo,ugrd=',maxval(nemsio_data),
+! & minval(nemsio_data),'iret=',iret,'k=',k
+ call split2d(nemsio_data,buffo,global_lats_a)
+ CALL interpred(1,kmsk,buffo,uug(1,1,k),global_lats_a,lonsperlat)
+ enddo
+! Read v
+ do k=1,levs
+ if(w3rlkind==8) then
+ call nemsio_readrecv(gfile_in,'vgrd','mid layer',k,nemsio_data,
+ & iret=iret)
+ elseif(w3rlkind==4) then
+ call nemsio_readrecvw34(gfile_in,'vgrd','mid layer',k,
+ & nemsio_data,iret=iret)
+ endif
+ call split2d(nemsio_data,buffo,global_lats_a)
+ CALL interpred(1,kmsk,buffo,vvg(1,1,k),global_lats_a,lonsperlat)
+ enddo
+! Read T -- this is real temperature
+ do k=1,levs
+ if(w3rlkind==8) then
+ call nemsio_readrecv(gfile_in,'tmp','mid layer',k,nemsio_data,
+ & iret=iret)
+ elseif(w3rlkind==4) then
+ call nemsio_readrecvw34(gfile_in,'tmp','mid layer',k,
+ & nemsio_data,iret=iret)
+ endif
+ call split2d(nemsio_data,buffo,global_lats_a)
+ CALL interpred(1,kmsk,buffo,ttg(1,1,k),global_lats_a,lonsperlat)
+ enddo
+
+! Initial Tracers with zero
+!
+!* rqg(:,:,:) = 0.0
+ rqg(:,:,:) = 1.e-20
+
+!! Generalized tracers:
+!! Loop through ntrac to read in met + chem tracers
+!*
+ do n = 1, ntrac
+ vname = trim(gfs_dyn_tracer%vname(n, 1))
+ if(me==0) print *,'LU_TRC: initialize ',n,vname
+ &,' w3rlkind=',w3rlkind
+ do k=1,levs
+ if(w3rlkind==8) then
+ call nemsio_readrecv(gfile_in,trim(vname),
+ & 'mid layer',k,nemsio_data,iret=iret)
+ elseif(w3rlkind==4) then
+ call nemsio_readrecvw34(gfile_in,trim(vname),
+ & 'mid layer',k,nemsio_data,iret=iret)
+ endif
+ if(iret == 0) then
+!* if(me==0) print *,'LU_TRC: tracer read in ok -',
+!* & gfs_dyn_tracer%vname(n, 1),k
+ if(me==0 .and. k==1)
+ & print *,'LU_TRC: tracer read in ok '
+! if (me == 0 .and. n == 1) then
+! write(0,*)' k=',k,' rqg=',minval(nemsio_data)
+! &,maxval(nemsio_data)
+! endif
+
+ call split2d(nemsio_data,buffo,global_lats_a)
+ CALL interpred(1,kmsk,buffo,rqg(1,1,k+(n-1)*levs),
+ & global_lats_a,lonsperlat)
+ else
+!* if(me==0) print *,'LU_TRC: tracer not found in input; ',
+!* & 'set chem tracer to default values',me,k
+ if(me==0 .and. k==1) print *,
+ & 'LU_TRC: tracer not found in input; set to default'
+ endif
+ enddo
+ enddo
+ etime = timef()
+ call nemsio_close(gfile_in,iret)
+ iprint = 0
+!!!!
+!!!! Added by PJP
+ deallocate (recnamei)
+ deallocate (reclevtypi)
+ deallocate (reclevi)
+ deallocate (nemsio_data)
+ RETURN
+ END
Index: checkout/dyn/compns_dynamics.f
===================================================================
--- checkout/dyn/compns_dynamics.f (revision 93160)
+++ checkout/dyn/compns_dynamics.f (revision 94274)
@@ -117,6 +117,7 @@
& fhout_hf,fhmax_hf,
& cont_eq_opt1,quamon,wgtm,time_extrap_etadot,iter_one_no_interp,
& opt1_3d_qcubic,dfilevs,yhalo,
+ & iau,iaufiles_fg,iaufiles_anl,iaufhrs,iau_delthrs,yhalo,iaulnp,
& sppt,sppt_tau,sppt_lscale,sppt_logit,
& vcamp,vc_lscale,vc_tau,vc_logit,
& iseed_shum,iseed_sppt,shum,shum_tau,shum_lscale,iseed_vc,
@@ -270,6 +271,11 @@
wam_ipe_coupling = .false.
height_dependent_g = .false.
!
+! iau parameters
+ iau = .false.
+ iaufhrs = -1
+ iaufiles_fg = ''
+ iaufiles_anl = ''
if (me == 0) print *,' nlunit=',nlunit,' gfs_dyn_namelist=',
& gfs_dyn_namelist
!$$$ read(5,nam_dyn)
Index: checkout/dyn/do_dynamics_slg_loop.f
===================================================================
--- checkout/dyn/do_dynamics_slg_loop.f (revision 93160)
+++ checkout/dyn/do_dynamics_slg_loop.f (revision 94274)
@@ -12,7 +12,7 @@
! & ANL_GR_A_1, ANL_GR_A_2,
& start_step,restart_step,reset_step,end_step,
& dfiend_step, pdryini,nblck,ZHOUR,
- & pwat,ptot,ptrc,nfcstdate7)
+ & pwat,ptot,ptrc,nfcstdate7,iniauinterval)
!!
use gfs_dyn_machine , only : kind_evod, kind_grid
use gfs_dyn_resol_def , only : latg,latg2,levh,levs,levp1,
@@ -68,6 +68,7 @@
INTEGER, SAVE :: ndfih, kdtdfi
CHARACTER(16) :: CFHOUR1
INTEGER,INTENT(IN), dimension(latg) :: LONSPERLAT, GLOBAL_LATS_A
+ LOGICAL,INTENT(IN) :: iniauinterval
!!
INTEGER,INTENT(IN) :: nfcstdate7(7)
REAL(KIND=KIND_EVOD), INTENT(in) :: deltim,PHOUR
@@ -337,23 +338,25 @@
endif
!----------------------------------------------------------
!
- do i=1,len_trie_ls
- sum_k_rqchange_ls(i,1) = trie_ls(i,1,p_q)
- sum_k_rqchange_ls(i,2) = trie_ls(i,2,p_q)
- enddo
- do i=1,len_trio_ls
- sum_k_rqchango_ls(i,1) = trio_ls(i,1,p_q)
- sum_k_rqchango_ls(i,2) = trio_ls(i,2,p_q)
- enddo
-
- do i=1,len_trie_ls
- trie_ls(i,1,p_q) = save_qe_ls(i,1)
- trie_ls(i,2,p_q) = save_qe_ls(i,2)
- enddo
- do i=1,len_trio_ls
- trio_ls(i,1,p_q) = save_qo_ls(i,1)
- trio_ls(i,2,p_q) = save_qo_ls(i,2)
- enddo
+ if (.not. iniauinterval) then ! skip if in IAU forcing interval
+ do i=1,len_trie_ls
+ sum_k_rqchange_ls(i,1) = trie_ls(i,1,p_q)
+ sum_k_rqchange_ls(i,2) = trie_ls(i,2,p_q)
+ enddo
+ do i=1,len_trio_ls
+ sum_k_rqchango_ls(i,1) = trio_ls(i,1,p_q)
+ sum_k_rqchango_ls(i,2) = trio_ls(i,2,p_q)
+ enddo
+
+ do i=1,len_trie_ls
+ trie_ls(i,1,p_q) = save_qe_ls(i,1)
+ trie_ls(i,2,p_q) = save_qe_ls(i,2)
+ enddo
+ do i=1,len_trio_ls
+ trio_ls(i,1,p_q) = save_qo_ls(i,1)
+ trio_ls(i,2,p_q) = save_qo_ls(i,2)
+ enddo
+ endif
!
!$omp parallel do private(k,kp_d,kp_z,kp_t,kp_x,kp_w,kp_y,j)
@@ -450,7 +453,8 @@
!
call do_dynamics_gridn2anl_slg(grid_gr,anl_gr_a_2
&, rcs2_a
- &, global_lats_a,lonsperlat)
+ &, global_lats_a,lonsperlat
+ &, iniauinterval)
!
! transform values in grid to spectral
!
@@ -462,6 +466,7 @@
!
!----------------------------------------------------------
!
+ if (.not. iniauinterval) then ! skip if in IAU forcing interval
do i=1,len_trie_ls
sum_k_rqchange_ls(i,1) = trie_ls(i,1,p_q)
sum_k_rqchange_ls(i,2) = trie_ls(i,2,p_q)
@@ -479,6 +484,7 @@
trio_ls(i,1,p_q) = save_qo_ls(i,1)
trio_ls(i,2,p_q) = save_qo_ls(i,2)
enddo
+ endif
! print *,' ----- do bfilter ------ '
@@ -827,6 +833,7 @@
!
! testing mass correction on sep 25
!!
+ if (.not. iniauinterval) then
if(gg_tracers)then
do i=1,len_trie_ls
trie_ls(i,1,p_q) = trie_ls(i,1,p_q)
@@ -856,6 +863,7 @@
enddo
enddo
endif
+ endif
! testing mass correction on sep 25
!
!$omp parallel do private(k,i,item,jtem,ktem,ltem,mtem)
Index: checkout/dyn/do_dynamics_mod.f
===================================================================
--- checkout/dyn/do_dynamics_mod.f (revision 93160)
+++ checkout/dyn/do_dynamics_mod.f (revision 94274)
@@ -281,7 +281,8 @@
! --------------------------------------------------------------
subroutine do_dynamics_gridn2anl_slg(grid_gr,anl_gr_a_2,
& rcs2,
- & global_lats_a,lonsperlat)
+ & global_lats_a,lonsperlat
+ & ,iniauinterval)
real(kind=kind_grid) grid_gr(lonf*lats_node_a_max,lotgr)
! real(kind=kind_evod) anl_gr_a_2(lonfx*lota,lats_dim_ext)
@@ -289,6 +290,7 @@
real(kind=kind_grid) rcs2(latg2)
integer,intent(in):: global_lats_a(latg)
integer,intent(in):: lonsperlat(latg)
+ logical,intent(in) :: iniauinterval
integer jlonf, ilan, kap, kau, kav, kat, kar
&, lan,lat,lon_dim,lons_lat,k,i,kk
real rcs_loc
@@ -319,21 +321,20 @@
enddo
enddo
!
- if (gg_tracers) then
-!$omp parallel do private(i,ilan)
+ if (gg_tracers .and. .not. iniauinterval) then
do i=1,lons_lat
ilan = i+jlonf
anl_gr_a_2(i,kap,lan) = grid_gr(ilan,g_rqtk)
enddo
else
-!$omp parallel do private(i)
do i=1,lons_lat
- anl_gr_a_2(i,kap,lan) = 0.0
+! anl_gr_a_2(i,kap,lan) = 0.0
+ ilan = i + jlonf
+ anl_gr_a_2(i,kap,lan) = grid_gr(ilan,g_zqp) ! surface pressure with IAU increment
enddo
endif
!
if( .not. ndslfv) then
-!$omp parallel do private(i,k,ilan,kk)
do k=1,levh
kk = kar-1+k
do i=1,lons_lat
Index: checkout/dyn/gfs_dynamics_internal_state_mod.f
===================================================================
--- checkout/dyn/gfs_dynamics_internal_state_mod.f (revision 93160)
+++ checkout/dyn/gfs_dynamics_internal_state_mod.f (revision 94274)
@@ -26,6 +26,7 @@
! Feb 2011 sarah lu, add thermodyn_id, sfcpress_id
! Sep 2011 weiyu yang, modified for using the ESMF 5.2.0r library.
! Sep 2012 jun wang, add nemsio_in to specify input files
+! Aug 2016 P Pegion add iniauinterval
!
! !interface:
!
@@ -76,7 +77,7 @@
character(32) ,allocatable :: filename_base(:)
integer :: ipt_lats_node_a
integer :: lats_node_a
- logical :: adiabatic
+ logical :: adiabatic,iniauinterval
logical :: slg_flag
!jwe
Index: checkout/dyn/wrtout_dynamics.f
===================================================================
--- checkout/dyn/wrtout_dynamics.f (revision 93160)
+++ checkout/dyn/wrtout_dynamics.f (revision 94274)
@@ -106,6 +106,10 @@
!!
character CFHOUR*40,CFORM*40,filename*255
integer jdate(4),nzsig,ndigyr,ndig,kh,ioproc
+! --- locals for iau:
+ integer :: jdate_iau(4),idat(8),jdat(8)
+ real (kind=kind_evod) :: zhour_iau,fhour_iau
+ real :: rinc(5)
!!
REAL (KIND=KIND_grid) pdryini
INTEGER, dimension(latg) :: GLOBAL_lats_a, lonsperlat
@@ -177,6 +181,28 @@
JDATE(4) = MOD(IDATE(4)-1,100) + 1
ENDIF
+ if ( iau .and. fhour >= 6.0 ) then
+ idat = 0
+ idat(1) = idate(4)
+ idat(2) = idate(2)
+ idat(3) = idate(3)
+ idat(5) = idate(1)
+ rinc = 0.
+ rinc(2) = 6.
+ call w3movdat(rinc,idat,jdat)
+ jdate_iau(4) = jdat(1)
+ jdate_iau(2) = jdat(2)
+ jdate_iau(3) = jdat(3)
+ jdate_iau(1) = jdat(5)
+ fhour_iau = fhour - 6.0
+ zhour_iau = zhour - 6.0
+ if ( zhour_iau < 0.0 ) zhour_iau = 0.0
+ else
+ jdate_iau = jdate
+ fhour_iau = fhour
+ zhour_iau = zhour
+ endif
+
!sela set lfnhr to false for writing one step output etc.
lfnhr = .true. ! no output
! lfnhr=3600*abs(fhour-nint(fhour)).le.1.or.phour.eq.0
@@ -406,10 +432,17 @@
if( .not. ndslfv ) then
- call twrites_hst(filename,ioproc,fhour,idate,
- & ls_nodes,max_ls_nodes,trie_ls,trio_ls,
- & trie_ls,trio_ls,
- & thermodyn_id_out,sfcpress_id_out,pdryini)
+ if(iau) then
+ call twrites_hst(filename,ioproc,fhour,jdate_iau,
+ & ls_nodes,max_ls_nodes,trie_ls,trio_ls,
+ & trie_ls,trio_ls,
+ & thermodyn_id_out,sfcpress_id_out,pdryini)
+ else
+ call twrites_hst(filename,ioproc,fhour,idate,
+ & ls_nodes,max_ls_nodes,trie_ls,trio_ls,
+ & trie_ls,trio_ls,
+ & thermodyn_id_out,sfcpress_id_out,pdryini)
+ end if
else
Index: checkout/dyn/spect_to_grid_iau.f
===================================================================
--- checkout/dyn/spect_to_grid_iau.f (nonexistent)
+++ checkout/dyn/spect_to_grid_iau.f (revision 94274)
@@ -0,0 +1,57 @@
+ subroutine spect_to_grid_iau(
+ & trie_ls,trio_ls,datag,
+ & ls_node,ls_nodes,max_ls_nodes,
+ & lats_nodes_a,global_lats_a,lonsperlar,
+ & plnev_a,plnod_a,nlevs,lon1,lon2)
+
+ use gfs_dyn_resol_def
+ use gfs_dyn_layout1
+ use gfs_dyn_vert_def
+ use gfs_dyn_coordinate_def
+ use gfs_dyn_mpi_def
+ implicit none
+ real(kind=kind_evod), intent(in) :: trie_ls(len_trie_ls,2,nlevs)
+ real(kind=kind_evod), intent(in) :: trio_ls(len_trio_ls,2,nlevs)
+ real(kind=kind_phys), intent(out) :: datag(lon2,lats_node_a,nlevs)
+ integer, intent(in) :: ls_node(ls_dim,3),ls_nodes(ls_dim,nodes),
+ & nlevs,max_ls_nodes(nodes),lats_nodes_a(nodes),
+ & global_lats_a(latg),lonsperlar(latg),lon1,lon2
+ real(kind=kind_evod),intent(in) :: plnev_a(len_trie_ls,latg2),
+ & plnod_a(len_trio_ls,latg2)
+! local vars
+ real(kind=kind_evod) for_gr_a_1(lon1,nlevs,lats_dim_a)
+ real(kind=kind_evod) for_gr_a_2(lon2,nlevs,lats_dim_a)
+ integer i,j,k
+ integer l,lan,lat
+ integer lons_lat
+
+ call sumfln_slg_gg(trie_ls,
+ & trio_ls,
+ & lat1s_a,
+ & plnev_a,plnod_a,
+ & nlevs,ls_node,latg2,
+ & lats_dim_a,nlevs,for_gr_a_1,
+ & ls_nodes,max_ls_nodes,
+ & lats_nodes_a,global_lats_a,
+ & lats_node_a,ipt_lats_node_a,lon1,
+ & lonsperlar,lon1,latg,0)
+
+ do lan=1,lats_node_a
+ lat = global_lats_a(ipt_lats_node_a-1+lan)
+ lons_lat = lonsperlar(lat)
+ CALL FOUR_TO_GRID(for_gr_a_1(1,1,lan),for_gr_a_2(1,1,lan),
+ & lon1,lon2,lons_lat,nlevs)
+ enddo
+
+ datag = 0.
+ do lan=1,lats_node_a
+ lat = global_lats_a(ipt_lats_node_a-1+lan)
+ lons_lat = lonsperlar(lat)
+ do k=1,nlevs
+ do i=1,lons_lat
+ datag(i,lan,k) = for_gr_a_2(i,k,lan)
+ enddo
+ enddo
+ enddo
+ end subroutine spect_to_grid_iau
+
Index: checkout/dyn/treadeo.io_iau.f
===================================================================
--- checkout/dyn/treadeo.io_iau.f (nonexistent)
+++ checkout/dyn/treadeo.io_iau.f (revision 94274)
@@ -0,0 +1,319 @@
+ SUBROUTINE TREADEO_IAU(IDATE,grid_gr,
+ X LS_NODE,LS_NODES,MAX_LS_NODES,
+ X IPRINT,
+ & global_lats_a,lats_nodes_a,lonsperlat, cfile,
+ & epse, epso, epsedn,epsodn, plnew_a, plnow_a,
+ & plnev_a, plnod_a,snnp1ev, snnp1od)
+!
+!-- revision history
+!
+! Sep 2012, J. Wang: change interface for nems gfs
+! Oct 2012, H. Juang: recover terrain laplacian preparation
+! May 2013 S. Moorthi - correct definition of idvt for WAM
+! Jun 2013, H. Juang: add dp and spectral moisture in case of ndsl
+! Sep 2013, H. Juang: add computation of dp, not only give dp space
+!
+ use gfs_dyn_resol_def
+ use gfs_dyn_layout1
+ use gfs_dyn_coordinate_def ! hmhj
+! use sig_io
+ use gfs_dyn_io_header
+ use gfs_dyn_iau_module, only: gq_iau,gu_iau,gv_iau,gt_iau,grq_iau
+ use namelist_dynamics_def
+! use gfs_dyn_namelist_def
+ use gfs_dyn_vert_def
+ use gfs_dyn_mpi_def
+ use gfs_dyn_physcons, rerth => con_rerth, grav => con_g,
+ & rkap => con_rocp, cpd => con_cp
+ USE gfs_dyn_date_def, ONLY: FHOUR
+ use gfs_dyn_tracer_config, ONLY : gfs_dyn_tracer,glbsum ! generalized tracer
+ use sigio_module
+ use sigio_r_module
+
+!
+ IMPLICIT NONE
+!
+ INTEGER,intent(inout) :: IDATE(4)
+ REAL(KIND=KIND_GRID),dimension(lonf,lats_node_a_max,gq_iau),
+ & intent(out) :: grid_gr
+ INTEGER,intent(in) :: LS_NODE(LS_DIM,3)
+ INTEGER,intent(in) :: LS_NODES(LS_DIM,NODES)
+ INTEGER,intent(in) :: MAX_LS_NODES(NODES)
+ integer,intent(in) :: lats_nodes_a(nodes)
+ INTEGER,intent(inout) :: IPRINT
+ integer,intent(in) :: global_lats_a(latg),lonsperlat(latg)
+ character*(*),intent(in) :: cfile
+ real(kind=kind_evod) :: trie_ls(len_trie_ls,2,gq_iau)
+ real(kind=kind_evod) :: trio_ls(len_trio_ls,2,gq_iau)
+ real(kind=kind_evod),dimension(len_trie_ls),intent(in) :: epse,
+ & epsedn
+ real(kind=kind_evod),dimension(len_trio_ls),intent(in) :: epso,
+ & epsodn
+!
+ real(kind=kind_evod),dimension(len_trie_ls,latg2),intent(in)
+ & :: plnew_a,plnev_a
+ real(kind=kind_evod),dimension(len_trio_ls,latg2),intent(in)
+ & :: plnow_a,plnod_a
+!
+ REAL(KIND=KIND_EVOD),intent(in) :: snnp1ev(LEN_TRIE_LS)
+ REAL(KIND=KIND_EVOD),intent(in) :: snnp1od(LEN_TRIO_LS)
+ integer lats_nodes_r(nodes)
+!
+!------local vars
+!
+ REAL(KIND=KIND_EVOD) QE(LEN_TRIE_LS,2,1)
+ &, TEE(LEN_TRIE_LS,2,LEVS)
+ &, DIE(LEN_TRIE_LS,2,LEVS)
+ &, ZEE(LEN_TRIE_LS,2,LEVS)
+ &, UUE(LEN_TRIE_LS,2,LEVS)
+ &, VVE(LEN_TRIE_LS,2,LEVS)
+ &, RQE(LEN_TRIE_LS,2,LEVS)
+ &, QO(LEN_TRIO_LS,2,1)
+ &, TEO(LEN_TRIO_LS,2,LEVS)
+ &, DIO(LEN_TRIO_LS,2,LEVS)
+ &, ZEO(LEN_TRIO_LS,2,LEVS)
+ &, UUO(LEN_TRIO_LS,2,LEVS)
+ &, VVO(LEN_TRIO_LS,2,LEVS)
+ &, RQO(LEN_TRIO_LS,2,LEVS)
+ &, Z00
+!
+ real(kind=kind_evod), allocatable :: syn_gr_a_2(:,:,:)
+!
+!
+!cmr ls_node(1,1) ... ls_node(ls_max_node,1) : values of L
+!cmr ls_node(1,2) ... ls_node(ls_max_node,2) : values of jbasev
+!cmr ls_node(1,3) ... ls_node(ls_max_node,3) : values of jbasod
+!
+ INTEGER J,K,L,LOCL,N,kk,ilan,itrac
+ integer i,lan,lat,lons_lat,lon,njeff,nn,lon_dim
+ integer indev
+ integer indod
+ integer indev1,indev2
+ integer indod1,indod2
+ REAL(KIND=KIND_EVOD) WAVES
+ REAL(KIND=KIND_EVOD), target :: TRISCA(LNT2)
+ REAL(KIND=KIND_EVOD) sikp1(levp1)
+ INTEGER INDLSEV,JBASEV
+ INTEGER INDLSOD,JBASOD
+c$$$ REAL(KIND=KIND_IO4) Z(lnt2)
+!
+ type(sigio_head) head
+ type(sigio_dbti) dati
+!
+ integer iret, num_dta,nft,lon1,lon2
+ real(kind=kind_evod) psurfff
+ real(kind=kind_evod) pressk, tem, rkapi, rkapp1
+!
+ real(kind=kind_evod) teref(levp1),ck5p(levp1) ! hmhj
+
+
+ INCLUDE 'function2'
+!!
+!!
+ print *,' enter treadeo.io_fd ' ! hmhj
+ if (semilag) then
+ lon1=lon_dim_a
+ lon2=lonf
+ else
+ lon1=lonfx
+ lon2=lonfx
+ endif
+ allocate(syn_gr_a_2(lon2,lats_node_a,levs))
+
+ nft=152
+ call sigio_rropen(nft,cfile,iret)
+ call sigio_alhead(head,iret)
+ call sigio_rrhead(nft,head,iret)
+!
+ ivsinp = head%ivs
+ IF (me .eq. 0) THEN
+ print *,' In treadeo iret=',iret,' cfile=',cfile
+ &,' ivs=',head%ivs,' levs=',head%levs
+ ENDIF
+ IF (iret .ne. 0) THEN
+ print *,' unable to read from unit ',nft,' Job Aborted'
+ &,' iret=',iret,' me=',me
+ call mpi_quit(7777)
+ ENDIF
+!
+ idvc = head%idvc !idvc=3:sigma-theta and/or p, idvc=2:sigma-p, idvc=1:sigma files
+ idsl = head%idsl
+!
+ dati%i = 2 ! Surface pressure
+ dati%f => TRISCA
+ call sigio_rrdbti(nft,head,dati,iret)
+! IF (me == 0) print *,' SFCPRES=',trisca(1:10)
+ CALL TRISEORI(TRISCA,QE(1,1,1),QO(1,1,1),1,LS_NODE)
+ DO K=1,LEVS
+ dati%i = k + 2 ! Virtual Temperature or CpT
+ dati%f => TRISCA
+
+ call sigio_rrdbti(nft,head,dati,iret)
+
+ CALL TRISEORI(TRISCA,TEE(1,1,K),TEO(1,1,K),1,LS_NODE)
+ ENDDO
+!
+ DO K=1,LEVS
+ dati%i = levs + 2 + (k-1) * 2 + 1 ! Divergence
+ dati%f => TRISCA
+ call sigio_rrdbti(nft,head,dati,iret)
+ CALL TRISEORI(TRISCA,DIE(1,1,K),DIO(1,1,K),1,LS_NODE)
+!
+ dati%i = levs + 2 + (k-1) * 2 + 2 ! Vorticity
+ dati%f => TRISCA
+ call sigio_rrdbti(nft,head,dati,iret)
+ CALL TRISEORI(TRISCA,ZEE(1,1,K),ZEO(1,1,K),1,LS_NODE)
+ END DO
+csela print*,' levh=',levh
+!
+!
+!
+ DO itrac=1,ntrac
+ RQE=0.
+ RQO=0.
+ kk = 0
+ IF (itrac == 1) THEN
+ kk = 1
+ ELSEIF (itrac == ntoz) THEN
+ kk = ntoz
+ ELSEIF (itrac >= ntcw .and. itrac < ntcw+ncld-1) THEN
+ DO n=1,ncld
+ IF (itrac == ntcw+n-1) kk = ntcw+n-1
+ ENDDO
+ ELSE
+ kk = itrac
+ ENDIF
+!
+ DO k=1,levs
+ dati%i = levs * (2+itrac) + 2 + k ! Tracers starting with q
+ dati%f => TRISCA
+ call sigio_rrdbti(nft,head,dati,iret)
+ CALL TRISEORI(TRISCA,RQE(1,1,k),RQO(1,1,k),1,LS_NODE)
+ END DO
+ call spect_to_grid_iau(RQE,RQO,
+ & syn_gr_a_2,
+ & ls_node,ls_nodes,max_ls_nodes,
+ & lats_nodes_a,global_lats_a,lonsperlat,
+ & plnev_a,plnod_a,levs,lon1,lon2)
+!jw!$omp parallel do private(lan)
+!jw!$omp+private(lat,lon_dim,lons_lat,i,itrac)
+ DO lan=1,lats_node_a
+ lat = global_lats_a(ipt_lats_node_a-1+lan)
+ lon_dim = lon_dims_a(lan)
+ lons_lat = lonsperlat(lat)
+ DO k=1,levs
+ DO i=1,lons_lat
+ grid_gr(i,lan,grq_iau+(kk-1)*levs+k-1)=syn_gr_a_2(i,lan,k) ! make sure I have them in the right slots
+ ENDDO
+ ENDDO
+ ENDDO
+
+ END DO
+
+!
+ DO k=1,levs
+ call dezouv(DIE(1,1,k), ZEO(1,1,k),
+ x UUE(1,1,k), VVO(1,1,k),
+ x epsedn,epsodn,snnp1ev,snnp1od,ls_node)
+!
+ call dozeuv(DIO(1,1,k), ZEE(1,1,k),
+ x UUO(1,1,k), VVE(1,1,k),
+ x epsedn,epsodn,snnp1ev,snnp1od,ls_node)
+ ENDDO
+ call spect_to_grid_iau(QE,QO,
+ & syn_gr_a_2(:,:,1),
+ & ls_node,ls_nodes,max_ls_nodes,
+ & lats_nodes_a,global_lats_a,lonsperlat,
+ & plnev_a,plnod_a,1,lon1,lon2)
+!jw!$omp parallel do private(lan)
+!jw!$omp+private(lat,lon_dim,lons_lat,i)
+ DO lan=1,lats_node_a
+ lat = global_lats_a(ipt_lats_node_a-1+lan)
+ lon_dim = lon_dims_a(lan)
+ lons_lat = lonsperlat(lat)
+ DO i=1,lons_lat
+ grid_gr(i,lan,gq_iau)=syn_gr_a_2(i,lan,1)
+ ENDDO
+ ENDDO
+
+ call spect_to_grid_iau(UUE,UUO,
+ & syn_gr_a_2,
+ & ls_node,ls_nodes,max_ls_nodes,
+ & lats_nodes_a,global_lats_a,lonsperlat,
+ & plnev_a,plnod_a,levs,lon1,lon2)
+!jw!$omp parallel do private(lan)
+!jw!$omp+private(lat,lon_dim,lons_lat,i,k)
+ DO lan=1,lats_node_a
+ lat = global_lats_a(ipt_lats_node_a-1+lan)
+ lon_dim = lon_dims_a(lan)
+ lons_lat = lonsperlat(lat)
+ DO k=1,levs
+ DO i=1,lons_lat
+ grid_gr(i,lan,gu_iau+k-1)=syn_gr_a_2(i,lan,k)
+ ENDDO
+ ENDDO
+ ENDDO
+
+ call spect_to_grid_iau(VVE,VVO,
+ & syn_gr_a_2,
+ & ls_node,ls_nodes,max_ls_nodes,
+ & lats_nodes_a,global_lats_a,lonsperlat,
+ & plnev_a,plnod_a,levs,lon1,lon2)
+!jw!$omp parallel do private(lan)
+!jw!$omp+private(lat,lon_dim,lons_lat,i,k)
+ DO lan=1,lats_node_a
+ lat = global_lats_a(ipt_lats_node_a-1+lan)
+ lon_dim = lon_dims_a(lan)
+ lons_lat = lonsperlat(lat)
+ DO k=1,levs
+ DO i=1,lons_lat
+ grid_gr(i,lan,gv_iau+k-1)=syn_gr_a_2(i,lan,k)
+ ENDDO
+ ENDDO
+ ENDDO
+
+ call spect_to_grid_iau(TEE,TEO,
+ & syn_gr_a_2,
+ & ls_node,ls_nodes,max_ls_nodes,
+ & lats_nodes_a,global_lats_a,lonsperlat,
+ & plnev_a,plnod_a,levs,lon1,lon2)
+!jw!$omp parallel do private(lan)
+!jw!$omp+private(lat,lon_dim,lons_lat,i,k)
+ DO lan=1,lats_node_a
+ lat = global_lats_a(ipt_lats_node_a-1+lan)
+ lon_dim = lon_dims_a(lan)
+ lons_lat = lonsperlat(lat)
+ DO k=1,levs
+ DO i=1,lons_lat
+ grid_gr(i,lan,gt_iau+k-1)=syn_gr_a_2(i,lan,k)
+ ENDDO
+ ENDDO
+ ENDDO
+
+! IF(me==0) THEN
+! open (991,file='ps_grid',form='unformatted')
+! write(991)real(grid_gr(1:lonf,1:lats_node_a,gq_iau),kind=4)
+! close(991)
+! ENDIF
+! call mpi_barrier(MPI_COMM_ALL,i)
+! call mpi_abort(mpi_comm_all,j,i)
+! print *,'in treadeo.io, ps=',
+! & maxval(grid_gr(:,1:lats_node_a,gq_iau)),
+! & minval(grid_gr(:,1:lats_node_a,gq_iau))
+! print*,'in treadeo.io, t=',
+! & maxval(grid_gr(:,1:lats_node_a,gt_iau:gt_iau+levs-1)),
+! & minval(grid_gr(:,1:lats_node_a,gt_iau:gt_iau+levs-1))
+! print*,'in treadeo.io, tracer=',
+! & maxval(grid_gr(:,1:lats_node_a,grq_iau:grq_iau+levh-1)),
+! & minval(grid_gr(:,1:lats_node_a,grq_iau:grq_iau+levh-1))
+! print*,'in treadeo.io, u=',
+! & maxval(grid_gr(:,1:lats_node_a,gu_iau:gu_iau+levs-1)),
+! & minval(grid_gr(:,1:lats_node_a,gu_iau:gu_iau+levs-1))
+! print*,'in treadeo.io, v=',
+! & maxval(grid_gr(:,1:lats_node_a,gv_iau:gv_iau+levs-1)),
+! & minval(grid_gr(:,1:lats_node_a,gv_iau:gv_iau+levs-1))
+!
+! print *,' leave treadeo.io_fd '
+
+ RETURN
+ END
Property changes on: checkout/dyn/treadeo.io_iau.f
___________________________________________________________________
Added: svn:executable
## -0,0 +1 ##
+*
\ No newline at end of property
Index: checkout/dyn/makefile
===================================================================
--- checkout/dyn/makefile (revision 93160)
+++ checkout/dyn/makefile (revision 94274)
@@ -75,6 +75,7 @@
gfs_dyn_bfilt_def.o \
gfs_dyn_glats.o \
gfs_dyn_dfi_mod.o \
+ gfs_dyn_iau_mod.o \
nemsio_def.o \
wrtout_dynamics.o \
wrtout_dynamics_slg_nemsio.o \
@@ -163,6 +164,7 @@
grid_to_spect_rqt.o \
gridzz_to_spect.o \
spect_to_grid.o \
+ spect_to_grid_iau.o \
spect_to_grid_slg.o \
spect_to_grid_gz.o \
spect_to_grid_rqt.o \
@@ -224,12 +226,14 @@
sigio_r_module.o \
treadeo.io.o \
treadeo.io_slg.o \
+ treadeo.io_iau.o \
treadeo_nemsio.o \
treadeo_nemsio_slg.o \
treadg_nemsio.o \
treads_nemsio.o \
+ treadeo_nemsio_iau.o \
twrites_rst.o \
- twrites_rst_idea.o \
+ twrites_rst_idea.o \
twrites_hst.o \
twriteg_rst.o \
grid_to_spect_inp.o \
@@ -353,6 +357,9 @@
gfs_dynamics_grid_comp_mod.o: gfs_dynamics_grid_comp_mod.f90
$(FC) $(FFLAG90) -c gfs_dynamics_grid_comp_mod.f90
+gfs_dyn_iau_mod.o: gfs_dyn_iau_mod.f90
+ $(FC) $(FFLAG90) -c gfs_dyn_iau_mod.f90
+
#
#
#
Index: checkout/dyn/gfs_dynamics_initialize_mod.f
===================================================================
--- checkout/dyn/gfs_dynamics_initialize_mod.f (revision 93160)
+++ checkout/dyn/gfs_dynamics_initialize_mod.f (revision 94274)
@@ -33,6 +33,7 @@
! Jan 20 2013 J. Wang idea diffusion init interface change
! Jun 26 2014 S. Moorthi modified to read lonsperlat from a file
! Jul 21 2014 S. Moorthi removed num_reduce ; some cleanup
+! Aug 17 2016 P. Pegion add call for iau update init
!
!
! !interface:
@@ -48,11 +49,13 @@
use gfs_dyn_write_state, only : buff_mult_pieceg
use gfs_dyn_layout1, only : ipt_lats_node_a, lats_node_a_max
use gfs_dyn_resol_def, only : adiabatic, thermodyn_id, sfcpress_id
- use namelist_dynamics_def, only : fhrot,fhini,nemsio_in,do_shum,do_sppt,do_skeb,do_vc
+ use namelist_dynamics_def, only : fhrot,fhini,nemsio_in,do_shum,do_sppt,do_skeb,do_vc,iau
use gfs_dyn_tracer_config, only: gfs_dyn_tracer, tracer_config_init,gfs_dyn_tracer
use gfs_dyn_io_header, only: z_r,z
! stochastic perturbations
use gfs_dyn_stoch_data, only: init_stochdata
+! iau forcing
+ use gfs_dyn_iau_module, only: init_iau
!#ifndef IBM
! USE omp_lib
!#endif
@@ -782,6 +785,9 @@
endif
endif
!!
+ if (iau) then
+ call init_iau(gis_dyn)
+ endif
gis_dyn%zhour = fhour
rc = 0
!
More information about the WAM-IPE
mailing list