diff --git a/src/components/data_comps/datm/datm_comp_mod.F90 b/src/components/data_comps/datm/datm_comp_mod.F90
index a2af1a37cee..af768598810 100644
--- a/src/components/data_comps/datm/datm_comp_mod.F90
+++ b/src/components/data_comps/datm/datm_comp_mod.F90
@@ -77,6 +77,7 @@ module datm_comp_mod
integer(IN) :: krc,krl,ksc,ksl,kswndr,kswndf,kswvdr,kswvdf,kswnet
integer(IN) :: kanidr,kanidf,kavsdr,kavsdf
integer(IN) :: stbot,swind,sz,spbot,sshum,stdew,srh,slwdn,sswdn,sswdndf,sswdndr
+ integer(IN) :: stbot_noise, sprecn_noise, sswdn_noise, slwdn_noise
integer(IN) :: sprecc,sprecl,sprecn,sco2p,sco2d,sswup,sprec,starcf
! water isotopes / tracer input
@@ -161,7 +162,7 @@ module datm_comp_mod
! other fields used in calculations. Fields that are simply read and passed directly to
! the coupler do not need to be in these lists.
- integer(IN),parameter :: ktranss = 33
+ integer(IN),parameter :: ktranss = 37
character(16),parameter :: stofld(1:ktranss) = &
(/"strm_tbot ","strm_wind ","strm_z ","strm_pbot ", &
@@ -174,7 +175,9 @@ module datm_comp_mod
"strm_prec_af ","strm_u_af ","strm_v_af ","strm_tbot_af ", &
"strm_pbot_af ","strm_shum_af ","strm_swdn_af ","strm_lwdn_af ", &
"strm_rh_18O ","strm_rh_HDO ", &
- "strm_precn_16O ","strm_precn_18O ","strm_precn_HDO " &
+ "strm_precn_16O ","strm_precn_18O ","strm_precn_HDO ", &
+ ! add perturbations (Yorck)
+ "strm_tbot_noise ","strm_precn_noise","strm_swdn_noise ","strm_lwdn_noise " &
/)
character(16),parameter :: stifld(1:ktranss) = &
@@ -189,7 +192,9 @@ module datm_comp_mod
"pbot_af ","shum_af ","swdn_af ","lwdn_af ", &
! isotopic forcing
"rh_18O ","rh_HDO ", &
- "precn_16O ","precn_18O ","precn_HDO " &
+ "precn_16O ","precn_18O ","precn_HDO ", &
+ ! add perturbations (Yorck)
+ "tbot_noise ","precn_noise ","swdn_noise ","lwdn_noise " &
/)
character(CL), pointer :: ilist_av(:) ! input list for translation
@@ -444,6 +449,21 @@ subroutine datm_comp_init(Eclock, x2a, a2x, &
sprec = mct_aVect_indexRA(avstrm,'strm_prec' ,perrWith='quiet')
starcf = mct_aVect_indexRA(avstrm,'strm_tarcf' ,perrWith='quiet')
+
+ ! Yorck changes: when having an ensemble run, the memory consumption is too large if the
+ ! forcing data is perturbed one by one so that I get a file for each month for each member
+ ! idea: perturb the forcings in the CLM sourcecode with a noise file.
+ ! for each perturbed variable (temperature, precipitation, longwave and shortwave radiation),
+ ! a stream is introduced (or one stream for all)
+
+ stbot_noise = mct_aVect_indexRA(avstrm,'strm_tbot_noise' ,perrWith='quiet')
+ sprecn_noise = mct_aVect_indexRA(avstrm,'strm_precn_noise' ,perrWith='quiet')
+ sswdn_noise = mct_aVect_indexRA(avstrm,'strm_swdn_noise' ,perrWith='quiet')
+ slwdn_noise = mct_aVect_indexRA(avstrm,'strm_lwdn_noise' ,perrWith='quiet')
+
+
+
+
! anomaly forcing
sprecsf = mct_aVect_indexRA(avstrm,'strm_precsf' ,perrWith='quiet')
sprec_af = mct_aVect_indexRA(avstrm,'strm_prec_af' ,perrWith='quiet')
@@ -824,11 +844,11 @@ subroutine datm_comp_run(EClock, x2a, a2x, &
enddo ! lsize
- case('CORE_IAF_JRA', 'CORE_RYF_JRA')
+ case('CORE_IAF_JRA')
if (firstcall) then
if (sprec < 1 .or. sswdn < 1) then
- write(logunit,F00) 'ERROR: prec and swdn must be in streams for IAF/RYF JRA'
- call shr_sys_abort(trim(subname)//'ERROR: prec and swdn must be in streams for IAF/RYF JRA')
+ write(logunit,F00) 'ERROR: prec and swdn must be in streams for CORE_IAF_JRA'
+ call shr_sys_abort(trim(subname)//'ERROR: prec and swdn must be in streams for CORE_IAF_JRA')
endif
if (trim(factorFn) == 'null') then
windFactor = 1.0_R8
@@ -900,8 +920,14 @@ subroutine datm_comp_run(EClock, x2a, a2x, &
rtmp = maxval(avstrm%rAttr(stdew,:))
call shr_mpi_max(rtmp,tdewmax,mpicom,'datm_tdew',all=.true.)
endif
- if (my_task == master_task) &
- write(logunit,*) trim(subname),' max values = ',tbotmax,tdewmax,anidrmax
+ if (my_task == master_task) then
+ write(logunit,*) trim(subname),' max values = ',tbotmax,tdewmax,anidrmax
+ if (stbot_noise > 1 ) then
+ write(logunit,*) trim(subname),' Using noise from files '
+ end if
+ end if
+
+
endif
lsize = mct_avect_lsize(a2x)
do n = 1,lsize
@@ -910,8 +936,36 @@ subroutine datm_comp_run(EClock, x2a, a2x, &
!--- temperature ---
if (tbotmax < 50.0_R8) a2x%rAttr(ktbot,n) = a2x%rAttr(ktbot,n) + tkFrz
+
+ ! Temperature perturbations (Yorck)
+ if (stbot_noise > 0) then
+ if (my_task == master_task .and. n==1) then
+ !write(logunit,*) trim(subname),' temperature before: ', a2x%rAttr(ktbot,n)
+ !write(logunit,*) trim(subname),' perturbation: ', avstrm%rAttr(stbot_noise,n)
+ end if
+ a2x%rAttr(ktbot,n) = a2x%rAttr(ktbot,n) + avstrm%rAttr(stbot_noise,n)
+ if (my_task == master_task .and. n==1) then
+ !write(logunit,*) trim(subname),' temperature after: ', a2x%rAttr(ktbot,n)
+ end if
+ end if
+
+
! Limit very cold forcing to 180K
a2x%rAttr(ktbot,n) = max(180._r8, a2x%rAttr(ktbot,n))
+
+
+ ! already here perturbation of lwdn as I do not understand where this is used in the code (Yorck)
+ if (slwdn_noise > 0) then
+ if (my_task == master_task .and. n==1) then
+ !write(logunit,*) trim(subname),' lwdn before: ', avstrm%rAttr(slwdn,n)
+ !write(logunit,*) trim(subname),' perturbation: ', avstrm%rAttr(slwdn_noise,n)
+ end if
+ avstrm%rAttr(slwdn,n) = avstrm%rAttr(slwdn,n) + avstrm%rAttr(slwdn_noise,n)
+ if (my_task == master_task .and. n==1) then
+ !write(logunit,*) trim(subname),' lwdn after: ', avstrm%rAttr(slwdn,n)
+ end if
+ end if
+
a2x%rAttr(kptem,n) = a2x%rAttr(ktbot,n)
!--- pressure ---
@@ -973,18 +1027,44 @@ subroutine datm_comp_run(EClock, x2a, a2x, &
! relationship between incoming NIR or VIS radiation and ratio of
! direct to diffuse radiation calculated based on one year's worth of
! hourly CAM output from CAM version cam3_5_55
- swndr = avstrm%rAttr(sswdn,n) * 0.50_R8
+
+ ! add noise data --> multiplicative: (Yorck)
+ if (sswdn_noise>0) then
+ swndr = avstrm%rAttr(sswdn,n) * 0.50_R8 * avstrm%rAttr(sswdn_noise,n)
+ else
+ swndr = avstrm%rAttr(sswdn,n) * 0.50_R8
+ end if
ratio_rvrf = min(0.99_R8,max(0.29548_R8 + 0.00504_R8*swndr &
-1.4957e-05_R8*swndr**2 + 1.4881e-08_R8*swndr**3,0.01_R8))
a2x%rAttr(kswndr,n) = ratio_rvrf*swndr
- swndf = avstrm%rAttr(sswdn,n) * 0.50_R8
+
+ if (sswdn_noise>0) then
+ swndf = avstrm%rAttr(sswdn,n) * 0.50_R8 * avstrm%rAttr(sswdn_noise,n)
+ else
+ swndf = avstrm%rAttr(sswdn,n) * 0.50_R8
+ end if
a2x%rAttr(kswndf,n) = (1._R8 - ratio_rvrf)*swndf
- swvdr = avstrm%rAttr(sswdn,n) * 0.50_R8
+ if (sswdn_noise>0) then
+ if (my_task == master_task .and. n==1) then
+ !write(logunit,*) trim(subname),' swdn before: ', avstrm%rAttr(sswdn,n)
+ !write(logunit,*) trim(subname),' perturbation: ', avstrm%rAttr(sswdn_noise,n)
+ end if
+ swvdr = avstrm%rAttr(sswdn,n) * 0.50_R8 * avstrm%rAttr(sswdn_noise,n)
+ if (my_task == master_task .and. n==1) then
+ !write(logunit,*) trim(subname),' swvdr: ', swvdr
+ end if
+ else
+ swvdr = avstrm%rAttr(sswdn,n) * 0.50_R8
+ end if
ratio_rvrf = min(0.99_R8,max(0.17639_R8 + 0.00380_R8*swvdr &
-9.0039e-06_R8*swvdr**2 + 8.1351e-09_R8*swvdr**3,0.01_R8))
a2x%rAttr(kswvdr,n) = ratio_rvrf*swvdr
- swvdf = avstrm%rAttr(sswdn,n) * 0.50_R8
+ if (sswdn_noise>0) then
+ swvdf = avstrm%rAttr(sswdn,n) * 0.50_R8 * avstrm%rAttr(sswdn_noise,n)
+ else
+ swvdf = avstrm%rAttr(sswdn,n) * 0.50_R8
+ end if
a2x%rAttr(kswvdf,n) = (1._R8 - ratio_rvrf)*swvdf
else
call shr_sys_abort(subName//'ERROR: cannot compute short-wave down')
@@ -1005,8 +1085,21 @@ subroutine datm_comp_run(EClock, x2a, a2x, &
a2x%rAttr(krc,n) = avstrm%rAttr(sprecc,n)
a2x%rAttr(krl,n) = avstrm%rAttr(sprecl,n)
elseif (sprecn > 0) then
- a2x%rAttr(krc,n) = avstrm%rAttr(sprecn,n)*0.1_R8
- a2x%rAttr(krl,n) = avstrm%rAttr(sprecn,n)*0.9_R8
+ ! precipiation perturbation (Yorck)
+ if (sprecn_noise > 0) then
+ if (my_task == master_task .and. n==1) then
+ !write(logunit,*) trim(subname),' precn before: ', avstrm%rAttr(sprecn,n)
+ !write(logunit,*) trim(subname),' perturbation: ', avstrm%rAttr(sprecn_noise,n)
+ end if
+ a2x%rAttr(krc,n) = avstrm%rAttr(sprecn,n)*0.1_R8*avstrm%rAttr(sprecn_noise,n)
+ a2x%rAttr(krl,n) = avstrm%rAttr(sprecn,n)*0.9_R8*avstrm%rAttr(sprecn_noise,n)
+ if (my_task == master_task .and. n==1) then
+ !write(logunit,*) trim(subname),' precn after: ', avstrm%rAttr(sprecn,n)*avstrm%rAttr(sprecn_noise,n)
+ end if
+ else
+ a2x%rAttr(krc,n) = avstrm%rAttr(sprecn,n)*0.1_R8
+ a2x%rAttr(krl,n) = avstrm%rAttr(sprecn,n)*0.9_R8
+ end if
else
call shr_sys_abort(subName//'ERROR: cannot compute rain and snow')
endif
@@ -1040,7 +1133,6 @@ subroutine datm_comp_run(EClock, x2a, a2x, &
a2x%rAttr(ksl,n) = a2x%rAttr(ksl,n)*min(1.e2_r8,avstrm%rAttr(sprecsf,n))
a2x%rAttr(krc,n) = a2x%rAttr(krc,n)*min(1.e2_r8,avstrm%rAttr(sprecsf,n))
a2x%rAttr(krl,n) = a2x%rAttr(krl,n)*min(1.e2_r8,avstrm%rAttr(sprecsf,n))
-
end do
endif
diff --git a/src/drivers/mct/main/cime_comp_mod.F90 b/src/drivers/mct/main/cime_comp_mod.F90
index 6c6df452e55..a8662f9c794 100644
--- a/src/drivers/mct/main/cime_comp_mod.F90
+++ b/src/drivers/mct/main/cime_comp_mod.F90
@@ -40,6 +40,9 @@ module cime_comp_mod
use mct_mod ! mct_ wrappers for mct lib
use perf_mod
use ESMF
+#if defined(USE_OASIS)
+ use mod_oasis
+#endif
!----------------------------------------------------------------------------
! component model interfaces (init, run, final methods)
@@ -580,7 +583,11 @@ module cime_comp_mod
!*******************************************************************************
!===============================================================================
+#ifdef USE_PDAF
subroutine cime_pre_init1(esmf_log_option, pdaf_comm, pdaf_id, pdaf_max)
+#else
+ subroutine cime_pre_init1(esmf_log_option)
+#endif
use shr_pio_mod, only : shr_pio_init1, shr_pio_init2
use seq_comm_mct, only: num_inst_driver
!----------------------------------------------------------
@@ -588,9 +595,11 @@ subroutine cime_pre_init1(esmf_log_option, pdaf_comm, pdaf_id, pdaf_max)
!----------------------------------------------------------
character(CS), intent(out) :: esmf_log_option ! For esmf_logfile_kind
+#ifdef USE_PDAF
integer, optional, intent(in) :: pdaf_comm
integer, optional, intent(in) :: pdaf_id
integer, optional, intent(in) :: pdaf_max
+#endif
integer, dimension(num_inst_total) :: comp_id, comp_comm, comp_comm_iam
logical :: comp_iamin(num_inst_total)
@@ -598,19 +607,48 @@ subroutine cime_pre_init1(esmf_log_option, pdaf_comm, pdaf_id, pdaf_max)
integer :: it
integer :: driver_id
integer :: driver_comm
+#if defined(USE_OASIS)
+ integer :: oas_comp_id
+#endif
+
+#ifndef USE_PDAF
+ call mpi_init(ierr)
+ call shr_mpi_chkerr(ierr,subname//' mpi_init')
+#endif
+#if defined(USE_OASIS)
+ call oasis_init_comp (oas_comp_id, 'eCLM', ierr)
+ if (ierr /= 0) then
+ call oasis_abort(oas_comp_id, 'cime_pre_init1', 'oasis_init_comp error')
+ end if
+ call oasis_get_localcomm(global_comm, ierr)
+ if (ierr /= 0) then
+ call oasis_abort(oas_comp_id, 'cime_pre_init1', 'oasis_get_localcomm error')
+ end if
+#else
+#ifdef USE_PDAF
if (present(pdaf_comm)) then
global_comm = pdaf_comm
else
call mpi_comm_dup(MPI_COMM_WORLD, global_comm, ierr)
call shr_mpi_chkerr(ierr,subname//' mpi_comm_dup')
end if
+#else
+ call mpi_comm_dup(MPI_COMM_WORLD, global_comm, ierr)
+ call shr_mpi_chkerr(ierr,subname//' mpi_comm_dup')
+#endif
+#endif
+
comp_comm = MPI_COMM_NULL
time_brun = mpi_wtime()
!--- Initialize multiple driver instances, if requested ---
+#ifdef USE_PDAF
call cime_cpl_init(global_comm, driver_comm, num_inst_driver, driver_id, &
pdaf_id, pdaf_max)
+#else
+ call cime_cpl_init(global_comm, driver_comm, num_inst_driver, driver_id)
+#endif
call shr_pio_init1(num_inst_total,NLFileName, driver_comm)
!
@@ -622,9 +660,11 @@ subroutine cime_pre_init1(esmf_log_option, pdaf_comm, pdaf_id, pdaf_max)
if (num_inst_driver > 1) then
call seq_comm_init(global_comm, driver_comm, NLFileName, drv_comm_ID=driver_id)
write(cpl_inst_tag,'("_",i4.4)') driver_id
+#ifdef USE_PDAF
else if (present(pdaf_id) .and. present(pdaf_max)) then
call seq_comm_init(global_comm, driver_comm, NLFileName, &
pdaf_id=pdaf_id, pdaf_max=pdaf_max)
+#endif
else
call seq_comm_init(global_comm, driver_comm, NLFileName)
cpl_inst_tag = ''
@@ -4168,6 +4208,10 @@ subroutine cime_final()
mpicom=mpicom_GLOID)
endif
+#if defined(USE_OASIS)
+ call oasis_terminate (ierr)
+ call shr_mpi_chkerr(ierr,subname//' oasis_terminate')
+#endif
call t_finalizef()
end subroutine cime_final
@@ -4238,8 +4282,12 @@ subroutine cime_comp_barriers(mpicom, timer)
endif
end subroutine cime_comp_barriers
+#ifdef USE_PDAF
subroutine cime_cpl_init(comm_in, comm_out, num_inst_driver, id, &
pdaf_id, pdaf_max)
+#else
+ subroutine cime_cpl_init(comm_in, comm_out, num_inst_driver, id)
+#endif
!-----------------------------------------------------------------------
!
! Initialize multiple coupler instances, if requested
@@ -4252,8 +4300,9 @@ subroutine cime_cpl_init(comm_in, comm_out, num_inst_driver, id, &
integer , intent(out) :: comm_out
integer , intent(out) :: num_inst_driver
integer , intent(out) :: id ! instance ID, starts from 1
+#ifdef USE_PDAF
integer , intent(in), optional :: pdaf_id, pdaf_max
-
+#endif
!
! Local variables
!
@@ -4295,10 +4344,14 @@ subroutine cime_cpl_init(comm_in, comm_out, num_inst_driver, id, &
' : Total PE number must be a multiple of coupler instance number')
end if
+#ifdef USE_PDAF
if (pdaf_max > 1) then
call mpi_comm_split(comm_in, pdaf_id, 0, comm_out, ierr)
call shr_mpi_chkerr(ierr,subname//' mpi_comm_split')
else if (num_inst_driver == 1) then
+#else
+ if (num_inst_driver == 1) then
+#endif
call mpi_comm_dup(comm_in, comm_out, ierr)
call shr_mpi_chkerr(ierr,subname//' mpi_comm_dup')
else
diff --git a/src/drivers/mct/shr/seq_comm_mct.F90 b/src/drivers/mct/shr/seq_comm_mct.F90
index 2f0d346ffa3..5da00275cfb 100644
--- a/src/drivers/mct/shr/seq_comm_mct.F90
+++ b/src/drivers/mct/shr/seq_comm_mct.F90
@@ -198,8 +198,12 @@ integer function seq_comm_get_ncomps()
seq_comm_get_ncomps = ncomps
end function seq_comm_get_ncomps
- subroutine seq_comm_init(global_comm_in,driver_comm_in,nmlfile,drv_comm_id,&
+#ifdef USE_PDAF
+ subroutine seq_comm_init(global_comm_in, driver_comm_in, nmlfile, drv_comm_id,&
pdaf_id, pdaf_max)
+#else
+ subroutine seq_comm_init(global_comm_in, driver_comm_in, nmlfile, drv_comm_id)
+#endif
!----------------------------------------------------------
!
! Arguments
@@ -208,8 +212,10 @@ subroutine seq_comm_init(global_comm_in,driver_comm_in,nmlfile,drv_comm_id,&
integer, intent(in) :: driver_comm_in
character(len=*), intent(IN) :: nmlfile
integer, intent(in), optional :: drv_comm_id
+#ifdef USE_PDAF
integer, intent(in), optional :: pdaf_id
integer, intent(in), optional :: pdaf_max
+#endif
!
! Local variables
!
@@ -407,6 +413,7 @@ subroutine seq_comm_init(global_comm_in,driver_comm_in,nmlfile,drv_comm_id,&
call mpi_bcast(pelist, size(pelist), MPI_INTEGER, 0, DRIVER_COMM, ierr)
call seq_comm_setcomm(CPLID,pelist,nthreads=cpl_nthreads,iname='CPL')
+#ifdef USE_PDAF
call comp_comm_init(driver_comm, atm_rootpe, atm_nthreads, atm_layout, atm_ntasks, atm_pestride, num_inst_atm, &
CPLID, ATMID, CPLATMID, ALLATMID, CPLALLATMID, 'ATM', count, drv_comm_id, pdaf_id, pdaf_max)
call comp_comm_init(driver_comm, lnd_rootpe, lnd_nthreads, lnd_layout, lnd_ntasks, lnd_pestride, num_inst_lnd, &
@@ -423,6 +430,24 @@ subroutine seq_comm_init(global_comm_in,driver_comm_in,nmlfile,drv_comm_id,&
CPLID, WAVID, CPLWAVID, ALLWAVID, CPLALLWAVID, 'WAV', count, drv_comm_id, pdaf_id, pdaf_max)
call comp_comm_init(driver_comm, esp_rootpe, esp_nthreads, esp_layout, esp_ntasks, esp_pestride, num_inst_esp, &
CPLID, ESPID, CPLESPID, ALLESPID, CPLALLESPID, 'ESP', count, drv_comm_id, pdaf_id, pdaf_max)
+#else
+ call comp_comm_init(driver_comm, atm_rootpe, atm_nthreads, atm_layout, atm_ntasks, atm_pestride, num_inst_atm, &
+ CPLID, ATMID, CPLATMID, ALLATMID, CPLALLATMID, 'ATM', count, drv_comm_id)
+ call comp_comm_init(driver_comm, lnd_rootpe, lnd_nthreads, lnd_layout, lnd_ntasks, lnd_pestride, num_inst_lnd, &
+ CPLID, LNDID, CPLLNDID, ALLLNDID, CPLALLLNDID, 'LND', count, drv_comm_id)
+ call comp_comm_init(driver_comm, ice_rootpe, ice_nthreads, ice_layout, ice_ntasks, ice_pestride, num_inst_ice, &
+ CPLID, ICEID, CPLICEID, ALLICEID, CPLALLICEID, 'ICE', count, drv_comm_id)
+ call comp_comm_init(driver_comm, ocn_rootpe, ocn_nthreads, ocn_layout, ocn_ntasks, ocn_pestride, num_inst_ocn, &
+ CPLID, OCNID, CPLOCNID, ALLOCNID, CPLALLOCNID, 'OCN', count, drv_comm_id)
+ call comp_comm_init(driver_comm, rof_rootpe, rof_nthreads, rof_layout, rof_ntasks, rof_pestride, num_inst_rof, &
+ CPLID, ROFID, CPLROFID, ALLROFID, CPLALLROFID, 'ROF', count, drv_comm_id)
+ call comp_comm_init(driver_comm, glc_rootpe, glc_nthreads, glc_layout, glc_ntasks, glc_pestride, num_inst_glc, &
+ CPLID, GLCID, CPLGLCID, ALLGLCID, CPLALLGLCID, 'GLC', count, drv_comm_id)
+ call comp_comm_init(driver_comm, wav_rootpe, wav_nthreads, wav_layout, wav_ntasks, wav_pestride, num_inst_wav, &
+ CPLID, WAVID, CPLWAVID, ALLWAVID, CPLALLWAVID, 'WAV', count, drv_comm_id)
+ call comp_comm_init(driver_comm, esp_rootpe, esp_nthreads, esp_layout, esp_ntasks, esp_pestride, num_inst_esp, &
+ CPLID, ESPID, CPLESPID, ALLESPID, CPLALLESPID, 'ESP', count, drv_comm_id)
+#endif
if (count /= ncomps) then
write(logunit,*) trim(subname),' ERROR in ID count ',count,ncomps
@@ -496,10 +521,16 @@ subroutine seq_comm_init(global_comm_in,driver_comm_in,nmlfile,drv_comm_id,&
end subroutine seq_comm_init
+#ifdef USE_PDAF
subroutine comp_comm_init(driver_comm, comp_rootpe, comp_nthreads, comp_layout, &
comp_ntasks, comp_pestride, num_inst_comp, &
CPLID, COMPID, CPLCOMPID, ALLCOMPID, CPLALLCOMPID, name, count, drv_comm_id, &
pdaf_id, pdaf_max)
+#else
+ subroutine comp_comm_init(driver_comm, comp_rootpe, comp_nthreads, comp_layout, &
+ comp_ntasks, comp_pestride, num_inst_comp, &
+ CPLID, COMPID, CPLCOMPID, ALLCOMPID, CPLALLCOMPID, name, count, drv_comm_id)
+#endif
integer, intent(in) :: driver_comm
integer, intent(in) :: comp_rootpe
integer, intent(in) :: comp_nthreads
@@ -514,8 +545,10 @@ subroutine comp_comm_init(driver_comm, comp_rootpe, comp_nthreads, comp_layout,
integer, intent(out) :: CPLALLCOMPID
integer, intent(inout) :: count
integer, intent(in), optional :: drv_comm_id
+#ifdef USE_PDAF
integer, intent(in), optional :: pdaf_id
integer, intent(in), optional :: pdaf_max
+#endif
character(len=*), intent(in) :: name
character(len=*), parameter :: subname = "comp_comm_init"
@@ -578,11 +611,15 @@ subroutine comp_comm_init(driver_comm, comp_rootpe, comp_nthreads, comp_layout,
pelist(2,1) = cmax(n)
pelist(3,1) = cstr(n)
endif
-
call mpi_bcast(pelist, size(pelist), MPI_INTEGER, 0, DRIVER_COMM, ierr)
+
+#ifdef USE_PDAF
if (present(pdaf_id) .and. present(pdaf_max)) then
call seq_comm_setcomm(COMPID(n),pelist,nthreads=comp_nthreads,iname=name,inst=pdaf_id,tinst=pdaf_max)
else if (present(drv_comm_id)) then
+#else
+ if (present(drv_comm_id)) then
+#endif
call seq_comm_setcomm(COMPID(n),pelist,nthreads=comp_nthreads,iname=name,inst=drv_comm_id)
else
call seq_comm_setcomm(COMPID(n),pelist,nthreads=comp_nthreads,iname=name,inst=n,tinst=num_inst_comp)
diff --git a/src/share/streams/shr_dmodel_mod.F90 b/src/share/streams/shr_dmodel_mod.F90
index 011cb8ccd01..394d0d52af0 100644
--- a/src/share/streams/shr_dmodel_mod.F90
+++ b/src/share/streams/shr_dmodel_mod.F90
@@ -764,6 +764,12 @@ subroutine shr_dmodel_readstrm(stream, pio_subsystem, pio_iotype, pio_iodesc, gs
logical :: fileopen
character(CL) :: currfile
+ ! Yorck:
+ character(CL) :: mfldName
+ integer :: position ! integer to look if string contains noise
+ integer :: caseId ! which ensemble member?
+ integer :: dt ! time sampling of forcings
+
integer(in) :: ndims
integer(in),pointer :: dimid(:)
type(file_desc_t) :: pioid
@@ -910,9 +916,37 @@ subroutine shr_dmodel_readstrm(stream, pio_subsystem, pio_iotype, pio_iodesc, gs
if (my_task == master_task) then
call shr_stream_getFileFieldName(stream,k,sfldName)
endif
+ ! Yorck adjustment
+ if (my_task == master_task) then
+ call shr_stream_getModelFieldName(stream,k,mfldName)
+ end if
+ call shr_mpi_bcast(mfldName,mpicom,'mfldName')
call shr_mpi_bcast(sfldName,mpicom,'sfldName')
rcode = pio_inq_varid(pioid,trim(sfldName),varid)
- frame = nt
+ !------------------------------------------------------------------------------------
+ ! Yorck adjustments :
+ ! if noise paramter, the frame has to be adjusted depending on the ensemble member
+ ! it is nt + caseId*24/dt, so that for case 0 (first ensemble member) it takes for
+ ! the first time step the first frame, for ens mem 1 the 8th (3 hourly data) and
+ ! so on. With this we ensure different forcings in time and time correlations.
+ ! For now, the time is hard coded. Not sure how to work around this...
+ ! idea: plug time and ensemble member into stream
+
+ if (INDEX(mfldName, "noise") == 0) then
+ frame = nt
+ else
+ if (my_task == master_task) then
+ call shr_stream_get_dt_and_caseId(stream,dt,caseId)
+ end if
+
+ call shr_mpi_bcast(caseId,mpicom)
+ call shr_mpi_bcast(dt,mpicom)
+
+ if (dt == 0) then
+ call shr_sys_abort(subname//"dt is zero, time resolution of forcing data has to be provided in noise streams")
+ end if
+ frame = nt+caseId*24/dt
+ end if
call pio_setframe(pioid,varid,frame)
call pio_read_darray(pioid,varid,pio_iodesc,av%rattr(k,:),rcode)
enddo
diff --git a/src/share/streams/shr_strdata_mod.F90 b/src/share/streams/shr_strdata_mod.F90
index 3b104d58e7d..60f2abc9d7b 100644
--- a/src/share/streams/shr_strdata_mod.F90
+++ b/src/share/streams/shr_strdata_mod.F90
@@ -1385,17 +1385,17 @@ end subroutine shr_strdata_pioinit_oldway
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_strdata_create_newway(SDAT, name, mpicom, compid, gsmap, ggrid, nxg, nyg, &
- !--- streams stuff required ---
- yearFirst, yearLast, yearAlign, offset, &
- domFilePath, domFileName, &
- domTvarName, domXvarName, domYvarName, domAreaName, domMaskName, &
- filePath, filename, fldListFile, fldListModel, &
- !--- strdata optional ---
- nzg, domZvarName, &
- taxMode, dtlimit, tintalgo, readmode, &
- fillalgo, fillmask, fillread, fillwrite, &
- mapalgo, mapmask, mapread, mapwrite, &
- calendar)
+ !--- streams stuff required ---
+ yearFirst, yearLast, yearAlign, offset, &
+ domFilePath, domFileName, &
+ domTvarName, domXvarName, domYvarName, domAreaName, domMaskName, &
+ filePath, filename, fldListFile, fldListModel, &
+ !--- strdata optional ---
+ caseId, dt, numEns, nzg, domZvarName, &
+ taxMode, dtlimit, tintalgo, readmode, &
+ fillalgo, fillmask, fillread, fillwrite, &
+ mapalgo, mapmask, mapread, mapwrite, &
+ calendar)
implicit none
@@ -1441,6 +1441,11 @@ subroutine shr_strdata_create_newway(SDAT, name, mpicom, compid, gsmap, ggrid, n
character(*),optional ,intent(in) :: readmode ! file read mode
character(*),optional, intent(in) :: calendar
+ ! Yorck
+ integer(IN) ,optional ,intent(in) :: caseId
+ integer(IN) ,optional ,intent(in) :: dt
+ integer(IN) ,optional ,intent(in) :: numEns
+
!EOP
! --- local ---
@@ -1503,20 +1508,37 @@ subroutine shr_strdata_create_newway(SDAT, name, mpicom, compid, gsmap, ggrid, n
endif
!---- Backwards compatibility requires Z be optional
+ ! Yorck adjustments: optionality of caseId, dt and numEns
if (present(domZvarName)) then
- call shr_stream_set(SDAT%stream(1),yearFirst,yearLast,yearAlign,offset,taxMode, &
- fldListFile,fldListModel,domFilePath,domFileName, &
- domTvarName,domXvarName,domYvarName,domZvarName, &
- domAreaName,domMaskName, &
- filePath,filename,trim(name))
- else
- call shr_stream_set(SDAT%stream(1),yearFirst,yearLast,yearAlign,offset,taxMode, &
- fldListFile,fldListModel,domFilePath,domFileName, &
- domTvarName,domXvarName,domYvarName,'undefined', &
- domAreaName,domMaskName, &
- filePath,filename,trim(name))
- endif
+ if (present(caseId)) then
+ call shr_stream_set(SDAT%stream(1),yearFirst,yearLast,yearAlign,offset,caseId,dt,numEns, &
+ taxMode,fldListFile,fldListModel,domFilePath,domFileName, &
+ domTvarName,domXvarName,domYvarName,domZvarName, &
+ domAreaName,domMaskName, &
+ filePath,filename,trim(name))
+ else
+ call shr_stream_set(SDAT%stream(1),yearFirst,yearLast,yearAlign,offset,0,0,0, &
+ taxMode,fldListFile,fldListModel,domFilePath,domFileName, &
+ domTvarName,domXvarName,domYvarName,domZvarName, &
+ domAreaName,domMaskName, &
+ filePath,filename,trim(name))
+ end if
+ else
+ if (present(caseId)) then
+ call shr_stream_set(SDAT%stream(1),yearFirst,yearLast,yearAlign,offset,caseId,dt,numEns, &
+ taxMode,fldListFile,fldListModel,domFilePath,domFileName, &
+ domTvarName,domXvarName,domYvarName,'undefined', &
+ domAreaName,domMaskName, &
+ filePath,filename,trim(name))
+ else
+ call shr_stream_set(SDAT%stream(1),yearFirst,yearLast,yearAlign,offset,0,0,0, &
+ taxMode,fldListFile,fldListModel,domFilePath,domFileName, &
+ domTvarName,domXvarName,domYvarName,'undefined', &
+ domAreaName,domMaskName, &
+ filePath,filename,trim(name))
+ end if
+ endif
if (present(nzg)) then
call shr_strdata_init(SDAT, mpicom, compid, &
@@ -1536,7 +1558,7 @@ subroutine shr_strdata_create_oldway(SDAT, name, mpicom, compid, gsmap, ggrid, n
filePath, filename, fldListFile, fldListModel, &
pio_subsystem, pio_iotype, &
!--- strdata optional ---
- nzg, domZvarName, &
+ caseId, dt, numEns, nzg, domZvarName, &
taxMode, dtlimit, tintalgo, readmode, &
fillalgo, fillmask, fillread, fillwrite, &
mapalgo, mapmask, mapread, mapwrite, &
@@ -1589,12 +1611,17 @@ subroutine shr_strdata_create_oldway(SDAT, name, mpicom, compid, gsmap, ggrid, n
character(*),optional ,intent(in) :: readmode ! file read mode
character(*),optional, intent(in) :: calendar
+ ! Yorck
+ integer(IN) ,optional ,intent(in) :: caseId
+ integer(IN) ,optional ,intent(in) :: dt
+ integer(IN) ,optional ,intent(in) :: numEns
+
!EOP
! pio variables are already in SDAT no need to copy them
call shr_strdata_create_newway(SDAT, name, mpicom, compid, gsmap, ggrid, nxg, nyg,&
yearFirst, yearLast, yearAlign, offset, domFilePath, domFilename, domTvarName, &
domXvarName, domYvarName, domAreaName, domMaskName, &
- filePath, filename, fldListFile, fldListModel, nzg=nzg, domZvarName=domZvarName,&
+ filePath, filename, fldListFile, fldListModel, caseId=caseId, dt=dt, numEns=numEns, nzg=nzg, domZvarName=domZvarName,&
taxMode=taxMode,dtlimit=dtlimit,tintalgo=tintalgo,readmode=readmode,&
fillalgo=fillalgo,fillmask=fillmask,fillread=fillread,fillwrite=fillwrite,mapalgo=mapalgo,&
mapmask=mapmask,mapread=mapread,mapwrite=mapwrite,calendar=calendar)
diff --git a/src/share/streams/shr_stream_mod.F90 b/src/share/streams/shr_stream_mod.F90
index 1e8d5c0bedf..f99b016f527 100644
--- a/src/share/streams/shr_stream_mod.F90
+++ b/src/share/streams/shr_stream_mod.F90
@@ -77,6 +77,9 @@ module shr_stream_mod
public :: shr_stream_setAbort ! set internal shr_stream abort flag
public :: shr_stream_getDebug ! get internal shr_stream debug level
public :: shr_stream_isInit ! check if stream is initialized
+
+ ! Yorck
+ public :: shr_stream_get_dt_and_caseId ! return dt and caseId
! public :: shr_stream_bcast ! broadcast a stream (untested)
! !PUBLIC DATA MEMBERS:
@@ -151,6 +154,14 @@ module shr_stream_mod
character(SHR_KIND_CS) :: domAreaName ! domain file: area var name
character(SHR_KIND_CS) :: domMaskName ! domain file: mask var name
+ ! Yorck: in the stream, also the time resolution of the forcing data and the caseId
+ ! is saved. This is done to get the right frame of the noise file later.
+ ! Also: the number of ensemble members the noise data was produced for has
+ ! to be included to make sure that the time information is right
+ integer(SHR_KIND_IN) :: caseId
+ integer(SHR_KIND_IN) :: dt
+ integer(SHR_KIND_IN) :: numEns
+
character(SHR_KIND_CS) :: tInterpAlgo ! Algorithm to use for time interpolation
character(SHR_KIND_CL) :: calendar ! stream calendar
end type shr_stream_streamType
@@ -642,6 +653,76 @@ subroutine shr_stream_init(strm,infoFile,yearFirst,yearLast,yearAlign,taxMode,rc
close(nUnit)
+ !-----------------------------------------------------------------------------
+ ! Yorck adaptations: read in time resolution of forcings and caseId, numEns
+ !-----------------------------------------------------------------------------
+
+ !--- find start tag ---
+ open(nUnit,file=infoFile,STATUS='OLD',FORM='FORMATTED',ACTION='READ')
+ startTag = ""
+ endTag = ""
+ call shr_stream_readUpToTag(nUnit,startTag,rc=rCode2)
+ if (rCode2 /= 0)then
+ rCode = rCode2
+ goto 999
+ end if
+ startTag = ""
+ endTag = ""
+ call shr_stream_readUpToTag(nUnit,startTag,optionalTag=.true.,rc=rCode2)
+ if (rCode2 == 0) then
+ !--- read data ---
+ read(nUnit,*,END=999) int
+ strm%caseId = int
+ else
+ strm%caseId = 0
+ end if
+
+ close(nUnit)
+
+ !--- find start tag ---
+ open(nUnit,file=infoFile,STATUS='OLD',FORM='FORMATTED',ACTION='READ')
+ startTag = ""
+ endTag = ""
+ call shr_stream_readUpToTag(nUnit,startTag,rc=rCode2)
+ if (rCode2 /= 0)then
+ rCode = rCode2
+ goto 999
+ end if
+ startTag = ""
+ endTag = ""
+ call shr_stream_readUpToTag(nUnit,startTag,optionalTag=.true.,rc=rCode2)
+ if (rCode2 == 0) then
+ !--- read data ---
+ read(nUnit,*,END=999) int
+ strm%dt = int
+ else
+ strm%dt = 0
+ end if
+
+ close(nUnit)
+
+ !--- find start tag ---
+ open(nUnit,file=infoFile,STATUS='OLD',FORM='FORMATTED',ACTION='READ')
+ startTag = ""
+ endTag = ""
+ call shr_stream_readUpToTag(nUnit,startTag,rc=rCode2)
+ if (rCode2 /= 0)then
+ rCode = rCode2
+ goto 999
+ end if
+ startTag = ""
+ endTag = ""
+ call shr_stream_readUpToTag(nUnit,startTag,optionalTag=.true.,rc=rCode2)
+ if (rCode2 == 0) then
+ !--- read data ---
+ read(nUnit,*,END=999) int
+ strm%numEns = int
+ else
+ strm%numEns = 0
+ end if
+
+ close(nUnit)
+
!-----------------------------------------------------------------------------
! get initial calendar value
!-----------------------------------------------------------------------------
@@ -681,11 +762,11 @@ end subroutine shr_stream_init
!
! !INTERFACE: ------------------------------------------------------------------
- subroutine shr_stream_set(strm,yearFirst,yearLast,yearAlign,offset,taxMode, &
- fldListFile,fldListModel,domFilePath,domFileName, &
- domTvarName,domXvarName,domYvarName,domZvarName, &
- domAreaName,domMaskName, &
- filePath,filename,dataSource,rc)
+ subroutine shr_stream_set(strm,yearFirst,yearLast,yearAlign,offset, caseId, dt, &
+ numEns, taxMode, fldListFile,fldListModel,domFilePath,domFileName, &
+ domTvarName,domXvarName,domYvarName,domZvarName, &
+ domAreaName,domMaskName, &
+ filePath,filename,dataSource,rc)
! !INPUT/OUTPUT PARAMETERS:
@@ -710,6 +791,11 @@ subroutine shr_stream_set(strm,yearFirst,yearLast,yearAlign,offset,taxMode, &
character(*) ,optional,intent(in) :: dataSource ! comment line
integer (SHR_KIND_IN),optional,intent(out) :: rc ! return code
+ ! Yorck
+ integer (SHR_KIND_IN),optional,intent(in) :: caseId ! ensemble member case ID
+ integer (SHR_KIND_IN),optional,intent(in) :: dt ! time resolution of forcing data
+ integer (SHR_KIND_IN),optional,intent(in) :: numEns ! number of ensemble members set to produce noise data
+
!EOP
!----- local -----
@@ -791,6 +877,19 @@ subroutine shr_stream_set(strm,yearFirst,yearLast,yearAlign,offset,taxMode, &
call fileVec%move_out(strm%file)
endif
+ ! Yorck
+ if (present(caseId)) then
+ strm%caseId = caseId
+ end if
+
+ if (present(dt)) then
+ strm%dt = dt
+ end if
+
+ if (present(numEns)) then
+ strm%numEns = numEns
+ end if
+
!-----------------------------------------------------------------------------
! get initial calendar value
!-----------------------------------------------------------------------------
@@ -871,6 +970,10 @@ subroutine shr_stream_default(strm,rc)
strm%domAreaName = ' '
strm%domMaskName = ' '
+ strm%caseId = 0
+ strm%dt = 0
+ strm%numEns = 0
+
strm%calendar = shr_cal_noleap
if ( present(rc) ) rc = 0
@@ -1516,6 +1619,9 @@ subroutine shr_stream_readTCoord(strm,k,rc)
integer(SHR_KIND_IN) :: ndate ! calendar date of time value
real(SHR_KIND_R8) :: nsec ! elapsed secs on calendar date
real(SHR_KIND_R8),allocatable :: tvar(:)
+ character(SHR_KIND_CL) :: mfldName !Yorck
+ integer(SHR_KIND_IN) :: nt_ ! Yorck
+ real(SHR_KIND_R8),allocatable :: tvar_lok(:) ! Yorck
!----- formats -----
character(*),parameter :: subname = '(shr_stream_readTCoord) '
character(*),parameter :: F01 = "('(shr_stream_readTCoord) ',a,2i7)"
@@ -1539,6 +1645,27 @@ subroutine shr_stream_readTCoord(strm,k,rc)
if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_inquire_dimension')
deallocate(dids)
+ ! as the time information in the noise files is wrong, the last timesteps which only resemble the number of ensemble
+ ! members and the time resolution of the forcing data have to be deleted. It is not necessarily the number
+ ! of ensemble members the experiment is conducted with but the number of ensemble members the noise data was produced with.
+
+ ! After nt is adjusted, also the read in number of timesteps has to be adjusted. I hope this works...
+ ! this could lead to errors, the first setted variable inside the noise stream has to be a noise variable
+ ! but as noise is the only thing setted inside these streams, this should still work if no one changes this...
+
+ call shr_stream_getModelFieldName(strm,1,mfldName)
+ if (INDEX(mfldName, "noise") .ne. 0) then
+ nt_ = nt
+ if (strm%numEns == 0) then
+ call shr_sys_abort(subname//' ERROR: number of ensemble members that were used for producing noise data has to be set in stream')
+ end if
+ if (strm%dt == 0) then
+ call shr_sys_abort(subname//' ERROR: time resolution of forcing data has to be set in stream')
+ end if
+ ! (numEns-1)*(24/dt) is the number of extra time steps introduced in the noise preperation
+ nt = nt - (strm%numEns-1)*(24/strm%dt)
+ end if
+
allocate(strm%file(k)%date(nt),strm%file(k)%secs(nt))
strm%file(k)%nt = nt
@@ -1563,10 +1690,27 @@ subroutine shr_stream_readTCoord(strm,k,rc)
strm%calendar = trim(shr_cal_calendarName(trim(calendar)))
allocate(tvar(nt))
- rcode = nf90_get_var(fid,vid,tvar)
- if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_get_var')
- rCode = nf90_close(fid)
- if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_close')
+
+ ! Yorck
+ if (INDEX(mfldName, "noise") .ne. 0) then
+ allocate(tvar_lok(nt_))
+ rcode = nf90_get_var(fid,vid,tvar_lok)
+ if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_get_var')
+ rCode = nf90_close(fid)
+ if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_close')
+
+ do n = 1, nt
+ tvar(n) = tvar_lok(n)
+ end do
+ deallocate(tvar_lok)
+
+ else
+ rcode = nf90_get_var(fid,vid,tvar)
+ if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_get_var')
+ rCode = nf90_close(fid)
+ if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_close')
+ end if
+ ! End Yorck
do n = 1,nt
call shr_cal_advDate(tvar(n),bunits,bdate,bsec,ndate,nsec,calendar)
strm%file(k)%date(n) = ndate
@@ -2582,6 +2726,10 @@ subroutine shr_stream_restWrite(strm,fileName,caseName,caseDesc,nstrms,rc)
write(nUnit) strm(k)%domAreaName ! domain file: area var name
write(nUnit) strm(k)%domMaskName ! domain file: mask var name
+ write(nUnit) strm(k)%dt ! Yorck
+ write(nUnit) strm(k)%caseId ! Yorck
+ write(nUnit) strm(k)%numEns ! Yorck
+
end do
close(nUnit)
@@ -2885,6 +3033,10 @@ subroutine shr_stream_dataDump(strm)
write(s_logunit,F00) "domAreaName = ", trim(strm%domAreaName)
write(s_logunit,F00) "domMaskName = ", trim(strm%domMaskName)
+ write(s_logunit,F01) "dt = ", strm%dt ! Yorck
+ write(s_logunit,F01) "caseId = ", strm%caseId ! Yorck
+ write(s_logunit,F01) "numEns = ", strm%numEns !Yorck
+
end subroutine shr_stream_dataDump
!===============================================================================
@@ -3232,6 +3384,10 @@ subroutine shr_stream_bcast(stream,comm,rc)
call shr_mpi_bcast(stream%domZvarName ,comm,subName)
call shr_mpi_bcast(stream%domMaskName ,comm,subName)
call shr_mpi_bcast(stream%calendar ,comm,subName)
+
+ !Yorck
+ call shr_mpi_bcast(stream%dt ,comm,subName)
+ call shr_mpi_bcast(stream%caseId ,comm,subName)
if (pid /= 0) allocate(stream%file(stream%nFiles))
@@ -3247,6 +3403,29 @@ subroutine shr_stream_bcast(stream,comm,rc)
end subroutine shr_stream_bcast
+
+ subroutine shr_stream_get_dt_and_caseId(stream,dt,caseId)
+
+ implicit none
+
+ ! !INPUT/OUTPUT PARAMETERS:
+
+ type(shr_stream_streamType) ,intent(in) :: stream ! stream in question
+ integer(SHR_KIND_IN) ,intent(out) :: dt !
+ integer(SHR_KIND_IN) ,intent(out) :: caseId !
+
+
+
+ !-------------------------------------------------------------------------------
+ ! Notes:
+ !-------------------------------------------------------------------------------
+
+ dt = stream%dt
+ caseId = stream%caseId
+
+
+ end subroutine shr_stream_get_dt_and_caseId
+
!===============================================================================
end module shr_stream_mod
!===============================================================================