Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Changes from perturbationroutineclm5 #1

Open
wants to merge 5 commits into
base: clm5_0-tsmp-pdaf
Choose a base branch
from
Open

Changes from perturbationroutineclm5 #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
e4d681e
changes from perturbationroutineclm5
s7yoewer May 17, 2024
8e140f9
Merge branch 'clm5_0-tsmp-pdaf' into clm5_0-tsmp-pdaf-perturbationrou…
jjokella May 23, 2024
a05d2e7
TSMP-PDAF: import changes from eCLM
jjokella Jun 19, 2024
4a31d3d
Merge branch 'clm5_0-tsmp-pdaf' into clm5_0-tsmp-pdaf-perturbationrou…
jjokella Jun 19, 2024
c98f167
Merge branch 'clm5_0-tsmp-pdaf' into clm5_0-tsmp-pdaf-perturbationrou…
jjokella Jun 26, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
122 changes: 107 additions & 15 deletions src/components/data_comps/datm/datm_comp_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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 ", &
Expand All @@ -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) = &
Expand All @@ -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
Expand Down Expand Up @@ -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')
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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 ---
Expand Down Expand Up @@ -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')
Expand All @@ -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
Expand Down Expand Up @@ -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

Expand Down
55 changes: 54 additions & 1 deletion src/drivers/mct/main/cime_comp_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -580,37 +583,72 @@ 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
!----------------------------------------------------------
!| Initialize MCT and MPI communicators and IO
!----------------------------------------------------------

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)
character(len=seq_comm_namelen) :: comp_name(num_inst_total)
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)
!
Expand All @@ -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 = ''
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
!
Expand Down Expand Up @@ -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
Expand Down
Loading