From fa291b5fa3d556a2a6bd6337a4ce94476b984335 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 9 Jul 2024 11:29:58 +0200 Subject: [PATCH 01/21] Add in-code perturbation routine for eCLM-PDAF from --- README.md | 4 + src/csm_share/streams/shr_dmodel_mod.F90 | 36 +++- src/csm_share/streams/shr_strdata_mod.F90 | 77 ++++++--- src/csm_share/streams/shr_stream_mod.F90 | 197 +++++++++++++++++++++- src/datm/datm_comp_mod.F90 | 116 +++++++++++-- 5 files changed, 383 insertions(+), 47 deletions(-) diff --git a/README.md b/README.md index 48387d917c..dbfe5ac2b8 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,7 @@ +`dev-perturbation`-Branch: Generated for usage of eCLM-PDAF with +Perturbation routine from +https://gitlab.jsc.fz-juelich.de/detect/cluster-c/c01/perturbationroutineclm5 + # eCLM [![CI](https://github.com/HPSCTerrSys/eCLM/actions/workflows/CI.yml/badge.svg)](https://github.com/HPSCTerrSys/eCLM/actions/workflows/CI.yml) diff --git a/src/csm_share/streams/shr_dmodel_mod.F90 b/src/csm_share/streams/shr_dmodel_mod.F90 index 011cb8ccd0..394d0d52af 100644 --- a/src/csm_share/streams/shr_dmodel_mod.F90 +++ b/src/csm_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/csm_share/streams/shr_strdata_mod.F90 b/src/csm_share/streams/shr_strdata_mod.F90 index 3b104d58e7..60f2abc9d7 100644 --- a/src/csm_share/streams/shr_strdata_mod.F90 +++ b/src/csm_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/csm_share/streams/shr_stream_mod.F90 b/src/csm_share/streams/shr_stream_mod.F90 index 1e8d5c0bed..f99b016f52 100644 --- a/src/csm_share/streams/shr_stream_mod.F90 +++ b/src/csm_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 !=============================================================================== diff --git a/src/datm/datm_comp_mod.F90 b/src/datm/datm_comp_mod.F90 index 0b7744d1da..af76859881 100644 --- a/src/datm/datm_comp_mod.F90 +++ b/src/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') @@ -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 From e666707e730b24c1268462616e8942d023cd4521 Mon Sep 17 00:00:00 2001 From: Yorck Ewerdwalbesloh Date: Fri, 25 Jul 2025 10:57:13 +0200 Subject: [PATCH 02/21] GRACE-DA: additional arrays, updated perturbation routines --- src/clm5/biogeophys/SoilHydrologyType.F90 | 8 + .../biogeophys/SoilStateInitTimeConstMod.F90 | 110 ++++++++++- src/clm5/biogeophys/WaterStateType.F90 | 92 +++++++++ src/clm5/cpl/lnd_comp_mct.F90 | 2 +- src/clm5/main/atm2lndType.F90 | 10 + src/clm5/main/clm_driver.F90 | 2 +- src/clm5/main/clm_varcon.F90 | 8 + src/clm5/main/lnd2atmMod.F90 | 186 +++++++++++++++++- src/csm_share/streams/shr_dmodel_mod.F90 | 1 - src/csm_share/streams/shr_stream_mod.F90 | 6 +- src/datm/datm_comp_mod.F90 | 95 +++------ 11 files changed, 442 insertions(+), 78 deletions(-) diff --git a/src/clm5/biogeophys/SoilHydrologyType.F90 b/src/clm5/biogeophys/SoilHydrologyType.F90 index 33aee7efda..af6de9165e 100644 --- a/src/clm5/biogeophys/SoilHydrologyType.F90 +++ b/src/clm5/biogeophys/SoilHydrologyType.F90 @@ -53,6 +53,10 @@ Module SoilHydrologyType real(r8), pointer :: i_0_col (:) ! col VIC average saturation in top soil layers real(r8), pointer :: ice_col (:,:) ! col VIC soil ice (kg/m2) for VIC soil layers + ! Yorck + real(r8), pointer :: wa_col_inc (:) ! increment col water in the unconfined aquifer (mm) + real(r8), pointer :: wa_col_mean (:) ! mean col water in the unconfined aquifer (mm) + contains ! Public routines @@ -117,6 +121,10 @@ subroutine InitAllocate(this, bounds) allocate(this%zwts_col (begc:endc)) ; this%zwts_col (:) = nan allocate(this%wa_col (begc:endc)) ; this%wa_col (:) = nan + ! Yorck + allocate(this%wa_col_inc (begc:endc)) ; this%wa_col_inc (:) = nan + allocate(this%wa_col_mean (begc:endc)) ; this%wa_col_mean (:) = nan + allocate(this%qcharge_col (begc:endc)) ; this%qcharge_col (:) = nan allocate(this%fracice_col (begc:endc,nlevgrnd)) ; this%fracice_col (:,:) = nan allocate(this%icefrac_col (begc:endc,nlevgrnd)) ; this%icefrac_col (:,:) = nan diff --git a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 index a157e16bb6..4829a7a755 100644 --- a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 @@ -150,12 +150,23 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) real(r8) ,pointer :: gti (:) ! read in - fmax real(r8) ,pointer :: sand3d (:,:) ! read in - soil texture: percent sand (needs to be a pointer for use in ncdio) real(r8) ,pointer :: clay3d (:,:) ! read in - soil texture: percent clay (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: psis_sat (:,:) ! read in - soil parameter: sucsat (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: shape_param (:,:) ! read in - soil parameter: bsw (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: thetas (:,:) ! read in - soil parameter: watsat (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: ks (:,:) ! read in - soil parameter: xksat (needs to be a pointer for use in ncdio) + + real(r8) ,pointer :: psis_sat_adj (:,:) ! read in - soil parameter: sucsat (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: shape_param_adj (:,:) ! read in - soil parameter: bsw (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: thetas_adj (:,:) ! read in - soil parameter: watsat (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: ks_adj (:,:) ! read in - soil parameter: xksat (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: organic3d (:,:) ! read in - organic matter: kg/m3 (needs to be a pointer for use in ncdio) character(len=256) :: locfn ! local filename integer :: ipedof integer :: begp, endp integer :: begc, endc integer :: begg, endg + logical :: parameters_in_file, parameters_in_file_adj !----------------------------------------------------------------------- begp = bounds%begp; endp= bounds%endp @@ -225,6 +236,15 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) allocate(sand3d(begg:endg,nlevsoifl)) allocate(clay3d(begg:endg,nlevsoifl)) + allocate(thetas(begg:endg,nlevsoifl)) + allocate(shape_param(begg:endg,nlevsoifl)) + allocate(psis_sat(begg:endg,nlevsoifl)) + allocate(ks(begg:endg,nlevsoifl)) + + allocate(thetas_adj(begg:endg,nlevgrnd)) + allocate(shape_param_adj(begg:endg,nlevgrnd)) + allocate(psis_sat_adj(begg:endg,nlevgrnd)) + allocate(ks_adj(begg:endg,nlevgrnd)) ! Determine organic_max from parameter file @@ -262,6 +282,59 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) call endrun(msg=' ERROR: PCT_CLAY NOT on surfdata file'//errMsg(sourcefile, __LINE__)) end if + ! include option to also read hydraulic parameters from file. Keep it variable so that the code also works for surface files that were + ! generated without parameter perturbation and parameter as input variables + + call ncd_io(ncid=ncid, varname='THETAS', flag='read', data=thetas, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file = .False. + else + parameters_in_file = .True. + end if + + if (parameters_in_file) then ! read also other paramters, if one of them is not included, use sand and clay to compute parameters + ! via pedotransfer function + call ncd_io(ncid=ncid, varname='SHAPE_PARAM', flag='read', data=shape_param, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file = .False. + end if + + call ncd_io(ncid=ncid, varname='PSIS_SAT', flag='read', data=psis_sat, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file = .False. + end if + + call ncd_io(ncid=ncid, varname='KSAT', flag='read', data=ks, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file = .False. + end if + end if + + call ncd_io(ncid=ncid, varname='THETAS_adj', flag='read', data=thetas_adj, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file_adj = .False. + else + parameters_in_file_adj = .True. + end if + + if (parameters_in_file_adj) then ! read also other paramters, if one of them is not included, use sand and clay to compute parameters + ! via pedotransfer function + call ncd_io(ncid=ncid, varname='SHAPE_PARAM_adj', flag='read', data=shape_param_adj, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file_adj = .False. + end if + + call ncd_io(ncid=ncid, varname='PSIS_SAT_adj', flag='read', data=psis_sat_adj, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file_adj = .False. + end if + + call ncd_io(ncid=ncid, varname='KSAT_adj', flag='read', data=ks_adj, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file_adj = .False. + end if + end if + do p = begp,endp g = patch%gridcell(p) if ( sand3d(g,1)+clay3d(g,1) == 0.0_r8 )then @@ -454,6 +527,26 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) call pedotransf(ipedof, sand, clay, & soilstate_inst%watsat_col(c,lev), soilstate_inst%bsw_col(c,lev), soilstate_inst%sucsat_col(c,lev), xksat) + + ! if parameters are included in the file, watsat,... are overwritten with the values from there. If not, the pedotransfer + ! function is used + + if (parameters_in_file) then + if (lev <= nlevsoifl) then + ! Use values from the file for the soil layers + soilstate_inst%watsat_col(c,lev) = thetas(col%gridcell(c), lev) + soilstate_inst%bsw_col(c,lev) = shape_param(col%gridcell(c), lev) + soilstate_inst%sucsat_col(c,lev) = psis_sat(col%gridcell(c), lev) + xksat = ks(col%gridcell(c), lev) ! mm/s + else + ! Use the value from the 10th soil level as a default value + soilstate_inst%watsat_col(c,lev) = thetas(col%gridcell(c), 10) + soilstate_inst%bsw_col(c,lev) = shape_param(col%gridcell(c), 10) + soilstate_inst%sucsat_col(c,lev) = psis_sat(col%gridcell(c), 10) + xksat = ks(col%gridcell(c), 10) + end if + end if + om_watsat = max(0.93_r8 - 0.1_r8 *(zsoi(lev)/zsapric), 0.83_r8) om_b = min(2.7_r8 + 9.3_r8 *(zsoi(lev)/zsapric), 12.0_r8) om_sucsat = min(10.3_r8 - 0.2_r8 *(zsoi(lev)/zsapric), 10.1_r8) @@ -462,7 +555,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) soilstate_inst%bd_col(c,lev) = (1._r8 - soilstate_inst%watsat_col(c,lev))*2.7e3_r8 soilstate_inst%watsat_col(c,lev) = (1._r8 - om_frac) * soilstate_inst%watsat_col(c,lev) + om_watsat*om_frac tkm = (1._r8-om_frac) * (8.80_r8*sand+2.92_r8*clay)/(sand+clay)+om_tkm*om_frac ! W/(m K) - soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * (2.91_r8 + 0.159_r8*clay) + om_frac*om_b + soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * soilstate_inst%bsw_col(c,lev) + om_frac*om_b soilstate_inst%sucsat_col(c,lev) = (1._r8-om_frac) * soilstate_inst%sucsat_col(c,lev) + om_sucsat*om_frac soilstate_inst%hksat_min_col(c,lev) = xksat @@ -486,6 +579,17 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) end if soilstate_inst%hksat_col(c,lev) = uncon_frac*uncon_hksat + (perc_frac*om_frac)*om_hksat + + if (parameters_in_file_adj) then + ! Use values from the file for the soil layers + soilstate_inst%watsat_col(c,lev) = thetas_adj(col%gridcell(c), lev) + soilstate_inst%bsw_col(c,lev) = shape_param_adj(col%gridcell(c), lev) + soilstate_inst%sucsat_col(c,lev) = psis_sat_adj(col%gridcell(c), lev) + soilstate_inst%hksat_col(c,lev) = ks_adj(col%gridcell(c), lev) ! mm/s + end if + + + soilstate_inst%tkmg_col(c,lev) = tkm ** (1._r8- soilstate_inst%watsat_col(c,lev)) soilstate_inst%tksatu_col(c,lev) = soilstate_inst%tkmg_col(c,lev)*0.57_r8**soilstate_inst%watsat_col(c,lev) @@ -560,7 +664,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) tkm = (1._r8-om_frac)*(8.80_r8*sand+2.92_r8*clay)/(sand+clay) + om_tkm * om_frac ! W/(m K) - soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*(2.91_r8 + 0.159_r8*clay) + om_frac * om_b_lake + soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*soilstate_inst%bsw_col(c,lev) + om_frac * om_b_lake soilstate_inst%sucsat_col(c,lev) = (1._r8-om_frac)*soilstate_inst%sucsat_col(c,lev) + om_sucsat_lake * om_frac @@ -624,6 +728,8 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) deallocate(sand3d, clay3d, organic3d) deallocate(zisoifl, zsoifl, dzsoifl) + deallocate(thetas, shape_param, psis_sat, ks) + deallocate(thetas_adj, shape_param_adj, psis_sat_adj, ks_adj) end subroutine SoilStateInitTimeConst diff --git a/src/clm5/biogeophys/WaterStateType.F90 b/src/clm5/biogeophys/WaterStateType.F90 index 65944de240..2cd4a4a715 100644 --- a/src/clm5/biogeophys/WaterStateType.F90 +++ b/src/clm5/biogeophys/WaterStateType.F90 @@ -54,6 +54,57 @@ module WaterstateType real(r8), pointer :: ice1_grc (:) ! grc initial gridcell total h2o ice content real(r8), pointer :: ice2_grc (:) ! grc post land cover change total ice content real(r8), pointer :: tws_grc (:) ! grc total water storage (mm H2O) + + ! Yorck additions, variables for getting monthly means which can be compared to a GRACE measurements, other variables will just + ! provide instantaneous values + + real(r8), pointer :: h2osno_col_mean (:) ! col snow water (mm H2O) + real(r8), pointer :: h2osoi_liq_col_mean (:,:) ! col liquid water (kg/m2) (new) (-nlevsno+1:nlevgrnd) + real(r8), pointer :: h2osoi_ice_col_mean (:,:) ! col ice lens (kg/m2) (new) (-nlevsno+1:nlevgrnd) + real(r8), pointer :: h2ocan_patch_mean (:) ! patch canopy water (mm H2O) + real(r8), pointer :: snocan_patch_mean (:) ! patch canopy water (mm H2O) + real(r8), pointer :: h2osfc_col_mean (:) ! col surface water (mm H2O) + real(r8), pointer :: total_plant_stored_h2o_col_mean(:) ! col water that is bound in plants, including roots, sapwood, leaves, etc + ! in most cases, the vegetation scheme does not have a dynamic + ! water storage in plants, and thus 0.0 is a suitable for the trivial case. + ! When FATES is coupled in with plant hydraulics turned on, this storage + ! term is set to non-zero. (kg/m2 H2O) + + + real(r8), pointer :: tws_hactive (:) ! TWS for hydrological active columns + real(r8), pointer :: tws_hactive_mean (:) ! TWS for hydrological active columns + + + ! also increments for all variables are added + real(r8), pointer :: h2osno_col_inc (:) ! col snow water (mm H2O) + real(r8), pointer :: h2osoi_liq_col_inc (:,:) ! col liquid water (kg/m2) (new) (-nlevsno+1:nlevgrnd) + real(r8), pointer :: h2osoi_ice_col_inc (:,:) ! col ice lens (kg/m2) (new) (-nlevsno+1:nlevgrnd) + real(r8), pointer :: h2ocan_patch_inc (:) ! patch canopy water (mm H2O) + real(r8), pointer :: h2osfc_col_inc (:) ! col surface water (mm H2O) + real(r8), pointer :: total_plant_stored_h2o_col_inc (:) ! col water that is bound in plants, including roots, sapwood, leaves, etc + ! in most cases, the vegetation scheme does not have a dynamic + ! water storage in plants, and thus 0.0 is a suitable for the trivial case. + ! When FATES is coupled in with plant hydraulics turned on, this storage + ! term is set to non-zero. (kg/m2 H2O) + real(r8), pointer :: tws_grc_inc (:) ! grc total water storage (mm H2O) + + ! for analysis of the statevector, also variables for putting the state inside are initialized + real(r8), pointer :: tws_state_before (:) ! TWS state + real(r8), pointer :: h2osoi_liq_state_before (:,:) ! soil liq state + real(r8), pointer :: h2osoi_ice_state_before (:,:) ! soil ice state + real(r8), pointer :: h2osfc_state_before (:) ! surface water state + real(r8), pointer :: h2osno_state_before (:) ! snow state + real(r8), pointer :: h2ocan_state_before (:) ! canopy state + + real(r8), pointer :: tws_state_after (:) ! TWS state + real(r8), pointer :: h2osoi_liq_state_after (:,:) ! soil liq state + real(r8), pointer :: h2osoi_ice_state_after (:,:) ! soil ice state + real(r8), pointer :: h2osfc_state_after (:) ! surface water state + real(r8), pointer :: h2osno_state_after (:) ! snow state + real(r8), pointer :: h2ocan_state_after (:) ! canopy state + + ! END Yorck + #ifdef COUP_OAS_PFL real(r8), pointer :: pfl_psi_col (:,:) ! ParFlow pressure head COUP_OAS_PFL real(r8), pointer :: pfl_h2osoi_liq_col (:,:) ! ParFlow soil liquid COUP_OAS_PFL @@ -211,6 +262,41 @@ subroutine InitAllocate(this, bounds) allocate(this%ice2_grc (begg:endg)) ; this%ice2_grc (:) = nan allocate(this%tws_grc (begg:endg)) ; this%tws_grc (:) = nan + ! Yorck additions (see above) + allocate(this%h2osoi_ice_col_mean (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_ice_col_mean (:,:) = nan + allocate(this%h2osoi_liq_col_mean (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_liq_col_mean (:,:) = nan + allocate(this%h2osno_col_mean (begc:endc)) ; this%h2osno_col_mean (:) = nan + allocate(this%h2ocan_patch_mean (begp:endp)) ; this%h2ocan_patch_mean (:) = nan + allocate(this%snocan_patch_mean (begp:endp)) ; this%snocan_patch_mean (:) = nan + allocate(this%h2osfc_col_mean (begc:endc)) ; this%h2osfc_col_mean (:) = nan + allocate(this%total_plant_stored_h2o_col_mean(begc:endc)) ; this%total_plant_stored_h2o_col_mean(:) = nan + + allocate(this%tws_hactive (begg:endg)) ; this%tws_hactive (:) = nan + allocate(this%tws_hactive_mean (begg:endg)) ; this%tws_hactive_mean (:) = nan + + allocate(this%h2osoi_ice_col_inc (begc:endc,1:nlevsoi)) ; this%h2osoi_ice_col_inc (:,:) = nan + allocate(this%h2osoi_liq_col_inc (begc:endc,1:nlevsoi)) ; this%h2osoi_liq_col_inc (:,:) = nan + allocate(this%h2osno_col_inc (begc:endc)) ; this%h2osno_col_inc (:) = nan + allocate(this%h2ocan_patch_inc (begp:endp)) ; this%h2ocan_patch_inc (:) = nan + allocate(this%h2osfc_col_inc (begc:endc)) ; this%h2osfc_col_inc (:) = nan + allocate(this%total_plant_stored_h2o_col_inc (begc:endc)) ; this%total_plant_stored_h2o_col_inc (:) = nan + allocate(this%tws_grc_inc (begg:endg)) ; this%tws_grc_inc (:) = nan + + allocate(this%tws_state_before (begg:endg)) ; this%tws_state_before (:) = nan + allocate(this%h2osoi_liq_state_before(begg:endg,1:nlevsoi)) ; this%h2osoi_liq_state_before(:,:) = nan + allocate(this%h2osoi_ice_state_before(begg:endg,1:nlevsoi)) ; this%h2osoi_ice_state_before(:,:) = nan + allocate(this%h2osno_state_before (begg:endg)) ; this%h2osno_state_before (:) = nan + allocate(this%h2osfc_state_before (begg:endg)) ; this%h2osfc_state_before (:) = nan + allocate(this%h2ocan_state_before (begg:endg)) ; this%h2ocan_state_before (:) = nan + + allocate(this%tws_state_after (begg:endg)) ; this%tws_state_after (:) = nan + allocate(this%h2osoi_liq_state_after (begg:endg,1:nlevsoi)) ; this%h2osoi_liq_state_after (:,:) = nan + allocate(this%h2osoi_ice_state_after (begg:endg,1:nlevsoi)) ; this%h2osoi_ice_state_after (:,:) = nan + allocate(this%h2osno_state_after (begg:endg)) ; this%h2osno_state_after (:) = nan + allocate(this%h2osfc_state_after (begg:endg)) ; this%h2osfc_state_after (:) = nan + allocate(this%h2ocan_state_after (begg:endg)) ; this%h2ocan_state_after (:) = nan + ! END Yorck + allocate(this%total_plant_stored_h2o_col(begc:endc)) ; this%total_plant_stored_h2o_col(:) = nan allocate(this%snw_rds_col (begc:endc,-nlevsno+1:0)) ; this%snw_rds_col (:,:) = nan @@ -402,6 +488,12 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='total water storage', & ptr_lnd=this%tws_grc) + ! Yorck: add also TWS_hactive in history outputs + this%tws_hactive(begg:endg) = spval + call hist_addfld1d (fname='TWS_hactive', units='mm', & + avgflag='A', long_name='total water storage of hydrological active cells', & + ptr_lnd=this%tws_hactive) + ! (rgk 02-02-2017) There is intentionally no entry here for stored plant water ! I think that since the value is zero in all cases except ! for FATES plant hydraulics, it will be confusing for users diff --git a/src/clm5/cpl/lnd_comp_mct.F90 b/src/clm5/cpl/lnd_comp_mct.F90 index 56179beab6..3ccfc10483 100644 --- a/src/clm5/cpl/lnd_comp_mct.F90 +++ b/src/clm5/cpl/lnd_comp_mct.F90 @@ -725,7 +725,7 @@ subroutine lnd_handle_resume( infodata ) ! Otherwise restart was modified and we are resuming from data assimulation else resume_from_data_assim = .true. - write(iulog,*) 'resume_from_DA ', resume_from_data_assim + !write(iulog,*) 'resume_from_DA ', resume_from_data_assim end if if ( resume_from_data_assim ) call update_DA_nstep() diff --git a/src/clm5/main/atm2lndType.F90 b/src/clm5/main/atm2lndType.F90 index 524908fe3c..4f64d41348 100644 --- a/src/clm5/main/atm2lndType.F90 +++ b/src/clm5/main/atm2lndType.F90 @@ -147,6 +147,11 @@ module atm2lndType real(r8) , pointer :: t_mo_patch (:) => null() ! patch 30-day average temperature (Kelvin) real(r8) , pointer :: t_mo_min_patch (:) => null() ! patch annual min of t_mo (Kelvin) + + ! Yorck + real(r8), pointer :: volr_grc_inc (:) => null() ! rof volr total volume increment (m3) + real(r8), pointer :: volr_grc_mean (:) => null() ! rof volr total volume mean (m3) + contains procedure, public :: Init @@ -552,6 +557,11 @@ subroutine InitAllocate(this, bounds) allocate(this%pfl_psi_grc (begg:endg,1:nlevgrnd)); this%pfl_psi_grc (:,:) = ival allocate(this%pfl_h2osoi_liq_grc (begg:endg,1:nlevgrnd)); this%pfl_h2osoi_liq_grc (:,:) = ival #endif + + ! Yorck + allocate(this%volr_grc_inc (begg:endg)) ; this%volr_grc_inc (:) = ival + allocate(this%volr_grc_mean (begg:endg)) ; this%volr_grc_mean (:) = ival + ! anomaly forcing allocate(this%bc_precip_grc (begg:endg)) ; this%bc_precip_grc (:) = ival allocate(this%af_precip_grc (begg:endg)) ; this%af_precip_grc (:) = ival diff --git a/src/clm5/main/clm_driver.F90 b/src/clm5/main/clm_driver.F90 index 7fef829b08..e2c47fb181 100644 --- a/src/clm5/main/clm_driver.F90 +++ b/src/clm5/main/clm_driver.F90 @@ -1067,7 +1067,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro waterstate_inst, waterflux_inst, irrigation_inst, energyflux_inst, & solarabs_inst, drydepvel_inst, & vocemis_inst, fireemis_inst, dust_inst, ch4_inst, glc_behavior, & - lnd2atm_inst, & + lnd2atm_inst, soilhydrology_inst, soilstate_inst, & net_carbon_exchange_grc = net_carbon_exchange_grc(bounds_proc%begg:bounds_proc%endg)) deallocate(net_carbon_exchange_grc) call t_stopf('lnd2atm') diff --git a/src/clm5/main/clm_varcon.F90 b/src/clm5/main/clm_varcon.F90 index d0a2053568..dc4f776424 100644 --- a/src/clm5/main/clm_varcon.F90 +++ b/src/clm5/main/clm_varcon.F90 @@ -106,6 +106,14 @@ module clm_varcon ! Keep this negative to avoid conflicts with possible valid values integer , public, parameter :: ispval = -9999 ! special value for int data + + ! ------------------------------------------------------------------------ + ! Yorck variables + ! ------------------------------------------------------------------------ + + integer :: averaging_var = 0 ! averaging integer for averaging TWS variables according to GRACE interval + integer :: set_averaging_to_zero = ispval ! integer indicating at which point in time the averageTemp has to be set to zero due to missing observation values + ! ------------------------------------------------------------------------ ! These are tunable constants from clm2_3 ! ------------------------------------------------------------------------ diff --git a/src/clm5/main/lnd2atmMod.F90 b/src/clm5/main/lnd2atmMod.F90 index 950c8fd610..cfcf94e81e 100644 --- a/src/clm5/main/lnd2atmMod.F90 +++ b/src/clm5/main/lnd2atmMod.F90 @@ -12,7 +12,7 @@ module lnd2atmMod use shr_megan_mod , only : shr_megan_mechcomps_n use shr_fire_emis_mod , only : shr_fire_emis_mechcomps_n use clm_varpar , only : numrad, ndst, nlevgrnd, nlevsoi !ndst = number of dust bins. - use clm_varcon , only : rair, grav, cpair, hfus, tfrz, spval + use clm_varcon , only : rair, grav, cpair, hfus, tfrz, spval, set_averaging_to_zero, averaging_var, aquifer_water_baseline use clm_varctl , only : iulog, use_lch4 use seq_drydep_mod , only : n_drydep, drydep_method, DD_XLND use decompMod , only : bounds_type @@ -38,6 +38,11 @@ module lnd2atmMod use LandunitType , only : lun use GridcellType , only : grc use landunit_varcon , only : istice_mec, istsoil, istcrop + use clm_time_manager , only : get_nstep + use SoilHydrologyType , only : soilhydrology_type + use SoilStateType , only : soilstate_type + use PatchType , only : patch + ! ! !PUBLIC TYPES: implicit none @@ -126,7 +131,7 @@ subroutine lnd2atm(bounds, & waterstate_inst, waterflux_inst, irrigation_inst, energyflux_inst, & solarabs_inst, drydepvel_inst, & vocemis_inst, fireemis_inst, dust_inst, ch4_inst, glc_behavior, & - lnd2atm_inst, & + lnd2atm_inst, soilhydrology_inst, soilstate_inst, & net_carbon_exchange_grc) ! ! !DESCRIPTION: @@ -153,17 +158,32 @@ subroutine lnd2atm(bounds, & type(ch4_type) , intent(in) :: ch4_inst type(glc_behavior_type) , intent(in) :: glc_behavior type(lnd2atm_type) , intent(inout) :: lnd2atm_inst + ! Yorck + type(soilhydrology_type) , intent(inout) :: soilhydrology_inst + type(soilstate_type) , intent(inout) :: soilstate_inst + ! end Yorck real(r8) , intent(in) :: net_carbon_exchange_grc( bounds%begg: ) ! net carbon exchange between land and atmosphere, positive for source (gC/m2/s) ! ! !LOCAL VARIABLES: - integer :: c, g, j ! indices + integer :: c, g, j, p, l, index, counter ! indices real(r8) :: qflx_ice_runoff_col(bounds%begc:bounds%endc) ! total column-level ice runoff real(r8) :: eflx_sh_ice_to_liq_grc(bounds%begg:bounds%endg) ! sensible heat flux generated from the ice to liquid conversion, averaged to gridcell real(r8), parameter :: amC = 12.0_r8 ! Atomic mass number for Carbon real(r8), parameter :: amO = 16.0_r8 ! Atomic mass number for Oxygen real(r8), parameter :: amCO2 = amC + 2.0_r8*amO ! Atomic mass number for CO2 + real(r8), parameter :: m_per_mm = 1.e-3_r8 ! 0.001 meters per mm + real(r8), parameter :: sec_per_hr = 3600 ! 3600 s in 1 hour ! The following converts g of C to kg of CO2 real(r8), parameter :: convertgC2kgCO2 = 1.0e-3_r8 * (amCO2/amC) + integer :: nstep ! time step number + + REAL, allocatable :: h2osoi_liq_grc(:) + REAL, allocatable :: h2osoi_ice_grc(:) + REAL, allocatable :: h2osno_grc(:) + REAL, allocatable :: h2osfc_grc(:) + REAL, allocatable :: h2ocan_grc(:) + logical, allocatable :: found(:) + logical, allocatable :: found_patch(:) !------------------------------------------------------------------------ SHR_ASSERT_ALL((ubound(net_carbon_exchange_grc) == (/bounds%endg/)), errMsg(sourcefile, __LINE__)) @@ -442,6 +462,166 @@ subroutine lnd2atm(bounds, & enddo #endif + ! Yorck: update also TWS of hydrological active cells + + if (allocated(h2osoi_liq_grc)) deallocate(h2osoi_liq_grc) + allocate(h2osoi_liq_grc(bounds%begg:bounds%endg)) + h2osoi_liq_grc(bounds%begg:bounds%endg) = spval + + if (allocated(h2osoi_ice_grc)) deallocate(h2osoi_ice_grc) + allocate(h2osoi_ice_grc(bounds%begg:bounds%endg)) + h2osoi_ice_grc(bounds%begg:bounds%endg) = spval + + if (allocated(h2osno_grc)) deallocate(h2osno_grc) + allocate(h2osno_grc(bounds%begg:bounds%endg)) + h2osno_grc(bounds%begg:bounds%endg) = spval + + if (allocated(h2osfc_grc)) deallocate(h2osfc_grc) + allocate(h2osfc_grc(bounds%begg:bounds%endg)) + h2osfc_grc(bounds%begg:bounds%endg) = spval + + if (allocated(h2ocan_grc)) deallocate(h2ocan_grc) + allocate(h2ocan_grc(bounds%begg:bounds%endg)) + h2ocan_grc(bounds%begg:bounds%endg) = spval + + if (allocated(found)) deallocate(found) + allocate(found(bounds%begg:bounds%endg)) + found(bounds%begg:bounds%endg) = .false. + + if (allocated(found_patch)) deallocate(found_patch) + allocate(found_patch(bounds%begg:bounds%endg)) + found_patch(bounds%begg:bounds%endg) = .false. + + do g = bounds%begg, bounds%endg + counter = 0 + do c = bounds%begc, bounds%endc + if (col%hydrologically_active(c) .and. col%gridcell(c)==g) then + do j = 1,nlevsoi + if (j==1 .and. .not. found(g)) then + found(g) = .true. + h2osoi_liq_grc(g) = 0._r8 + h2osoi_ice_grc(g) = 0._r8 + h2osno_grc(g) = 0._r8 + h2osfc_grc(g) = 0._r8 + end if + if (j<=col%nbedrock(c)) then + h2osoi_liq_grc(g) = h2osoi_liq_grc(g) + waterstate_inst%h2osoi_liq_col(c,j) + h2osoi_ice_grc(g) = h2osoi_ice_grc(g) + waterstate_inst%h2osoi_ice_col(c,j) + end if + end do + + + + h2osno_grc(g) = h2osno_grc(g) + waterstate_inst%h2osno_col(c) + h2osfc_grc(g) = h2osfc_grc(g) + waterstate_inst%h2osfc_col(c) + counter= counter+1 + end if + end do + if (counter.ne.0) then + h2osoi_liq_grc(g) = h2osoi_liq_grc(g)/counter + h2osoi_ice_grc(g) = h2osoi_ice_grc(g)/counter + h2osno_grc(g) = h2osno_grc(g)/counter + h2osfc_grc(g) = h2osfc_grc(g)/counter + end if + + + counter = 0 + do p = bounds%begp, bounds%endp + c = patch%column(p) + if (col%hydrologically_active(c) .and. col%gridcell(c)==g .and. patch%active(p)) then + if (.not. found_patch(g)) then + found_patch(g) = .true. + h2ocan_grc(g) = 0 + end if + h2ocan_grc(g) = h2ocan_grc(g) + waterstate_inst%h2ocan_patch(p) + counter = counter+1 + end if + end do + if (counter.ne.0) then + h2ocan_grc(g) = h2ocan_grc(g)/counter + else + h2ocan_grc(g) = 0 + end if + + if (found(g)) then + waterstate_inst%tws_hactive(g) = h2osoi_liq_grc(g) + h2osoi_ice_grc(g) + h2osno_grc(g) + h2osfc_grc(g) + h2ocan_grc(g) + end if + end do + + if (allocated(h2osoi_liq_grc)) deallocate(h2osoi_liq_grc) + + if (allocated(h2osoi_ice_grc)) deallocate(h2osoi_ice_grc) + + if (allocated(h2osno_grc)) deallocate(h2osno_grc) + + if (allocated(h2osfc_grc)) deallocate(h2osfc_grc) + + if (allocated(h2ocan_grc)) deallocate(h2ocan_grc) + + if (allocated(found)) deallocate(found) + + if (allocated(found_patch)) deallocate(found_patch) + + + ! YORCK --> idea: when TWS is calculated at the end of a time step, update all running mean variables. + ! this is done for all compartments of TWS, the variables are initialized in the respective type modules + ! here because the update of all important variables inside CLM is finished! + nstep = get_nstep() + if (nstep.eq.set_averaging_to_zero) then + + averaging_var = 0 + + end if + + averaging_var = averaging_var+1 + + ! if averaging var was resetted, reset also running averages + if (averaging_var == 1) then + do c = bounds%begc, bounds%endc + waterstate_inst%h2osno_col_mean(c) = 0._r8 + soilhydrology_inst%wa_col_mean(c) = 0._r8 + waterstate_inst%h2osfc_col_mean(c) = 0._r8 + waterstate_inst%total_plant_stored_h2o_col_mean(c) = 0._r8 + do j = 1,nlevgrnd + waterstate_inst%h2osoi_liq_col_mean(c,j) = 0._r8 + waterstate_inst%h2osoi_ice_col_mean(c,j) = 0._r8 + end do + end do + do g = bounds%begg, bounds%endg + atm2lnd_inst%volr_grc_mean(g) = 0._r8 + waterstate_inst%tws_hactive_mean(g) = 0._r8 + end do + + do p = bounds%begp, bounds%endp + waterstate_inst%h2ocan_patch_mean(p) = 0._r8 + waterstate_inst%snocan_patch_mean(p) = 0._r8 + end do + end if + + ! Yorck + ! update running average with mean = (var + (count-1)*mean)/count for all variables that contribute to TWS + do c = bounds%begc, bounds%endc + waterstate_inst%h2osno_col_mean(c) = (waterstate_inst%h2osno_col(c)+(averaging_var-1)*waterstate_inst%h2osno_col_mean(c))/averaging_var + soilhydrology_inst%wa_col_mean(c) = (soilhydrology_inst%wa_col(c)+(averaging_var-1)*soilhydrology_inst%wa_col_mean(c))/averaging_var + waterstate_inst%h2osfc_col_mean(c) = (waterstate_inst%h2osfc_col(c)+(averaging_var-1)*waterstate_inst%h2osfc_col_mean(c))/averaging_var + waterstate_inst%total_plant_stored_h2o_col_mean(c) = (waterstate_inst%total_plant_stored_h2o_col(c)+(averaging_var-1) * & + waterstate_inst%total_plant_stored_h2o_col_mean(c))/averaging_var + do j = 1,nlevgrnd + waterstate_inst%h2osoi_liq_col_mean(c,j) = (waterstate_inst%h2osoi_liq_col(c,j)+(averaging_var-1)* waterstate_inst%h2osoi_liq_col_mean(c,j))/averaging_var + waterstate_inst%h2osoi_ice_col_mean(c,j) = (waterstate_inst%h2osoi_ice_col(c,j)+(averaging_var-1)* waterstate_inst%h2osoi_ice_col_mean(c,j))/averaging_var + end do + end do + + do g = bounds%begg, bounds%endg + atm2lnd_inst%volr_grc_mean(g) = (atm2lnd_inst%volr_grc(g)+(averaging_var-1)*atm2lnd_inst%volr_grc_mean(g))/averaging_var + waterstate_inst%tws_hactive_mean(g) = (waterstate_inst%tws_hactive(g)+(averaging_var-1)*waterstate_inst%tws_hactive_mean(g))/averaging_var + end do + + do p = bounds%begp, bounds%endp + waterstate_inst%h2ocan_patch_mean(p) = (waterstate_inst%h2ocan_patch(p)+(averaging_var-1)*waterstate_inst%h2ocan_patch_mean(p))/averaging_var + waterstate_inst%snocan_patch_mean(p) = (waterstate_inst%snocan_patch(p)+(averaging_var-1)*waterstate_inst%snocan_patch_mean(p))/averaging_var + end do + end subroutine lnd2atm !----------------------------------------------------------------------- diff --git a/src/csm_share/streams/shr_dmodel_mod.F90 b/src/csm_share/streams/shr_dmodel_mod.F90 index 394d0d52af..65c227ea0b 100644 --- a/src/csm_share/streams/shr_dmodel_mod.F90 +++ b/src/csm_share/streams/shr_dmodel_mod.F90 @@ -929,7 +929,6 @@ subroutine shr_dmodel_readstrm(stream, pio_subsystem, pio_iotype, pio_iodesc, gs ! 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 diff --git a/src/csm_share/streams/shr_stream_mod.F90 b/src/csm_share/streams/shr_stream_mod.F90 index f99b016f52..1717efbebb 100644 --- a/src/csm_share/streams/shr_stream_mod.F90 +++ b/src/csm_share/streams/shr_stream_mod.F90 @@ -1645,13 +1645,12 @@ subroutine shr_stream_readTCoord(strm,k,rc) if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_inquire_dimension') deallocate(dids) + ! Yorck ! 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... + ! After nt is adjusted, also the read in number of timesteps has to be adjusted. call shr_stream_getModelFieldName(strm,1,mfldName) if (INDEX(mfldName, "noise") .ne. 0) then @@ -3404,6 +3403,7 @@ subroutine shr_stream_bcast(stream,comm,rc) end subroutine shr_stream_bcast + ! Yorck subroutine shr_stream_get_dt_and_caseId(stream,dt,caseId) implicit none diff --git a/src/datm/datm_comp_mod.F90 b/src/datm/datm_comp_mod.F90 index af76859881..e68a141abe 100644 --- a/src/datm/datm_comp_mod.F90 +++ b/src/datm/datm_comp_mod.F90 @@ -937,35 +937,35 @@ 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)) + !!! begin perturbation block (Yorck) --> perturb temperature, long wave radiation, short wave radiation and precipitation + ! --> further variables can be introduced (changes also have to be made in reading in routines of streams for this, for now only + ! these variables are included) + ! tbot + if (stbot_noise > 0) then + a2x%rAttr(ktbot,n) = a2x%rAttr(ktbot,n) + avstrm%rAttr(stbot_noise,n) + end if - ! already here perturbation of lwdn as I do not understand where this is used in the code (Yorck) + ! lwdn 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 + ! swdn + if (sswdn_noise>0) then + avstrm%rAttr(sswdn,n) = avstrm%rAttr(sswdn,n) * avstrm%rAttr(sswdn_noise,n) + end if + + ! sprecn + if (sprecn_noise > 0) then + avstrm%rAttr(sprecn,n) = avstrm%rAttr(sprecn,n)*avstrm%rAttr(sprecn_noise,n) + end if + + !!! end perturbation block + + ! Limit very cold forcing to 180K + a2x%rAttr(ktbot,n) = max(180._r8, a2x%rAttr(ktbot,n)) a2x%rAttr(kptem,n) = a2x%rAttr(ktbot,n) !--- pressure --- @@ -1027,44 +1027,18 @@ 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 - - ! 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 + swndr = avstrm%rAttr(sswdn,n) * 0.50_R8 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 - - 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 + swndf = avstrm%rAttr(sswdn,n) * 0.50_R8 a2x%rAttr(kswndf,n) = (1._R8 - ratio_rvrf)*swndf - 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 + swvdr = avstrm%rAttr(sswdn,n) * 0.50_R8 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 - 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 + swvdf = avstrm%rAttr(sswdn,n) * 0.50_R8 a2x%rAttr(kswvdf,n) = (1._R8 - ratio_rvrf)*swvdf else call shr_sys_abort(subName//'ERROR: cannot compute short-wave down') @@ -1085,21 +1059,8 @@ 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 - ! 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 + a2x%rAttr(krc,n) = avstrm%rAttr(sprecn,n)*0.1_R8 + a2x%rAttr(krl,n) = avstrm%rAttr(sprecn,n)*0.9_R8 else call shr_sys_abort(subName//'ERROR: cannot compute rain and snow') endif From 58e64f19c49e69782e9995e1d607692a697c6ca9 Mon Sep 17 00:00:00 2001 From: Yorck Ewerdwalbesloh Date: Fri, 25 Jul 2025 11:01:30 +0200 Subject: [PATCH 03/21] perturbation routines refactored --- src/csm_share/streams/shr_dmodel_mod.F90 | 1 - src/csm_share/streams/shr_stream_mod.F90 | 6 +- src/datm/datm_comp_mod.F90 | 95 +++++++----------------- 3 files changed, 31 insertions(+), 71 deletions(-) diff --git a/src/csm_share/streams/shr_dmodel_mod.F90 b/src/csm_share/streams/shr_dmodel_mod.F90 index 394d0d52af..65c227ea0b 100644 --- a/src/csm_share/streams/shr_dmodel_mod.F90 +++ b/src/csm_share/streams/shr_dmodel_mod.F90 @@ -929,7 +929,6 @@ subroutine shr_dmodel_readstrm(stream, pio_subsystem, pio_iotype, pio_iodesc, gs ! 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 diff --git a/src/csm_share/streams/shr_stream_mod.F90 b/src/csm_share/streams/shr_stream_mod.F90 index f99b016f52..1717efbebb 100644 --- a/src/csm_share/streams/shr_stream_mod.F90 +++ b/src/csm_share/streams/shr_stream_mod.F90 @@ -1645,13 +1645,12 @@ subroutine shr_stream_readTCoord(strm,k,rc) if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_inquire_dimension') deallocate(dids) + ! Yorck ! 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... + ! After nt is adjusted, also the read in number of timesteps has to be adjusted. call shr_stream_getModelFieldName(strm,1,mfldName) if (INDEX(mfldName, "noise") .ne. 0) then @@ -3404,6 +3403,7 @@ subroutine shr_stream_bcast(stream,comm,rc) end subroutine shr_stream_bcast + ! Yorck subroutine shr_stream_get_dt_and_caseId(stream,dt,caseId) implicit none diff --git a/src/datm/datm_comp_mod.F90 b/src/datm/datm_comp_mod.F90 index af76859881..e68a141abe 100644 --- a/src/datm/datm_comp_mod.F90 +++ b/src/datm/datm_comp_mod.F90 @@ -937,35 +937,35 @@ 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)) + !!! begin perturbation block (Yorck) --> perturb temperature, long wave radiation, short wave radiation and precipitation + ! --> further variables can be introduced (changes also have to be made in reading in routines of streams for this, for now only + ! these variables are included) + ! tbot + if (stbot_noise > 0) then + a2x%rAttr(ktbot,n) = a2x%rAttr(ktbot,n) + avstrm%rAttr(stbot_noise,n) + end if - ! already here perturbation of lwdn as I do not understand where this is used in the code (Yorck) + ! lwdn 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 + ! swdn + if (sswdn_noise>0) then + avstrm%rAttr(sswdn,n) = avstrm%rAttr(sswdn,n) * avstrm%rAttr(sswdn_noise,n) + end if + + ! sprecn + if (sprecn_noise > 0) then + avstrm%rAttr(sprecn,n) = avstrm%rAttr(sprecn,n)*avstrm%rAttr(sprecn_noise,n) + end if + + !!! end perturbation block + + ! Limit very cold forcing to 180K + a2x%rAttr(ktbot,n) = max(180._r8, a2x%rAttr(ktbot,n)) a2x%rAttr(kptem,n) = a2x%rAttr(ktbot,n) !--- pressure --- @@ -1027,44 +1027,18 @@ 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 - - ! 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 + swndr = avstrm%rAttr(sswdn,n) * 0.50_R8 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 - - 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 + swndf = avstrm%rAttr(sswdn,n) * 0.50_R8 a2x%rAttr(kswndf,n) = (1._R8 - ratio_rvrf)*swndf - 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 + swvdr = avstrm%rAttr(sswdn,n) * 0.50_R8 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 - 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 + swvdf = avstrm%rAttr(sswdn,n) * 0.50_R8 a2x%rAttr(kswvdf,n) = (1._R8 - ratio_rvrf)*swvdf else call shr_sys_abort(subName//'ERROR: cannot compute short-wave down') @@ -1085,21 +1059,8 @@ 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 - ! 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 + a2x%rAttr(krc,n) = avstrm%rAttr(sprecn,n)*0.1_R8 + a2x%rAttr(krl,n) = avstrm%rAttr(sprecn,n)*0.9_R8 else call shr_sys_abort(subName//'ERROR: cannot compute rain and snow') endif From 0a556d4e04dcecfdf71dcb96fc79efd0a7b57949 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Wed, 15 Oct 2025 09:31:04 +0200 Subject: [PATCH 04/21] update SoilStateInitTimeConstMod for parameter perturbaton perturbation of the soil hydraulic parameters author: Yorck Ewerdwalbesloh --- .../biogeophys/SoilStateInitTimeConstMod.F90 | 110 +++++++++++++++++- 1 file changed, 108 insertions(+), 2 deletions(-) diff --git a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 index a157e16bb6..4829a7a755 100644 --- a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 @@ -150,12 +150,23 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) real(r8) ,pointer :: gti (:) ! read in - fmax real(r8) ,pointer :: sand3d (:,:) ! read in - soil texture: percent sand (needs to be a pointer for use in ncdio) real(r8) ,pointer :: clay3d (:,:) ! read in - soil texture: percent clay (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: psis_sat (:,:) ! read in - soil parameter: sucsat (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: shape_param (:,:) ! read in - soil parameter: bsw (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: thetas (:,:) ! read in - soil parameter: watsat (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: ks (:,:) ! read in - soil parameter: xksat (needs to be a pointer for use in ncdio) + + real(r8) ,pointer :: psis_sat_adj (:,:) ! read in - soil parameter: sucsat (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: shape_param_adj (:,:) ! read in - soil parameter: bsw (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: thetas_adj (:,:) ! read in - soil parameter: watsat (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: ks_adj (:,:) ! read in - soil parameter: xksat (needs to be a pointer for use in ncdio) + real(r8) ,pointer :: organic3d (:,:) ! read in - organic matter: kg/m3 (needs to be a pointer for use in ncdio) character(len=256) :: locfn ! local filename integer :: ipedof integer :: begp, endp integer :: begc, endc integer :: begg, endg + logical :: parameters_in_file, parameters_in_file_adj !----------------------------------------------------------------------- begp = bounds%begp; endp= bounds%endp @@ -225,6 +236,15 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) allocate(sand3d(begg:endg,nlevsoifl)) allocate(clay3d(begg:endg,nlevsoifl)) + allocate(thetas(begg:endg,nlevsoifl)) + allocate(shape_param(begg:endg,nlevsoifl)) + allocate(psis_sat(begg:endg,nlevsoifl)) + allocate(ks(begg:endg,nlevsoifl)) + + allocate(thetas_adj(begg:endg,nlevgrnd)) + allocate(shape_param_adj(begg:endg,nlevgrnd)) + allocate(psis_sat_adj(begg:endg,nlevgrnd)) + allocate(ks_adj(begg:endg,nlevgrnd)) ! Determine organic_max from parameter file @@ -262,6 +282,59 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) call endrun(msg=' ERROR: PCT_CLAY NOT on surfdata file'//errMsg(sourcefile, __LINE__)) end if + ! include option to also read hydraulic parameters from file. Keep it variable so that the code also works for surface files that were + ! generated without parameter perturbation and parameter as input variables + + call ncd_io(ncid=ncid, varname='THETAS', flag='read', data=thetas, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file = .False. + else + parameters_in_file = .True. + end if + + if (parameters_in_file) then ! read also other paramters, if one of them is not included, use sand and clay to compute parameters + ! via pedotransfer function + call ncd_io(ncid=ncid, varname='SHAPE_PARAM', flag='read', data=shape_param, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file = .False. + end if + + call ncd_io(ncid=ncid, varname='PSIS_SAT', flag='read', data=psis_sat, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file = .False. + end if + + call ncd_io(ncid=ncid, varname='KSAT', flag='read', data=ks, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file = .False. + end if + end if + + call ncd_io(ncid=ncid, varname='THETAS_adj', flag='read', data=thetas_adj, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file_adj = .False. + else + parameters_in_file_adj = .True. + end if + + if (parameters_in_file_adj) then ! read also other paramters, if one of them is not included, use sand and clay to compute parameters + ! via pedotransfer function + call ncd_io(ncid=ncid, varname='SHAPE_PARAM_adj', flag='read', data=shape_param_adj, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file_adj = .False. + end if + + call ncd_io(ncid=ncid, varname='PSIS_SAT_adj', flag='read', data=psis_sat_adj, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file_adj = .False. + end if + + call ncd_io(ncid=ncid, varname='KSAT_adj', flag='read', data=ks_adj, dim1name=grlnd, readvar=readvar) + if (.not. readvar) then + parameters_in_file_adj = .False. + end if + end if + do p = begp,endp g = patch%gridcell(p) if ( sand3d(g,1)+clay3d(g,1) == 0.0_r8 )then @@ -454,6 +527,26 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) call pedotransf(ipedof, sand, clay, & soilstate_inst%watsat_col(c,lev), soilstate_inst%bsw_col(c,lev), soilstate_inst%sucsat_col(c,lev), xksat) + + ! if parameters are included in the file, watsat,... are overwritten with the values from there. If not, the pedotransfer + ! function is used + + if (parameters_in_file) then + if (lev <= nlevsoifl) then + ! Use values from the file for the soil layers + soilstate_inst%watsat_col(c,lev) = thetas(col%gridcell(c), lev) + soilstate_inst%bsw_col(c,lev) = shape_param(col%gridcell(c), lev) + soilstate_inst%sucsat_col(c,lev) = psis_sat(col%gridcell(c), lev) + xksat = ks(col%gridcell(c), lev) ! mm/s + else + ! Use the value from the 10th soil level as a default value + soilstate_inst%watsat_col(c,lev) = thetas(col%gridcell(c), 10) + soilstate_inst%bsw_col(c,lev) = shape_param(col%gridcell(c), 10) + soilstate_inst%sucsat_col(c,lev) = psis_sat(col%gridcell(c), 10) + xksat = ks(col%gridcell(c), 10) + end if + end if + om_watsat = max(0.93_r8 - 0.1_r8 *(zsoi(lev)/zsapric), 0.83_r8) om_b = min(2.7_r8 + 9.3_r8 *(zsoi(lev)/zsapric), 12.0_r8) om_sucsat = min(10.3_r8 - 0.2_r8 *(zsoi(lev)/zsapric), 10.1_r8) @@ -462,7 +555,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) soilstate_inst%bd_col(c,lev) = (1._r8 - soilstate_inst%watsat_col(c,lev))*2.7e3_r8 soilstate_inst%watsat_col(c,lev) = (1._r8 - om_frac) * soilstate_inst%watsat_col(c,lev) + om_watsat*om_frac tkm = (1._r8-om_frac) * (8.80_r8*sand+2.92_r8*clay)/(sand+clay)+om_tkm*om_frac ! W/(m K) - soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * (2.91_r8 + 0.159_r8*clay) + om_frac*om_b + soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * soilstate_inst%bsw_col(c,lev) + om_frac*om_b soilstate_inst%sucsat_col(c,lev) = (1._r8-om_frac) * soilstate_inst%sucsat_col(c,lev) + om_sucsat*om_frac soilstate_inst%hksat_min_col(c,lev) = xksat @@ -486,6 +579,17 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) end if soilstate_inst%hksat_col(c,lev) = uncon_frac*uncon_hksat + (perc_frac*om_frac)*om_hksat + + if (parameters_in_file_adj) then + ! Use values from the file for the soil layers + soilstate_inst%watsat_col(c,lev) = thetas_adj(col%gridcell(c), lev) + soilstate_inst%bsw_col(c,lev) = shape_param_adj(col%gridcell(c), lev) + soilstate_inst%sucsat_col(c,lev) = psis_sat_adj(col%gridcell(c), lev) + soilstate_inst%hksat_col(c,lev) = ks_adj(col%gridcell(c), lev) ! mm/s + end if + + + soilstate_inst%tkmg_col(c,lev) = tkm ** (1._r8- soilstate_inst%watsat_col(c,lev)) soilstate_inst%tksatu_col(c,lev) = soilstate_inst%tkmg_col(c,lev)*0.57_r8**soilstate_inst%watsat_col(c,lev) @@ -560,7 +664,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) tkm = (1._r8-om_frac)*(8.80_r8*sand+2.92_r8*clay)/(sand+clay) + om_tkm * om_frac ! W/(m K) - soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*(2.91_r8 + 0.159_r8*clay) + om_frac * om_b_lake + soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*soilstate_inst%bsw_col(c,lev) + om_frac * om_b_lake soilstate_inst%sucsat_col(c,lev) = (1._r8-om_frac)*soilstate_inst%sucsat_col(c,lev) + om_sucsat_lake * om_frac @@ -624,6 +728,8 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) deallocate(sand3d, clay3d, organic3d) deallocate(zisoifl, zsoifl, dzsoifl) + deallocate(thetas, shape_param, psis_sat, ks) + deallocate(thetas_adj, shape_param_adj, psis_sat_adj, ks_adj) end subroutine SoilStateInitTimeConst From 9c4761d3d59c077bbfc1721439f0274ced44a0e8 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Mon, 10 Nov 2025 13:27:21 +0100 Subject: [PATCH 05/21] Set perturbation routines under USE_PDAF preparation for PR into master --- .../biogeophys/SoilStateInitTimeConstMod.F90 | 31 +++-- src/csm_share/streams/shr_dmodel_mod.F90 | 13 +- src/csm_share/streams/shr_strdata_mod.F90 | 123 ++++++++++++------ src/csm_share/streams/shr_stream_mod.F90 | 54 ++++++-- src/datm/datm_comp_mod.F90 | 45 ++++--- 5 files changed, 184 insertions(+), 82 deletions(-) diff --git a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 index 4829a7a755..7806bd3c38 100644 --- a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 @@ -150,6 +150,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) real(r8) ,pointer :: gti (:) ! read in - fmax real(r8) ,pointer :: sand3d (:,:) ! read in - soil texture: percent sand (needs to be a pointer for use in ncdio) real(r8) ,pointer :: clay3d (:,:) ! read in - soil texture: percent clay (needs to be a pointer for use in ncdio) +#ifdef USE_PDAF real(r8) ,pointer :: psis_sat (:,:) ! read in - soil parameter: sucsat (needs to be a pointer for use in ncdio) real(r8) ,pointer :: shape_param (:,:) ! read in - soil parameter: bsw (needs to be a pointer for use in ncdio) real(r8) ,pointer :: thetas (:,:) ! read in - soil parameter: watsat (needs to be a pointer for use in ncdio) @@ -159,14 +160,16 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) real(r8) ,pointer :: shape_param_adj (:,:) ! read in - soil parameter: bsw (needs to be a pointer for use in ncdio) real(r8) ,pointer :: thetas_adj (:,:) ! read in - soil parameter: watsat (needs to be a pointer for use in ncdio) real(r8) ,pointer :: ks_adj (:,:) ! read in - soil parameter: xksat (needs to be a pointer for use in ncdio) - +#endif real(r8) ,pointer :: organic3d (:,:) ! read in - organic matter: kg/m3 (needs to be a pointer for use in ncdio) character(len=256) :: locfn ! local filename integer :: ipedof integer :: begp, endp integer :: begc, endc integer :: begg, endg +#ifdef USE_PDAF logical :: parameters_in_file, parameters_in_file_adj +#endif !----------------------------------------------------------------------- begp = bounds%begp; endp= bounds%endp @@ -236,6 +239,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) allocate(sand3d(begg:endg,nlevsoifl)) allocate(clay3d(begg:endg,nlevsoifl)) +#ifdef USE_PDAF allocate(thetas(begg:endg,nlevsoifl)) allocate(shape_param(begg:endg,nlevsoifl)) allocate(psis_sat(begg:endg,nlevsoifl)) @@ -245,6 +249,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) allocate(shape_param_adj(begg:endg,nlevgrnd)) allocate(psis_sat_adj(begg:endg,nlevgrnd)) allocate(ks_adj(begg:endg,nlevgrnd)) +#ifdef ! Determine organic_max from parameter file @@ -282,6 +287,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) call endrun(msg=' ERROR: PCT_CLAY NOT on surfdata file'//errMsg(sourcefile, __LINE__)) end if +#ifdef USE_PDAF ! include option to also read hydraulic parameters from file. Keep it variable so that the code also works for surface files that were ! generated without parameter perturbation and parameter as input variables @@ -334,7 +340,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) parameters_in_file_adj = .False. end if end if - +#endif do p = begp,endp g = patch%gridcell(p) if ( sand3d(g,1)+clay3d(g,1) == 0.0_r8 )then @@ -527,7 +533,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) call pedotransf(ipedof, sand, clay, & soilstate_inst%watsat_col(c,lev), soilstate_inst%bsw_col(c,lev), soilstate_inst%sucsat_col(c,lev), xksat) - +#ifdef USE_PDAF ! if parameters are included in the file, watsat,... are overwritten with the values from there. If not, the pedotransfer ! function is used @@ -546,7 +552,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) xksat = ks(col%gridcell(c), 10) end if end if - +#endif om_watsat = max(0.93_r8 - 0.1_r8 *(zsoi(lev)/zsapric), 0.83_r8) om_b = min(2.7_r8 + 9.3_r8 *(zsoi(lev)/zsapric), 12.0_r8) om_sucsat = min(10.3_r8 - 0.2_r8 *(zsoi(lev)/zsapric), 10.1_r8) @@ -555,7 +561,11 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) soilstate_inst%bd_col(c,lev) = (1._r8 - soilstate_inst%watsat_col(c,lev))*2.7e3_r8 soilstate_inst%watsat_col(c,lev) = (1._r8 - om_frac) * soilstate_inst%watsat_col(c,lev) + om_watsat*om_frac tkm = (1._r8-om_frac) * (8.80_r8*sand+2.92_r8*clay)/(sand+clay)+om_tkm*om_frac ! W/(m K) +#ifndef USE_PDAF + soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * (2.91_r8 + 0.159_r8*clay) + om_frac*om_b +#else soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * soilstate_inst%bsw_col(c,lev) + om_frac*om_b +#endif soilstate_inst%sucsat_col(c,lev) = (1._r8-om_frac) * soilstate_inst%sucsat_col(c,lev) + om_sucsat*om_frac soilstate_inst%hksat_min_col(c,lev) = xksat @@ -579,7 +589,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) end if soilstate_inst%hksat_col(c,lev) = uncon_frac*uncon_hksat + (perc_frac*om_frac)*om_hksat - +#ifdef USE_PDAF if (parameters_in_file_adj) then ! Use values from the file for the soil layers soilstate_inst%watsat_col(c,lev) = thetas_adj(col%gridcell(c), lev) @@ -587,9 +597,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) soilstate_inst%sucsat_col(c,lev) = psis_sat_adj(col%gridcell(c), lev) soilstate_inst%hksat_col(c,lev) = ks_adj(col%gridcell(c), lev) ! mm/s end if - - - +#endif soilstate_inst%tkmg_col(c,lev) = tkm ** (1._r8- soilstate_inst%watsat_col(c,lev)) soilstate_inst%tksatu_col(c,lev) = soilstate_inst%tkmg_col(c,lev)*0.57_r8**soilstate_inst%watsat_col(c,lev) @@ -663,8 +671,11 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) soilstate_inst%watsat_col(c,lev) = (1._r8 - om_frac)*soilstate_inst%watsat_col(c,lev) + om_watsat_lake * om_frac tkm = (1._r8-om_frac)*(8.80_r8*sand+2.92_r8*clay)/(sand+clay) + om_tkm * om_frac ! W/(m K) - +#ifndef USE_PDAF + soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*(2.91_r8 + 0.159_r8*clay) + om_frac * om_b_lake +#else soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*soilstate_inst%bsw_col(c,lev) + om_frac * om_b_lake +#endif soilstate_inst%sucsat_col(c,lev) = (1._r8-om_frac)*soilstate_inst%sucsat_col(c,lev) + om_sucsat_lake * om_frac @@ -728,8 +739,10 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) deallocate(sand3d, clay3d, organic3d) deallocate(zisoifl, zsoifl, dzsoifl) +#ifdef USE_PDAF deallocate(thetas, shape_param, psis_sat, ks) deallocate(thetas_adj, shape_param_adj, psis_sat_adj, ks_adj) +#endif end subroutine SoilStateInitTimeConst diff --git a/src/csm_share/streams/shr_dmodel_mod.F90 b/src/csm_share/streams/shr_dmodel_mod.F90 index 65c227ea0b..92809378fa 100644 --- a/src/csm_share/streams/shr_dmodel_mod.F90 +++ b/src/csm_share/streams/shr_dmodel_mod.F90 @@ -764,12 +764,13 @@ subroutine shr_dmodel_readstrm(stream, pio_subsystem, pio_iotype, pio_iodesc, gs logical :: fileopen character(CL) :: currfile - ! Yorck: +#ifdef USE_PDAF + ! Yorck: character(CL) :: mfldName integer :: position ! integer to look if string contains noise integer :: caseId ! which ensemble member? integer :: dt ! time sampling of forcings - +#endif integer(in) :: ndims integer(in),pointer :: dimid(:) type(file_desc_t) :: pioid @@ -916,13 +917,16 @@ 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 +#ifdef USE_PDAF ! Yorck adjustment if (my_task == master_task) then call shr_stream_getModelFieldName(stream,k,mfldName) end if call shr_mpi_bcast(mfldName,mpicom,'mfldName') +#endif call shr_mpi_bcast(sfldName,mpicom,'sfldName') rcode = pio_inq_varid(pioid,trim(sfldName),varid) +#ifdef USE_PDAF !------------------------------------------------------------------------------------ ! Yorck adjustments : ! if noise paramter, the frame has to be adjusted depending on the ensemble member @@ -932,7 +936,9 @@ subroutine shr_dmodel_readstrm(stream, pio_subsystem, pio_iotype, pio_iodesc, gs ! idea: plug time and ensemble member into stream if (INDEX(mfldName, "noise") == 0) then - frame = nt +#endif + frame = nt +#ifdef USE_PDAF else if (my_task == master_task) then call shr_stream_get_dt_and_caseId(stream,dt,caseId) @@ -946,6 +952,7 @@ subroutine shr_dmodel_readstrm(stream, pio_subsystem, pio_iotype, pio_iodesc, gs end if frame = nt+caseId*24/dt end if +#endif call pio_setframe(pioid,varid,frame) call pio_read_darray(pioid,varid,pio_iodesc,av%rattr(k,:),rcode) enddo diff --git a/src/csm_share/streams/shr_strdata_mod.F90 b/src/csm_share/streams/shr_strdata_mod.F90 index 60f2abc9d7..83c83c8dfc 100644 --- a/src/csm_share/streams/shr_strdata_mod.F90 +++ b/src/csm_share/streams/shr_strdata_mod.F90 @@ -1385,17 +1385,20 @@ 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 --- - caseId, dt, numEns, 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 --- +#ifdef USE_PDAF + caseId, dt, numEns, & +#endif + nzg, domZvarName, & + taxMode, dtlimit, tintalgo, readmode, & + fillalgo, fillmask, fillread, fillwrite, & + mapalgo, mapmask, mapread, mapwrite, & + calendar) implicit none @@ -1441,11 +1444,12 @@ 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 +#ifdef USE_PDAF + ! Yorck integer(IN) ,optional ,intent(in) :: caseId integer(IN) ,optional ,intent(in) :: dt integer(IN) ,optional ,intent(in) :: numEns - +#endif !EOP ! --- local --- @@ -1508,37 +1512,61 @@ 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 +#ifndef USE_PDAF + 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 + ! Yorck adjustments: optionality of caseId, dt and numEns 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)) + 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)) + 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 +#endif + else +#ifndef USE_PDAF + 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)) +#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)) + 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)) + 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 +#endif + endif if (present(nzg)) then call shr_strdata_init(SDAT, mpicom, compid, & @@ -1558,7 +1586,10 @@ subroutine shr_strdata_create_oldway(SDAT, name, mpicom, compid, gsmap, ggrid, n filePath, filename, fldListFile, fldListModel, & pio_subsystem, pio_iotype, & !--- strdata optional --- - caseId, dt, numEns, nzg, domZvarName, & +#ifdef USE_PDAF + caseId, dt, numEns, & +#endif + nzg, domZvarName, & taxMode, dtlimit, tintalgo, readmode, & fillalgo, fillmask, fillread, fillwrite, & mapalgo, mapmask, mapread, mapwrite, & @@ -1611,21 +1642,33 @@ 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 +#ifdef USE_PDAF ! Yorck integer(IN) ,optional ,intent(in) :: caseId integer(IN) ,optional ,intent(in) :: dt integer(IN) ,optional ,intent(in) :: numEns - +#endif !EOP ! pio variables are already in SDAT no need to copy them +#ifndef USE_PDAF 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, caseId=caseId, dt=dt, numEns=numEns, nzg=nzg, domZvarName=domZvarName,& + filePath, filename, fldListFile, fldListModel, 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) - +#else + 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, & + 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) +#endif end subroutine shr_strdata_create_oldway !=============================================================================== diff --git a/src/csm_share/streams/shr_stream_mod.F90 b/src/csm_share/streams/shr_stream_mod.F90 index 1717efbebb..4fc371f3a9 100644 --- a/src/csm_share/streams/shr_stream_mod.F90 +++ b/src/csm_share/streams/shr_stream_mod.F90 @@ -77,9 +77,10 @@ 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 - +#ifdef USE_PDAF ! Yorck public :: shr_stream_get_dt_and_caseId ! return dt and caseId +#endif ! public :: shr_stream_bcast ! broadcast a stream (untested) ! !PUBLIC DATA MEMBERS: @@ -154,13 +155,15 @@ 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 +#ifdef USE_PDAF + ! 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 +#endif character(SHR_KIND_CS) :: tInterpAlgo ! Algorithm to use for time interpolation character(SHR_KIND_CL) :: calendar ! stream calendar @@ -653,6 +656,7 @@ subroutine shr_stream_init(strm,infoFile,yearFirst,yearLast,yearAlign,taxMode,rc close(nUnit) +#ifdef USE_PDAF !----------------------------------------------------------------------------- ! Yorck adaptations: read in time resolution of forcings and caseId, numEns !----------------------------------------------------------------------------- @@ -722,6 +726,7 @@ subroutine shr_stream_init(strm,infoFile,yearFirst,yearLast,yearAlign,taxMode,rc end if close(nUnit) +#endif !----------------------------------------------------------------------------- ! get initial calendar value @@ -762,12 +767,19 @@ end subroutine shr_stream_init ! ! !INTERFACE: ------------------------------------------------------------------ +#ifdef USE_PDAF 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) - +#else + subroutine shr_stream_set(strm,yearFirst,yearLast,yearAlign,offset,taxMode, & + fldListFile,fldListModel,domFilePath,domFileName, & + domTvarName,domXvarName,domYvarName,domZvarName, & + domAreaName,domMaskName, & + filePath,filename,dataSource,rc) +#endif ! !INPUT/OUTPUT PARAMETERS: type(shr_stream_streamType) ,intent(inout) :: strm ! data stream @@ -791,11 +803,12 @@ subroutine shr_stream_set(strm,yearFirst,yearLast,yearAlign,offset, caseId, dt, character(*) ,optional,intent(in) :: dataSource ! comment line integer (SHR_KIND_IN),optional,intent(out) :: rc ! return code - ! Yorck +#ifdef USE_PDAF + ! 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 - +#endif !EOP !----- local ----- @@ -877,6 +890,7 @@ subroutine shr_stream_set(strm,yearFirst,yearLast,yearAlign,offset, caseId, dt, call fileVec%move_out(strm%file) endif +#ifdef USE_PDAF ! Yorck if (present(caseId)) then strm%caseId = caseId @@ -889,6 +903,7 @@ subroutine shr_stream_set(strm,yearFirst,yearLast,yearAlign,offset, caseId, dt, if (present(numEns)) then strm%numEns = numEns end if +#endif !----------------------------------------------------------------------------- ! get initial calendar value @@ -970,9 +985,11 @@ subroutine shr_stream_default(strm,rc) strm%domAreaName = ' ' strm%domMaskName = ' ' +#ifdef USE_PDAF strm%caseId = 0 strm%dt = 0 strm%numEns = 0 +#endif strm%calendar = shr_cal_noleap @@ -1619,9 +1636,11 @@ 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(:) +#ifdef USE_PDAF character(SHR_KIND_CL) :: mfldName !Yorck integer(SHR_KIND_IN) :: nt_ ! Yorck real(SHR_KIND_R8),allocatable :: tvar_lok(:) ! Yorck +#endif !----- formats ----- character(*),parameter :: subname = '(shr_stream_readTCoord) ' character(*),parameter :: F01 = "('(shr_stream_readTCoord) ',a,2i7)" @@ -1645,6 +1664,7 @@ subroutine shr_stream_readTCoord(strm,k,rc) if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_inquire_dimension') deallocate(dids) +#ifdef USE_PDAF ! Yorck ! 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 @@ -1664,6 +1684,7 @@ subroutine shr_stream_readTCoord(strm,k,rc) ! (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 +#endif allocate(strm%file(k)%date(nt),strm%file(k)%secs(nt)) strm%file(k)%nt = nt @@ -1689,7 +1710,7 @@ subroutine shr_stream_readTCoord(strm,k,rc) strm%calendar = trim(shr_cal_calendarName(trim(calendar))) allocate(tvar(nt)) - +#ifdef USE_PDAF ! Yorck if (INDEX(mfldName, "noise") .ne. 0) then allocate(tvar_lok(nt_)) @@ -1704,12 +1725,15 @@ subroutine shr_stream_readTCoord(strm,k,rc) 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') +#endif + 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') +#ifdef USE_PDAF end if ! End Yorck +#endif do n = 1,nt call shr_cal_advDate(tvar(n),bunits,bdate,bsec,ndate,nsec,calendar) strm%file(k)%date(n) = ndate @@ -2725,10 +2749,11 @@ 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 +#ifdef USE_PDAF write(nUnit) strm(k)%dt ! Yorck write(nUnit) strm(k)%caseId ! Yorck write(nUnit) strm(k)%numEns ! Yorck - +#endif end do close(nUnit) @@ -3032,9 +3057,11 @@ subroutine shr_stream_dataDump(strm) write(s_logunit,F00) "domAreaName = ", trim(strm%domAreaName) write(s_logunit,F00) "domMaskName = ", trim(strm%domMaskName) +#ifdef USE_PDAF write(s_logunit,F01) "dt = ", strm%dt ! Yorck write(s_logunit,F01) "caseId = ", strm%caseId ! Yorck write(s_logunit,F01) "numEns = ", strm%numEns !Yorck +#endif end subroutine shr_stream_dataDump @@ -3384,9 +3411,11 @@ subroutine shr_stream_bcast(stream,comm,rc) call shr_mpi_bcast(stream%domMaskName ,comm,subName) call shr_mpi_bcast(stream%calendar ,comm,subName) +#ifdef USE_PDAF !Yorck call shr_mpi_bcast(stream%dt ,comm,subName) call shr_mpi_bcast(stream%caseId ,comm,subName) +#endif if (pid /= 0) allocate(stream%file(stream%nFiles)) @@ -3402,7 +3431,7 @@ subroutine shr_stream_bcast(stream,comm,rc) end subroutine shr_stream_bcast - +#ifdef USE_PDAF ! Yorck subroutine shr_stream_get_dt_and_caseId(stream,dt,caseId) @@ -3425,6 +3454,7 @@ subroutine shr_stream_get_dt_and_caseId(stream,dt,caseId) end subroutine shr_stream_get_dt_and_caseId +#endif !=============================================================================== end module shr_stream_mod diff --git a/src/datm/datm_comp_mod.F90 b/src/datm/datm_comp_mod.F90 index e68a141abe..f7393d0c9e 100644 --- a/src/datm/datm_comp_mod.F90 +++ b/src/datm/datm_comp_mod.F90 @@ -77,7 +77,9 @@ 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 +#ifdef USE_PDAF integer(IN) :: stbot_noise, sprecn_noise, sswdn_noise, slwdn_noise +#endif integer(IN) :: sprecc,sprecl,sprecn,sco2p,sco2d,sswup,sprec,starcf ! water isotopes / tracer input @@ -162,7 +164,11 @@ 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. +#ifdef USE_PDAF integer(IN),parameter :: ktranss = 37 +#else + integer(IN),parameter :: ktranss = 33 +#endif character(16),parameter :: stofld(1:ktranss) = & (/"strm_tbot ","strm_wind ","strm_z ","strm_pbot ", & @@ -175,9 +181,11 @@ 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 ", & - ! add perturbations (Yorck) - "strm_tbot_noise ","strm_precn_noise","strm_swdn_noise ","strm_lwdn_noise " & + "strm_precn_16O ","strm_precn_18O ","strm_precn_HDO " & +#ifdef USE_PDAF + ! add perturbations (Yorck) + ,"strm_tbot_noise ","strm_precn_noise","strm_swdn_noise ","strm_lwdn_noise " & +#endif /) character(16),parameter :: stifld(1:ktranss) = & @@ -192,9 +200,11 @@ 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 " & +#ifdef USE_PDAF ! add perturbations (Yorck) - "tbot_noise ","precn_noise ","swdn_noise ","lwdn_noise " & + ,"tbot_noise ","precn_noise ","swdn_noise ","lwdn_noise " & +#endif /) character(CL), pointer :: ilist_av(:) ! input list for translation @@ -449,7 +459,7 @@ 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') - +#ifdef USE_PDAF ! 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. @@ -460,9 +470,7 @@ subroutine datm_comp_init(Eclock, x2a, a2x, & 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') - - - +#endif ! anomaly forcing sprecsf = mct_aVect_indexRA(avstrm,'strm_precsf' ,perrWith='quiet') @@ -920,14 +928,15 @@ 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 +#ifdef USE_PDAF 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 + if (stbot_noise > 1 ) then + write(logunit,*) trim(subname),' Using noise from files ' + end if end if - - +#endif endif lsize = mct_avect_lsize(a2x) do n = 1,lsize @@ -936,8 +945,7 @@ subroutine datm_comp_run(EClock, x2a, a2x, & !--- temperature --- if (tbotmax < 50.0_R8) a2x%rAttr(ktbot,n) = a2x%rAttr(ktbot,n) + tkFrz - - +#ifdef USE_PDAF !!! begin perturbation block (Yorck) --> perturb temperature, long wave radiation, short wave radiation and precipitation ! --> further variables can be introduced (changes also have to be made in reading in routines of streams for this, for now only ! these variables are included) @@ -963,7 +971,7 @@ subroutine datm_comp_run(EClock, x2a, a2x, & end if !!! end perturbation block - +#endif ! Limit very cold forcing to 180K a2x%rAttr(ktbot,n) = max(180._r8, a2x%rAttr(ktbot,n)) a2x%rAttr(kptem,n) = a2x%rAttr(ktbot,n) @@ -1094,6 +1102,7 @@ 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 From 559494a83a2c58386f3d96c51704b80a5ee8a351 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Mon, 10 Nov 2025 13:32:07 +0100 Subject: [PATCH 06/21] indentation --- src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 index 7806bd3c38..082ffa08cc 100644 --- a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 @@ -314,9 +314,9 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) if (.not. readvar) then parameters_in_file = .False. end if - end if + end if - call ncd_io(ncid=ncid, varname='THETAS_adj', flag='read', data=thetas_adj, dim1name=grlnd, readvar=readvar) + call ncd_io(ncid=ncid, varname='THETAS_adj', flag='read', data=thetas_adj, dim1name=grlnd, readvar=readvar) if (.not. readvar) then parameters_in_file_adj = .False. else @@ -339,7 +339,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) if (.not. readvar) then parameters_in_file_adj = .False. end if - end if + end if #endif do p = begp,endp g = patch%gridcell(p) From 8eaa7b1c499dcb6071bae67fb1ed49f9e2621967 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Mon, 10 Nov 2025 13:41:11 +0100 Subject: [PATCH 07/21] whitespace --- src/datm/datm_comp_mod.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/datm/datm_comp_mod.F90 b/src/datm/datm_comp_mod.F90 index f7393d0c9e..0563ea5366 100644 --- a/src/datm/datm_comp_mod.F90 +++ b/src/datm/datm_comp_mod.F90 @@ -201,10 +201,10 @@ module datm_comp_mod ! isotopic forcing "rh_18O ","rh_HDO ", & "precn_16O ","precn_18O ","precn_HDO " & -#ifdef USE_PDAF +#ifdef USE_PDAF ! add perturbations (Yorck) ,"tbot_noise ","precn_noise ","swdn_noise ","lwdn_noise " & -#endif +#endif /) character(CL), pointer :: ilist_av(:) ! input list for translation From a2c0c2f50fd4a1be666118ef40ae3b4d0de49093 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Mon, 10 Nov 2025 13:41:21 +0100 Subject: [PATCH 08/21] replace hard-coded 10 by nlevsoifl --- src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 index 082ffa08cc..f623a1722e 100644 --- a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 @@ -545,11 +545,11 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) soilstate_inst%sucsat_col(c,lev) = psis_sat(col%gridcell(c), lev) xksat = ks(col%gridcell(c), lev) ! mm/s else - ! Use the value from the 10th soil level as a default value - soilstate_inst%watsat_col(c,lev) = thetas(col%gridcell(c), 10) - soilstate_inst%bsw_col(c,lev) = shape_param(col%gridcell(c), 10) - soilstate_inst%sucsat_col(c,lev) = psis_sat(col%gridcell(c), 10) - xksat = ks(col%gridcell(c), 10) + ! Use the value from the 10th (default of nlevsoifl) soil level as a default value + soilstate_inst%watsat_col(c,lev) = thetas(col%gridcell(c), nlevsoifl) + soilstate_inst%bsw_col(c,lev) = shape_param(col%gridcell(c), nlevsoifl) + soilstate_inst%sucsat_col(c,lev) = psis_sat(col%gridcell(c), nlevsoifl) + xksat = ks(col%gridcell(c), nlevsoifl) end if end if #endif From 10225087e856251e0f2427b50b5e749a0c0dc843 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Mon, 10 Nov 2025 14:11:02 +0100 Subject: [PATCH 09/21] docs: Add perturbation documentation uses AI --- README.md | 4 - docs/_toc.yml | 2 + docs/users_guide/introduction/perturbation.md | 74 +++++++++++++++++++ 3 files changed, 76 insertions(+), 4 deletions(-) create mode 100644 docs/users_guide/introduction/perturbation.md diff --git a/README.md b/README.md index 2e6e9bf7a1..3af2d8e249 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,3 @@ -`dev-perturbation`-Branch: Generated for usage of eCLM-PDAF with -Perturbation routine from -https://gitlab.jsc.fz-juelich.de/detect/cluster-c/c01/perturbationroutineclm5 - # eCLM [![CI](https://github.com/HPSCTerrSys/eCLM/actions/workflows/build_and_run_eCLM.yml/badge.svg)](https://github.com/HPSCTerrSys/eCLM/actions/workflows/build_and_run_eCLM.yml) diff --git a/docs/_toc.yml b/docs/_toc.yml index 7f74453e61..659415e6bb 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -9,6 +9,8 @@ parts: title: Installing eCLM - file: users_guide/introduction/introduction title: Scientific Background + - file: users_guide/introduction/perturbation + title: Perturbation routines - caption: User's Guide chapters: diff --git a/docs/users_guide/introduction/perturbation.md b/docs/users_guide/introduction/perturbation.md new file mode 100644 index 0000000000..95872d47ba --- /dev/null +++ b/docs/users_guide/introduction/perturbation.md @@ -0,0 +1,74 @@ +## Perturbation routines + +This section describes **perturbation capabilities for eCLM-PDAF**. + +The implementation of the perturbation routines was first introduced +by Yorck Ewerdwalbesloh in +https://gitlab.jsc.fz-juelich.de/detect/cluster-c/c01/perturbationroutineclm5. + +### Key Components + +#### 1. Soil Parameter Perturbation (`SoilStateInitTimeConstMod.F90`) + +Allows reading perturbed hydraulic parameters from input files instead of computing them via pedotransfer functions. + +**Parameters:** +- `THETAS` - saturated water content +- `SHAPE_PARAM` - Brooks-Corey shape parameter (`bsw`) +- `PSIS_SAT` - saturated matric potential (`sucsat`) +- `KSAT` - saturated hydraulic conductivity (`xksat`) + +**Implementation:** +- Reads both original parameters (for `nlevsoifl=10` soil layers) and adjusted parameters with `_adj` suffix (for all `nlevgrnd` layers) +- Falls back to pedotransfer functions if parameters aren't in the file +- Modifies organic matter mixing to preserve perturbed parameter values + +#### 2. Noise-Based Forcing Perturbation + +Adds spatiotemporal noise to atmospheric forcing data for ensemble data assimilation. + +**Modified files:** +- `shr_stream_mod.F90` +- `shr_strdata_mod.F90` +- `shr_dmodel_mod.F90` + +**Key Features:** +- Stores ensemble metadata in stream structure: + - `caseId` - ensemble member ID + - `dt` - forcing time resolution + - `numEns` - total ensemble size +- **Time-shifting mechanism**: Different ensemble members read different temporal frames from the same noise file + - Formula: `frame = nt + caseId * 24/dt` + - Example: For 3-hourly data (`dt=3`), member 0 starts at frame `nt`, member 1 at `nt+8`, member 2 at `nt+16`, etc. +- Detects noise fields by checking if model field name contains `"noise"` +- Adjusts time coordinate reading to account for ensemble-extended noise files + +#### 3. DATM Integration (`datm_comp_mod.F90`) + +Passes ensemble information (`inst_index`, `dt_option`, `ninst`) from +the data atmosphere (DATM) component to the stream infrastructure, +enabling the noise perturbation mechanism to operate correctly within +CESM's multi-instance framework. + +### Design Rationale + +1. **Dual perturbation approach**: + - Parameter perturbation for soil properties (offline preprocessing) + - Forcing perturbation via temporal noise patterns (runtime) + +2. **Backward compatibility**: All changes are guarded by `USE_PDAF` preprocessor directives + +3. **Efficiency**: Uses a single noise file for all ensemble members by time-shifting, avoiding storage duplication + +4. **Flexibility**: Soil parameters can be perturbed or computed traditionally based on file contents + +### Modified Files + +- `src/clm5/biogeophys/SoilStateInitTimeConstMod.F90` - 121 additions +- `src/csm_share/streams/shr_stream_mod.F90` - 211 additions +- `src/csm_share/streams/shr_strdata_mod.F90` - 72 additions +- `src/csm_share/streams/shr_dmodel_mod.F90` - 40 additions +- `src/datm/datm_comp_mod.F90` - 62 additions +- `README.md` - documentation + + From c488fc421b3de23b60df8d7c80b7e9d9574d85c5 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Wed, 12 Nov 2025 17:47:12 +0100 Subject: [PATCH 10/21] docs: general info to docs, specifics into source code --- docs/users_guide/introduction/perturbation.md | 10 ++++++- src/datm/datm_comp_mod.F90 | 30 +++++++++++-------- 2 files changed, 27 insertions(+), 13 deletions(-) diff --git a/docs/users_guide/introduction/perturbation.md b/docs/users_guide/introduction/perturbation.md index 95872d47ba..faca86eaaf 100644 --- a/docs/users_guide/introduction/perturbation.md +++ b/docs/users_guide/introduction/perturbation.md @@ -25,7 +25,15 @@ Allows reading perturbed hydraulic parameters from input files instead of comput #### 2. Noise-Based Forcing Perturbation -Adds spatiotemporal noise to atmospheric forcing data for ensemble data assimilation. +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. + +Adds spatiotemporal noise to atmospheric forcing data for ensemble +data assimilation. **Modified files:** - `shr_stream_mod.F90` diff --git a/src/datm/datm_comp_mod.F90 b/src/datm/datm_comp_mod.F90 index 0563ea5366..25f4b6b018 100644 --- a/src/datm/datm_comp_mod.F90 +++ b/src/datm/datm_comp_mod.F90 @@ -78,7 +78,10 @@ module datm_comp_mod integer(IN) :: kanidr,kanidf,kavsdr,kavsdf integer(IN) :: stbot,swind,sz,spbot,sshum,stdew,srh,slwdn,sswdn,sswdndf,sswdndr #ifdef USE_PDAF - integer(IN) :: stbot_noise, sprecn_noise, sswdn_noise, slwdn_noise + integer(IN) :: stbot_noise ! temperature noise + integer(IN) :: sprecn_noise ! precipitation noise + integer(IN) :: sswdn_noise ! shortwave radiation noise + integer(IN) :: slwdn_noise ! longwave radiation noise #endif integer(IN) :: sprecc,sprecl,sprecn,sco2p,sco2d,sswup,sprec,starcf @@ -183,7 +186,7 @@ module datm_comp_mod "strm_rh_18O ","strm_rh_HDO ", & "strm_precn_16O ","strm_precn_18O ","strm_precn_HDO " & #ifdef USE_PDAF - ! add perturbations (Yorck) + ! add DATM-internal stream variable names for noise fields (Yorck) ,"strm_tbot_noise ","strm_precn_noise","strm_swdn_noise ","strm_lwdn_noise " & #endif /) @@ -202,7 +205,7 @@ module datm_comp_mod "rh_18O ","rh_HDO ", & "precn_16O ","precn_18O ","precn_HDO " & #ifdef USE_PDAF - ! add perturbations (Yorck) + ! add NetCDF input file variable names for noise fields (Yorck) ,"tbot_noise ","precn_noise ","swdn_noise ","lwdn_noise " & #endif /) @@ -460,12 +463,7 @@ subroutine datm_comp_init(Eclock, x2a, a2x, & starcf = mct_aVect_indexRA(avstrm,'strm_tarcf' ,perrWith='quiet') #ifdef USE_PDAF - ! 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) - + ! Get Noise Field Indices from Stream (noise streams are optional) 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') @@ -931,6 +929,7 @@ subroutine datm_comp_run(EClock, x2a, a2x, & if (my_task == master_task) & write(logunit,*) trim(subname),' max values = ',tbotmax,tdewmax,anidrmax #ifdef USE_PDAF + ! Log if temperature noise read from stream file if (my_task == master_task) then if (stbot_noise > 1 ) then write(logunit,*) trim(subname),' Using noise from files ' @@ -946,27 +945,34 @@ subroutine datm_comp_run(EClock, x2a, a2x, & !--- temperature --- if (tbotmax < 50.0_R8) a2x%rAttr(ktbot,n) = a2x%rAttr(ktbot,n) + tkFrz #ifdef USE_PDAF - !!! begin perturbation block (Yorck) --> perturb temperature, long wave radiation, short wave radiation and precipitation - ! --> further variables can be introduced (changes also have to be made in reading in routines of streams for this, for now only - ! these variables are included) + !!! begin perturbation block (Yorck) --> perturb + !!! temperature, long wave radiation, short wave radiation + !!! and precipitation + ! --> further variables can be introduced (changes also have + ! to be made in reading in routines of streams for this, for + ! now only these variables are included) ! tbot if (stbot_noise > 0) then + ! Additive noise a2x%rAttr(ktbot,n) = a2x%rAttr(ktbot,n) + avstrm%rAttr(stbot_noise,n) end if ! lwdn if (slwdn_noise > 0) then + ! Additive noise avstrm%rAttr(slwdn,n) = avstrm%rAttr(slwdn,n) + avstrm%rAttr(slwdn_noise,n) end if ! swdn if (sswdn_noise>0) then + ! Multiplicative noise avstrm%rAttr(sswdn,n) = avstrm%rAttr(sswdn,n) * avstrm%rAttr(sswdn_noise,n) end if ! sprecn if (sprecn_noise > 0) then + ! Multiplicative noise avstrm%rAttr(sprecn,n) = avstrm%rAttr(sprecn,n)*avstrm%rAttr(sprecn_noise,n) end if From 3e52efdfb6c8fc7f757ecd0cec837de52a2f1209 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Wed, 12 Nov 2025 20:33:23 +0100 Subject: [PATCH 11/21] bugfix wrong ifdef --- src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 index f623a1722e..a3cc6ece47 100644 --- a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 @@ -249,7 +249,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) allocate(shape_param_adj(begg:endg,nlevgrnd)) allocate(psis_sat_adj(begg:endg,nlevgrnd)) allocate(ks_adj(begg:endg,nlevgrnd)) -#ifdef +#endif ! Determine organic_max from parameter file From b80fe6d63942e22c4bc18ceb8ee5511b3c254e45 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 13 Nov 2025 10:36:53 +0100 Subject: [PATCH 12/21] Document the status of organic matter mixing in particular the difference in setting the Brooks-Corey Shape Parameter --- docs/users_guide/introduction/perturbation.md | 37 ++++++++++++++++++- 1 file changed, 36 insertions(+), 1 deletion(-) diff --git a/docs/users_guide/introduction/perturbation.md b/docs/users_guide/introduction/perturbation.md index faca86eaaf..82805639b8 100644 --- a/docs/users_guide/introduction/perturbation.md +++ b/docs/users_guide/introduction/perturbation.md @@ -19,10 +19,45 @@ Allows reading perturbed hydraulic parameters from input files instead of comput - `KSAT` - saturated hydraulic conductivity (`xksat`) **Implementation:** -- Reads both original parameters (for `nlevsoifl=10` soil layers) and adjusted parameters with `_adj` suffix (for all `nlevgrnd` layers) +- Reads both original parameters (for `nlevsoifl=10` soil layers) and + applies organic matter mixing +- OR: Reads adjusted parameters with `_adj` suffix (for all `nlevgrnd` + layers) and **overwrites** parameters from organic matter mixing. - Falls back to pedotransfer functions if parameters aren't in the file - Modifies organic matter mixing to preserve perturbed parameter values +##### Brooks-Corey Shape Parameter Mixing with Organic Matter + +**Backward Compatibility Note:** In the branch `dev-perturbation` it +is modified how the Brooks-Corey shape parameter (`bsw`) is mixed with +organic matter, which causes numerical differences compared to the +standard CLM configuration. + +**Standard CLM behavior** +```fortran +bsw = (1-om_frac) * (2.91 + 0.159*clay) + om_frac*om_b +``` + +`dev-perturbation` **behavior** +```fortran +bsw = (1-om_frac) * bsw + om_frac*om_b +``` + +**Impact:** In standard CLM, the mineral soil contribution to `bsw` in +organic matter mixing always uses the hard-coded Cosby et al. (1984) +Table 5 formula (`2.91 + 0.159*clay`), regardless of which +pedotransfer function was used to initially compute `bsw`. In +`dev-perturbation`, the mixing uses the actual `bsw` value computed by +the selected pedotransfer function (or read from file if parameters +are provided). + +This change affects both soil and lake columns (see +`SoilStateInitTimeConstMod.F90:565,567,675,677`). + +When comparing test cases to numerical precision, this will produce +differences even when not using perturbed parameters from input files, +since the organic matter mixing calculation differs. + #### 2. Noise-Based Forcing Perturbation When having an ensemble run, the memory consumption is too large if From e1eee04d4d09b952bf5e50fcca07f9b7583f287e Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 13 Nov 2025 11:08:58 +0100 Subject: [PATCH 13/21] backward compatible shape parameter handling --- docs/users_guide/introduction/perturbation.md | 36 +++---------------- .../biogeophys/SoilStateInitTimeConstMod.F90 | 12 +++++-- 2 files changed, 15 insertions(+), 33 deletions(-) diff --git a/docs/users_guide/introduction/perturbation.md b/docs/users_guide/introduction/perturbation.md index 82805639b8..a32ce9759e 100644 --- a/docs/users_guide/introduction/perturbation.md +++ b/docs/users_guide/introduction/perturbation.md @@ -26,37 +26,11 @@ Allows reading perturbed hydraulic parameters from input files instead of comput - Falls back to pedotransfer functions if parameters aren't in the file - Modifies organic matter mixing to preserve perturbed parameter values -##### Brooks-Corey Shape Parameter Mixing with Organic Matter - -**Backward Compatibility Note:** In the branch `dev-perturbation` it -is modified how the Brooks-Corey shape parameter (`bsw`) is mixed with -organic matter, which causes numerical differences compared to the -standard CLM configuration. - -**Standard CLM behavior** -```fortran -bsw = (1-om_frac) * (2.91 + 0.159*clay) + om_frac*om_b -``` - -`dev-perturbation` **behavior** -```fortran -bsw = (1-om_frac) * bsw + om_frac*om_b -``` - -**Impact:** In standard CLM, the mineral soil contribution to `bsw` in -organic matter mixing always uses the hard-coded Cosby et al. (1984) -Table 5 formula (`2.91 + 0.159*clay`), regardless of which -pedotransfer function was used to initially compute `bsw`. In -`dev-perturbation`, the mixing uses the actual `bsw` value computed by -the selected pedotransfer function (or read from file if parameters -are provided). - -This change affects both soil and lake columns (see -`SoilStateInitTimeConstMod.F90:565,567,675,677`). - -When comparing test cases to numerical precision, this will produce -differences even when not using perturbed parameters from input files, -since the organic matter mixing calculation differs. +##### Note about the Brooks-Corey Shape Parameter + +When perturbed soil parameters are read from input files, the organic +matter mixing for `bsw` uses the file-read value instead of the +hard-coded Cosby et al. (1984) Table 5 formula (`2.91 + 0.159*clay`). #### 2. Noise-Based Forcing Perturbation diff --git a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 index a3cc6ece47..b3054d47dd 100644 --- a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 @@ -564,7 +564,11 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) #ifndef USE_PDAF soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * (2.91_r8 + 0.159_r8*clay) + om_frac*om_b #else - soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * soilstate_inst%bsw_col(c,lev) + om_frac*om_b + if (parameters_in_file) then + soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * soilstate_inst%bsw_col(c,lev) + om_frac*om_b + else + soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac) * (2.91_r8 + 0.159_r8*clay) + om_frac*om_b + end if #endif soilstate_inst%sucsat_col(c,lev) = (1._r8-om_frac) * soilstate_inst%sucsat_col(c,lev) + om_sucsat*om_frac soilstate_inst%hksat_min_col(c,lev) = xksat @@ -674,7 +678,11 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) #ifndef USE_PDAF soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*(2.91_r8 + 0.159_r8*clay) + om_frac * om_b_lake #else - soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*soilstate_inst%bsw_col(c,lev) + om_frac * om_b_lake + if (parameters_in_file .or. parameters_in_file_adj) then + soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*soilstate_inst%bsw_col(c,lev) + om_frac * om_b_lake + else + soilstate_inst%bsw_col(c,lev) = (1._r8-om_frac)*(2.91_r8 + 0.159_r8*clay) + om_frac * om_b_lake + end if #endif soilstate_inst%sucsat_col(c,lev) = (1._r8-om_frac)*soilstate_inst%sucsat_col(c,lev) + om_sucsat_lake * om_frac From ff88af2c2e85c03b231272bdae1af762fae3dbc0 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 13 Nov 2025 13:00:04 +0100 Subject: [PATCH 14/21] docs: clarify `numEns` after check from Yorck Ewerdwalbesloh --- docs/users_guide/introduction/perturbation.md | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/docs/users_guide/introduction/perturbation.md b/docs/users_guide/introduction/perturbation.md index a32ce9759e..c43970d267 100644 --- a/docs/users_guide/introduction/perturbation.md +++ b/docs/users_guide/introduction/perturbation.md @@ -53,10 +53,15 @@ data assimilation. - Stores ensemble metadata in stream structure: - `caseId` - ensemble member ID - `dt` - forcing time resolution - - `numEns` - total ensemble size + - `numEns` - ensemble sizethat the perturbation file was created + for. Example: Running with 50 ensemble members with a noise file + created for an ensemble of up to 200 members. Then, `numEns` + should be set to 200 in the stream file for right usage of + perturbation dimension. - **Time-shifting mechanism**: Different ensemble members read different temporal frames from the same noise file - Formula: `frame = nt + caseId * 24/dt` - - Example: For 3-hourly data (`dt=3`), member 0 starts at frame `nt`, member 1 at `nt+8`, member 2 at `nt+16`, etc. + - Example: For 3-hourly data (`dt=3`), member 0 starts at frame + `nt`, member 1 at `nt+8`, member 2 at `nt+16`, etc. - Detects noise fields by checking if model field name contains `"noise"` - Adjusts time coordinate reading to account for ensemble-extended noise files From 2347dcffdbc2211fb19afcab1b757158b1a01a0d Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 13 Nov 2025 13:22:23 +0100 Subject: [PATCH 15/21] docs: fix and updates --- docs/users_guide/introduction/perturbation.md | 29 +++++++++---------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/docs/users_guide/introduction/perturbation.md b/docs/users_guide/introduction/perturbation.md index c43970d267..2a747327c8 100644 --- a/docs/users_guide/introduction/perturbation.md +++ b/docs/users_guide/introduction/perturbation.md @@ -1,4 +1,4 @@ -## Perturbation routines +# Perturbation routines # This section describes **perturbation capabilities for eCLM-PDAF**. @@ -6,9 +6,9 @@ The implementation of the perturbation routines was first introduced by Yorck Ewerdwalbesloh in https://gitlab.jsc.fz-juelich.de/detect/cluster-c/c01/perturbationroutineclm5. -### Key Components +## Key Components ## -#### 1. Soil Parameter Perturbation (`SoilStateInitTimeConstMod.F90`) +### 1. Soil Parameter Perturbation (`SoilStateInitTimeConstMod.F90`) ### Allows reading perturbed hydraulic parameters from input files instead of computing them via pedotransfer functions. @@ -26,13 +26,13 @@ Allows reading perturbed hydraulic parameters from input files instead of comput - Falls back to pedotransfer functions if parameters aren't in the file - Modifies organic matter mixing to preserve perturbed parameter values -##### Note about the Brooks-Corey Shape Parameter +#### Note about the Brooks-Corey Shape Parameter #### When perturbed soil parameters are read from input files, the organic matter mixing for `bsw` uses the file-read value instead of the hard-coded Cosby et al. (1984) Table 5 formula (`2.91 + 0.159*clay`). -#### 2. Noise-Based Forcing Perturbation +### 2. Noise-Based Forcing Perturbation ### 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 @@ -53,7 +53,7 @@ data assimilation. - Stores ensemble metadata in stream structure: - `caseId` - ensemble member ID - `dt` - forcing time resolution - - `numEns` - ensemble sizethat the perturbation file was created + - `numEns` - ensemble size that the perturbation file was created for. Example: Running with 50 ensemble members with a noise file created for an ensemble of up to 200 members. Then, `numEns` should be set to 200 in the stream file for right usage of @@ -65,14 +65,14 @@ data assimilation. - Detects noise fields by checking if model field name contains `"noise"` - Adjusts time coordinate reading to account for ensemble-extended noise files -#### 3. DATM Integration (`datm_comp_mod.F90`) +### 3. DATM Integration (`datm_comp_mod.F90`) ### Passes ensemble information (`inst_index`, `dt_option`, `ninst`) from the data atmosphere (DATM) component to the stream infrastructure, enabling the noise perturbation mechanism to operate correctly within CESM's multi-instance framework. -### Design Rationale +## Design Rationale ## 1. **Dual perturbation approach**: - Parameter perturbation for soil properties (offline preprocessing) @@ -84,13 +84,12 @@ CESM's multi-instance framework. 4. **Flexibility**: Soil parameters can be perturbed or computed traditionally based on file contents -### Modified Files +## Modified Source Code ## -- `src/clm5/biogeophys/SoilStateInitTimeConstMod.F90` - 121 additions -- `src/csm_share/streams/shr_stream_mod.F90` - 211 additions -- `src/csm_share/streams/shr_strdata_mod.F90` - 72 additions -- `src/csm_share/streams/shr_dmodel_mod.F90` - 40 additions -- `src/datm/datm_comp_mod.F90` - 62 additions -- `README.md` - documentation +- `src/clm5/biogeophys/SoilStateInitTimeConstMod.F90` +- `src/csm_share/streams/shr_stream_mod.F90` +- `src/csm_share/streams/shr_strdata_mod.F90` +- `src/csm_share/streams/shr_dmodel_mod.F90` +- `src/datm/datm_comp_mod.F90` From 39745b5d5ff36a4db88087915a5f7e070dc1ec5d Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 13 Nov 2025 14:00:46 +0100 Subject: [PATCH 16/21] GRACE-DA: Everything under USE_PDAF --- src/clm5/biogeophys/SoilHydrologyType.F90 | 6 ++++-- src/clm5/biogeophys/WaterStateType.F90 | 10 +++++---- src/clm5/cpl/lnd_comp_mct.F90 | 2 +- src/clm5/main/atm2lndType.F90 | 6 ++++-- src/clm5/main/clm_driver.F90 | 5 ++++- src/clm5/main/clm_varcon.F90 | 4 ++-- src/clm5/main/lnd2atmMod.F90 | 26 +++++++++++++++++++---- 7 files changed, 43 insertions(+), 16 deletions(-) diff --git a/src/clm5/biogeophys/SoilHydrologyType.F90 b/src/clm5/biogeophys/SoilHydrologyType.F90 index af6de9165e..c62ec465f1 100644 --- a/src/clm5/biogeophys/SoilHydrologyType.F90 +++ b/src/clm5/biogeophys/SoilHydrologyType.F90 @@ -53,10 +53,11 @@ Module SoilHydrologyType real(r8), pointer :: i_0_col (:) ! col VIC average saturation in top soil layers real(r8), pointer :: ice_col (:,:) ! col VIC soil ice (kg/m2) for VIC soil layers +#ifdef USE_PDAF ! Yorck real(r8), pointer :: wa_col_inc (:) ! increment col water in the unconfined aquifer (mm) real(r8), pointer :: wa_col_mean (:) ! mean col water in the unconfined aquifer (mm) - +#endif contains ! Public routines @@ -121,10 +122,11 @@ subroutine InitAllocate(this, bounds) allocate(this%zwts_col (begc:endc)) ; this%zwts_col (:) = nan allocate(this%wa_col (begc:endc)) ; this%wa_col (:) = nan +#ifdef USE_PDAF ! Yorck allocate(this%wa_col_inc (begc:endc)) ; this%wa_col_inc (:) = nan allocate(this%wa_col_mean (begc:endc)) ; this%wa_col_mean (:) = nan - +#endif allocate(this%qcharge_col (begc:endc)) ; this%qcharge_col (:) = nan allocate(this%fracice_col (begc:endc,nlevgrnd)) ; this%fracice_col (:,:) = nan allocate(this%icefrac_col (begc:endc,nlevgrnd)) ; this%icefrac_col (:,:) = nan diff --git a/src/clm5/biogeophys/WaterStateType.F90 b/src/clm5/biogeophys/WaterStateType.F90 index 2cd4a4a715..fdec7cf1b3 100644 --- a/src/clm5/biogeophys/WaterStateType.F90 +++ b/src/clm5/biogeophys/WaterStateType.F90 @@ -54,7 +54,7 @@ module WaterstateType real(r8), pointer :: ice1_grc (:) ! grc initial gridcell total h2o ice content real(r8), pointer :: ice2_grc (:) ! grc post land cover change total ice content real(r8), pointer :: tws_grc (:) ! grc total water storage (mm H2O) - +#ifdef USE_PDAF ! Yorck additions, variables for getting monthly means which can be compared to a GRACE measurements, other variables will just ! provide instantaneous values @@ -104,7 +104,7 @@ module WaterstateType real(r8), pointer :: h2ocan_state_after (:) ! canopy state ! END Yorck - +#endif #ifdef COUP_OAS_PFL real(r8), pointer :: pfl_psi_col (:,:) ! ParFlow pressure head COUP_OAS_PFL real(r8), pointer :: pfl_h2osoi_liq_col (:,:) ! ParFlow soil liquid COUP_OAS_PFL @@ -262,6 +262,7 @@ subroutine InitAllocate(this, bounds) allocate(this%ice2_grc (begg:endg)) ; this%ice2_grc (:) = nan allocate(this%tws_grc (begg:endg)) ; this%tws_grc (:) = nan +#ifdef USE_PDAF ! Yorck additions (see above) allocate(this%h2osoi_ice_col_mean (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_ice_col_mean (:,:) = nan allocate(this%h2osoi_liq_col_mean (begc:endc,-nlevsno+1:nlevgrnd)) ; this%h2osoi_liq_col_mean (:,:) = nan @@ -296,7 +297,7 @@ subroutine InitAllocate(this, bounds) allocate(this%h2osfc_state_after (begg:endg)) ; this%h2osfc_state_after (:) = nan allocate(this%h2ocan_state_after (begg:endg)) ; this%h2ocan_state_after (:) = nan ! END Yorck - +#endif allocate(this%total_plant_stored_h2o_col(begc:endc)) ; this%total_plant_stored_h2o_col(:) = nan allocate(this%snw_rds_col (begc:endc,-nlevsno+1:0)) ; this%snw_rds_col (:,:) = nan @@ -488,12 +489,13 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='total water storage', & ptr_lnd=this%tws_grc) +#ifdef USE_PDAF ! Yorck: add also TWS_hactive in history outputs this%tws_hactive(begg:endg) = spval call hist_addfld1d (fname='TWS_hactive', units='mm', & avgflag='A', long_name='total water storage of hydrological active cells', & ptr_lnd=this%tws_hactive) - +#endif ! (rgk 02-02-2017) There is intentionally no entry here for stored plant water ! I think that since the value is zero in all cases except ! for FATES plant hydraulics, it will be confusing for users diff --git a/src/clm5/cpl/lnd_comp_mct.F90 b/src/clm5/cpl/lnd_comp_mct.F90 index 3ccfc10483..56179beab6 100644 --- a/src/clm5/cpl/lnd_comp_mct.F90 +++ b/src/clm5/cpl/lnd_comp_mct.F90 @@ -725,7 +725,7 @@ subroutine lnd_handle_resume( infodata ) ! Otherwise restart was modified and we are resuming from data assimulation else resume_from_data_assim = .true. - !write(iulog,*) 'resume_from_DA ', resume_from_data_assim + write(iulog,*) 'resume_from_DA ', resume_from_data_assim end if if ( resume_from_data_assim ) call update_DA_nstep() diff --git a/src/clm5/main/atm2lndType.F90 b/src/clm5/main/atm2lndType.F90 index 4f64d41348..ea52e6f35e 100644 --- a/src/clm5/main/atm2lndType.F90 +++ b/src/clm5/main/atm2lndType.F90 @@ -147,11 +147,11 @@ module atm2lndType real(r8) , pointer :: t_mo_patch (:) => null() ! patch 30-day average temperature (Kelvin) real(r8) , pointer :: t_mo_min_patch (:) => null() ! patch annual min of t_mo (Kelvin) - +#ifdef USE_PDAF ! Yorck real(r8), pointer :: volr_grc_inc (:) => null() ! rof volr total volume increment (m3) real(r8), pointer :: volr_grc_mean (:) => null() ! rof volr total volume mean (m3) - +#endif contains procedure, public :: Init @@ -558,9 +558,11 @@ subroutine InitAllocate(this, bounds) allocate(this%pfl_h2osoi_liq_grc (begg:endg,1:nlevgrnd)); this%pfl_h2osoi_liq_grc (:,:) = ival #endif +#ifdef USE_PDAF ! Yorck allocate(this%volr_grc_inc (begg:endg)) ; this%volr_grc_inc (:) = ival allocate(this%volr_grc_mean (begg:endg)) ; this%volr_grc_mean (:) = ival +#endif ! anomaly forcing allocate(this%bc_precip_grc (begg:endg)) ; this%bc_precip_grc (:) = ival diff --git a/src/clm5/main/clm_driver.F90 b/src/clm5/main/clm_driver.F90 index e2c47fb181..b7f36c920d 100644 --- a/src/clm5/main/clm_driver.F90 +++ b/src/clm5/main/clm_driver.F90 @@ -1067,7 +1067,10 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro waterstate_inst, waterflux_inst, irrigation_inst, energyflux_inst, & solarabs_inst, drydepvel_inst, & vocemis_inst, fireemis_inst, dust_inst, ch4_inst, glc_behavior, & - lnd2atm_inst, soilhydrology_inst, soilstate_inst, & + lnd2atm_inst, & +#ifdef USE_PDAF + soilhydrology_inst, soilstate_inst, & +#endif net_carbon_exchange_grc = net_carbon_exchange_grc(bounds_proc%begg:bounds_proc%endg)) deallocate(net_carbon_exchange_grc) call t_stopf('lnd2atm') diff --git a/src/clm5/main/clm_varcon.F90 b/src/clm5/main/clm_varcon.F90 index dc4f776424..3ef51256fc 100644 --- a/src/clm5/main/clm_varcon.F90 +++ b/src/clm5/main/clm_varcon.F90 @@ -106,14 +106,14 @@ module clm_varcon ! Keep this negative to avoid conflicts with possible valid values integer , public, parameter :: ispval = -9999 ! special value for int data - +#ifdef USE_PDAF ! ------------------------------------------------------------------------ ! Yorck variables ! ------------------------------------------------------------------------ integer :: averaging_var = 0 ! averaging integer for averaging TWS variables according to GRACE interval integer :: set_averaging_to_zero = ispval ! integer indicating at which point in time the averageTemp has to be set to zero due to missing observation values - +#endif ! ------------------------------------------------------------------------ ! These are tunable constants from clm2_3 ! ------------------------------------------------------------------------ diff --git a/src/clm5/main/lnd2atmMod.F90 b/src/clm5/main/lnd2atmMod.F90 index cfcf94e81e..7d1d68602d 100644 --- a/src/clm5/main/lnd2atmMod.F90 +++ b/src/clm5/main/lnd2atmMod.F90 @@ -12,7 +12,10 @@ module lnd2atmMod use shr_megan_mod , only : shr_megan_mechcomps_n use shr_fire_emis_mod , only : shr_fire_emis_mechcomps_n use clm_varpar , only : numrad, ndst, nlevgrnd, nlevsoi !ndst = number of dust bins. - use clm_varcon , only : rair, grav, cpair, hfus, tfrz, spval, set_averaging_to_zero, averaging_var, aquifer_water_baseline + use clm_varcon , only : rair, grav, cpair, hfus, tfrz, spval +#ifdef USE_PDAF + use clm_varcon , only : set_averaging_to_zero, averaging_var, aquifer_water_baseline +#endif use clm_varctl , only : iulog, use_lch4 use seq_drydep_mod , only : n_drydep, drydep_method, DD_XLND use decompMod , only : bounds_type @@ -38,11 +41,12 @@ module lnd2atmMod use LandunitType , only : lun use GridcellType , only : grc use landunit_varcon , only : istice_mec, istsoil, istcrop +#ifdef USE_PDAF use clm_time_manager , only : get_nstep use SoilHydrologyType , only : soilhydrology_type use SoilStateType , only : soilstate_type use PatchType , only : patch - +#endif ! ! !PUBLIC TYPES: implicit none @@ -131,7 +135,10 @@ subroutine lnd2atm(bounds, & waterstate_inst, waterflux_inst, irrigation_inst, energyflux_inst, & solarabs_inst, drydepvel_inst, & vocemis_inst, fireemis_inst, dust_inst, ch4_inst, glc_behavior, & - lnd2atm_inst, soilhydrology_inst, soilstate_inst, & + lnd2atm_inst, & +#ifdef USE_PDAF + soilhydrology_inst, soilstate_inst, & +#endif net_carbon_exchange_grc) ! ! !DESCRIPTION: @@ -158,23 +165,31 @@ subroutine lnd2atm(bounds, & type(ch4_type) , intent(in) :: ch4_inst type(glc_behavior_type) , intent(in) :: glc_behavior type(lnd2atm_type) , intent(inout) :: lnd2atm_inst +#ifdef USE_PDAF ! Yorck type(soilhydrology_type) , intent(inout) :: soilhydrology_inst type(soilstate_type) , intent(inout) :: soilstate_inst ! end Yorck +#endif real(r8) , intent(in) :: net_carbon_exchange_grc( bounds%begg: ) ! net carbon exchange between land and atmosphere, positive for source (gC/m2/s) ! ! !LOCAL VARIABLES: - integer :: c, g, j, p, l, index, counter ! indices + integer :: c, g, j ! indices +#ifdef USE_PDAF + integer :: p, l, index, counter ! indices +#endif real(r8) :: qflx_ice_runoff_col(bounds%begc:bounds%endc) ! total column-level ice runoff real(r8) :: eflx_sh_ice_to_liq_grc(bounds%begg:bounds%endg) ! sensible heat flux generated from the ice to liquid conversion, averaged to gridcell real(r8), parameter :: amC = 12.0_r8 ! Atomic mass number for Carbon real(r8), parameter :: amO = 16.0_r8 ! Atomic mass number for Oxygen real(r8), parameter :: amCO2 = amC + 2.0_r8*amO ! Atomic mass number for CO2 +#ifdef USE_PDAF real(r8), parameter :: m_per_mm = 1.e-3_r8 ! 0.001 meters per mm real(r8), parameter :: sec_per_hr = 3600 ! 3600 s in 1 hour +#endif ! The following converts g of C to kg of CO2 real(r8), parameter :: convertgC2kgCO2 = 1.0e-3_r8 * (amCO2/amC) +#ifdef USE_PDAF integer :: nstep ! time step number REAL, allocatable :: h2osoi_liq_grc(:) @@ -184,6 +199,7 @@ subroutine lnd2atm(bounds, & REAL, allocatable :: h2ocan_grc(:) logical, allocatable :: found(:) logical, allocatable :: found_patch(:) +#endif !------------------------------------------------------------------------ SHR_ASSERT_ALL((ubound(net_carbon_exchange_grc) == (/bounds%endg/)), errMsg(sourcefile, __LINE__)) @@ -462,6 +478,7 @@ subroutine lnd2atm(bounds, & enddo #endif +#ifdef USE_PDAF ! Yorck: update also TWS of hydrological active cells if (allocated(h2osoi_liq_grc)) deallocate(h2osoi_liq_grc) @@ -621,6 +638,7 @@ subroutine lnd2atm(bounds, & waterstate_inst%h2ocan_patch_mean(p) = (waterstate_inst%h2ocan_patch(p)+(averaging_var-1)*waterstate_inst%h2ocan_patch_mean(p))/averaging_var waterstate_inst%snocan_patch_mean(p) = (waterstate_inst%snocan_patch(p)+(averaging_var-1)*waterstate_inst%snocan_patch_mean(p))/averaging_var end do +#endif end subroutine lnd2atm From 916b72fa2873c9c2b753852392349ab785403db2 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 18 Nov 2025 13:54:48 +0100 Subject: [PATCH 17/21] transfer information on noise perturbation from internal repository - README description - correlatedNoise.py script for generating NetCDF files - example stream files --- docs/users_guide/introduction/perturbation.md | 222 +++++++++ tools/correlatedNoise.py | 433 ++++++++++++++++++ 2 files changed, 655 insertions(+) create mode 100644 tools/correlatedNoise.py diff --git a/docs/users_guide/introduction/perturbation.md b/docs/users_guide/introduction/perturbation.md index 2a747327c8..2041900c8a 100644 --- a/docs/users_guide/introduction/perturbation.md +++ b/docs/users_guide/introduction/perturbation.md @@ -6,6 +6,27 @@ The implementation of the perturbation routines was first introduced by Yorck Ewerdwalbesloh in https://gitlab.jsc.fz-juelich.de/detect/cluster-c/c01/perturbationroutineclm5. +## Description + +Temperature, precipitation as well as long and short wave radiation +are perturbed. The variables are correlated to each other to ensure +physical consistency (e.g. positive perturbation of incoming shortwave +radiation is related to a negative perturbation of longwave radiation +and a positive perturbation of temperature). For this, the same cross +correlations were used as in Reichle et al. (2007) and Han et +al. (2014). The perturbations are further correlated in space and in +time to ensure that it is still kind of "smooth" and does not look +like random noise. Spatial correlation between two grid cells was +generated by assuming an isotropic correlation structure depending on +the distance between the two grid points. The correlation distance was +set differently for precipitation than for the other variables as +precipitation is much more local (however, the values can be easily +adjusted). The variables are also correlated temporal using a +first-order auto-regressive model following Evensen (2009). There the +temporal persistence was set so that the de-correlation time is 1 day +(again easily adjustable). The perturbation values are then read in +CLM5/eCLM and the variables are perturbed inside the model. + ## Key Components ## ### 1. Soil Parameter Perturbation (`SoilStateInitTimeConstMod.F90`) ### @@ -92,4 +113,205 @@ CESM's multi-instance framework. - `src/csm_share/streams/shr_dmodel_mod.F90` - `src/datm/datm_comp_mod.F90` +## Preparing atmospheric forcing noise + +Information on how to prepare atmospherice forcing noise. + +### Python script generating atmospheric forcing noise in NetCDF files + +See `tools/correlatedNoise.py`. + +Information on `correlatedNoise.py`: +- generates one noise file per month in which all ensemble members are considered, the number has to be specified +- accounts for correlations in space and time +- paths have to be adjusted +- dt has to be adjusted (3 for 3 hourly forcing data, 1 for 1 hourly forcing data,...) +- ds has to be adjusted (12.5 for 12.5 km grid, 3 for 3 km grid) +- nn could also be adjusted, major limiation: compute power (number of points which are used to generate noise, from those, the values are upscaled to the grid) +- rho could be adjusted depending on time resolution of forcings +- covariance fields and correlations can be adjusted + + +### Example stream files for atmospheric forcing noise ### + +Three example streamfiles for noise are found in the following subsections. It is +important to specify the variables at the bottom + +- timeInformation: which time resolution do the forcing files have for which the noise files are generated --> e.g. 3 for 3 hours +- caseId: caseId for the ensemble member that the noise is for --> e.g. 0 for the first ensemble member +- numEns: the number of ensemble members for which the perturbation file was generated +- variables for perturbation: tbot_noise, lwdn_noise, precn_noise, swdn_noise (temperaure, longwave radiation, precipitation, shortwave radiation) + + +#### user_datm.streams.precip_noise.stream_0000.txt #### + +```xml + + GENERIC + + + + time time + xc lon + yc lat + area area + mask mask + + + /path/to/domain/file + + + domain.lnd.EUR-11_EUR-11.230216_mask.nc + + + + + PRECTmms precn_noise + + + /path/to/noise/files + + +2003-01.nc +2003-02.nc +2003-03.nc +2003-04.nc +2003-05.nc +2003-06.nc +2003-07.nc +2003-08.nc +2003-09.nc +2003-10.nc +2003-11.nc +2003-12.nc + + + 0 + + + 3 + + + 0 + + + 64 + + +``` + +#### user_datm.streams.solar_noise.stream_0000.txt #### + +```xml + + GENERIC + + + + time time + xc lon + yc lat + area area + mask mask + + + /path/to/domain/file + + + domain.lnd.EUR-11_EUR-11.230216_mask.nc + + + + + FSDS swdn_noise + + + /path/to/noise/files + + +2003-01.nc +2003-02.nc +2003-03.nc +2003-04.nc +2003-05.nc +2003-06.nc +2003-07.nc +2003-08.nc +2003-09.nc +2003-10.nc +2003-11.nc +2003-12.nc + + + 0 + + + 3 + + + 0 + + + 64 + + +``` + +#### user_datm.streams.other_noise.stream_0000.txt #### + +```xml + + GENERIC + + + + time time + xc lon + yc lat + area area + mask mask + + + /path/to/domain/file + + + domain.lnd.EUR-11_EUR-11.230216_mask.nc + + + + + TBOT tbot_noise + FLDS lwdn_noise + + + /path/to/noise/files + + +2003-01.nc +2003-02.nc +2003-03.nc +2003-04.nc +2003-05.nc +2003-06.nc +2003-07.nc +2003-08.nc +2003-09.nc +2003-10.nc +2003-11.nc +2003-12.nc + + + 0 + + + 3 + + + 0 + + + 64 + + +``` diff --git a/tools/correlatedNoise.py b/tools/correlatedNoise.py new file mode 100644 index 0000000000..8d580c8047 --- /dev/null +++ b/tools/correlatedNoise.py @@ -0,0 +1,433 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Tue Feb 28 14:17:32 2023 + +@author: Yorck Ewerdwalbesloh +""" +# %% Description +# This script generates spatially and temporally correlated noise and correlation of the variables +# Four variables are disturbed: +# - Surface precipitation, log-normal multiplicative +# - temperature at lowest atmospheric level, additive +# - solar radiation, log-normal multilicative +# - long wave radiation, additive + + +# %% Modules +import numpy as np +import scipy.io +import netCDF4 as nc +from scipy.interpolate import griddata +import os +import datetime +from scipy.spatial import distance +import calendar +import json +import pdb + + +def copy_attr_dim(src, dst, usr=None): + # copy attributes + for name in src.ncattrs(): + dst.setncattr("original_attribute_" + name, src.getncattr(name)) + # copy dimensions + for name, dimension in src.dimensions.items(): + dst.createDimension(name, len(dimension)) + # Additional attribute + if usr is None: + usr = os.getlogin() + dst.setncattr("perturbed_by", usr) + dst.setncattr("perturbed_on_date", datetime.datetime.today().strftime("%d.%m.%y")) + + +def compute_distances(e, f): + # Create a 2D array of (x, y) coordinates + points = np.column_stack((e, f)) + dist = distance.cdist(points, points, "euclidean") + + return dist + + +def rnd_state_serialize(): + tmp_state = np.random.get_state() + save_state = () + for i in tmp_state: + if type(i) is np.ndarray: + save_state = save_state + (i.tolist(),) + else: + save_state = save_state + (i,) + json.dump(save_state, open("rnd_state.json", "w")) + + +def rnd_state_deserialize(): + tmp_state = json.load(open("rnd_state.json", "r")) + load_state = () + for i in tmp_state: + if type(i) is list: + load_state = load_state + (np.array(i),) + else: + load_state = load_state + (i,) + np.random.set_state(load_state) + + +# Log normal to normal and vice versa +# for standard deviation and mean (formula from wikipedia) +def ln_to_n(sd_ln, mean_ln): + term = 1.0 + sd_ln * sd_ln / mean_ln / mean_ln + return (np.sqrt(np.log(term)), np.log(mean_ln / np.sqrt(term))) + + +def n_to_ln(sd_n, mean_n): + return ( + (np.exp(sd_n * sd_n) - 1.0) * np.exp(2.0 * mean_n + sd_n * sd_n), + np.exp(mean_n + sd_n * sd_n / 2.0), + ) + + +# %% Parameters +outputpath = "/path/to/noise/" +if outputpath == "/path/to/noise/": + raise ValueError("ERROR: Please set 'outputpath' to your own output directory path!") + +dt = 3 # time interval --> 3 hours for ERA5 data --> if different, this has to be adjusted + +n_ens = 64 # number of ensemble members + +# Correlation Length +L_P = 80 # km for precipitation +L_T = 250 # km for temperature, solar radiation and long wave radiation + +# Persistance parameter +rho = ( + 7 / 8 +) # decorrelation time 24h. time step of forcings 3h --> decorrelation after one days which means 3h are 1/8 --> 7/8 decorrelation time step + +# Covariance matrix for fields +C_var = np.array( + [[1, 0, -0.8, 0.5], [0, 1, 0.4, 0.4], [-0.8, 0.4, 1, -0.5], [0.5, 0.4, -0.5, 1]] +) + +# standard deviations and mean, for precipitation and short wave radiations, log normal distributions are there. + +std = np.array([ln_to_n(0.3, 1.0)[0], 2, ln_to_n(0.3, 1.0)[0], 30]) +mean = np.array([ln_to_n(0.3, 1.0)[1], 0, ln_to_n(0.3, 1.0)[1], 0]) + + +# min and max for disturbance +trunc_PREC = np.array([0.3, 1.7]) +trunc_TBOT = np.array([-5, 5]) +trunc_FSDS = np.array([0.3, 1.7]) +trunc_FLDS = np.array([-70, 70]) + +# Distance between grid points +ds = 12.5 + +# number of points for generating noise +nn = 100 + +# Specify the range of years (e.g., 2003 to 2020) +start_year = 2002 +end_year = 2022 + +# Create the year and month arrays +years = [year for year in range(start_year, end_year + 1) for month in range(1, 13)] +months = [month for year in range(start_year, end_year + 1) for month in range(1, 13)] + + +rnd_state_file = "/path/to/random/state/rnd_state.json" +if rnd_state_file == "/path/to/random/state/rnd_state.json": + raise ValueError("ERROR: Please set 'rnd_state_file' to your own random state file path!") + +force_seed = True +# Either seed random number generator or continue with existing state +if not os.path.isfile(rnd_state_file) or force_seed: + np.random.seed(42) +else: + rnd_state_deserialize() + +print("Starting") + + +# %% Cholesky decomposition of covariance matrix for variable correlation +R_var = np.linalg.cholesky(C_var).T + +# %% Grid of atmospheric data +fname = "/path/to/forcings/era5/2011-01.nc" +if fname == "/path/to/forcings/era5/2011-01.nc": + raise ValueError("ERROR: Please set 'fname' to your own forcing data path!") + +with nc.Dataset(fname) as src: + dim_time = src.dimensions["time"].size + dim_lat = src.dimensions["lat"].size + dim_lon = src.dimensions["lon"].size + + +# %% grid of input data (not real longitude and latitude but only grid on which the noise is simulated) +X = np.arange(0, dim_lon * ds, ds) +Y = np.arange(0, dim_lat * ds, ds) + +c1, d1 = np.meshgrid(X, Y) +c = np.transpose(c1).ravel() +d = np.transpose(d1).ravel() + + +# %% grid for generating noise +X = np.linspace(np.min(c), np.max(c), nn) # grid points +Y = np.linspace(np.min(d), np.max(d), nn) + +e1, f1 = np.meshgrid(X, Y) # raster +e = np.transpose(e1).ravel() +f = np.transpose(f1).ravel() + +# %% generate Matrix for spatial correlation, distances between all grid points have to be computed + +dist = compute_distances(e, f) + +# precipitation --> shorter scales (weather) + +C_P = np.exp(-(dist / L_P)) +R_P = np.linalg.cholesky(C_P).T + +# temperature +C_T = np.exp(-(dist / L_T)) +R_T = np.linalg.cholesky(C_T).T + +# %% Generate random noise for each variable +e_PREC = np.random.randn(len(e)) +e_TBOT = np.random.randn(len(e)) +e_FSDS = np.random.randn(len(e)) +e_FLDS = np.random.randn(len(e)) + +# %% Correlations between variables +temp = np.array((e_PREC, e_TBOT, e_FSDS, e_FLDS)).T.dot(R_var) +e_PREC = temp[:, 0] +e_TBOT = temp[:, 1] +e_FSDS = temp[:, 2] +e_FLDS = temp[:, 3] +del temp + +# %% Correlate random noise +e_corrSpat_PREC = R_P.T.dot(e_PREC) +e_corrSpat_TBOT = R_T.T.dot(e_TBOT) +e_corrSpat_FSDS = R_T.T.dot(e_FSDS) +e_corrSpat_FLDS = R_T.T.dot(e_FLDS) + + +# %% Extend to large grid +e_corrSpat_PREC = griddata((e, f), e_corrSpat_PREC, (c, d), "linear") +e_add_PREC = 0.2 * np.random.randn(dim_lon * dim_lat, 1) +e_corrSpat_PREC = e_corrSpat_PREC.ravel() + e_add_PREC.ravel() + +e_corrSpat_TBOT = griddata((e, f), e_corrSpat_TBOT, (c, d), "linear") +e_add_TBOT = 0.2 * np.random.randn(dim_lon * dim_lat, 1) +e_corrSpat_TBOT = e_corrSpat_TBOT.ravel() + e_add_TBOT.ravel() + +e_corrSpat_FSDS = griddata((e, f), e_corrSpat_FSDS, (c, d), "linear") +e_add_FSDS = 0.2 * np.random.randn(dim_lon * dim_lat, 1) +e_corrSpat_FSDS = e_corrSpat_FSDS.ravel() + e_add_FSDS.ravel() + +e_corrSpat_FLDS = griddata((e, f), e_corrSpat_FLDS, (c, d), "linear") +e_add_FLDS = 0.2 * np.random.randn(dim_lon * dim_lat, 1) +e_corrSpat_FLDS = e_corrSpat_FLDS.ravel() + e_add_FLDS.ravel() + +if not force_seed: + rnd_state_serialize() + +e_prev_PREC = e_corrSpat_PREC +e_prev_TBOT = e_corrSpat_TBOT +e_prev_FSDS = e_corrSpat_FSDS +e_prev_FLDS = e_corrSpat_FLDS + +del e_corrSpat_PREC, e_corrSpat_TBOT, e_corrSpat_FSDS, e_corrSpat_FLDS +del e_add_PREC, e_add_TBOT, e_add_FSDS, e_add_FLDS + +PREC_old = 0 +TBOT_old = 0 +FSDS_old = 0 +FLDS_old = 0 +print("Looping over years") + +i = 0 +for year, month in zip(years, months * len(years)): + print(str(year) + "-" + str(month).zfill(2)) + filename = outputpath + str(year) + "-" + str(month).zfill(2) + ".nc" + + os.makedirs(os.path.dirname(filename), exist_ok=True) + fname = ( + "/path/to/forcings/era5/" + + str(year) + + "-" + + str(month).zfill(2) + + ".nc" + ) + if fname.startswith("/path/to/forcings/era5/"): + raise ValueError("ERROR: Please update the forcing data path in the loop to your own directory!") + + # number of extra entries due to ensemble run + timeExtra = (n_ens - 1) * (24 / dt) + + # number of time steps for files + timeNum = calendar.monthrange(year, month)[1] * (24 / dt) + timeExtra + # pdb.set_trace() + with nc.Dataset(filename, "w") as ncid: + dimid_lon = ncid.createDimension("lon", dim_lon) + dimid_lat = ncid.createDimension("lat", dim_lat) + dimid_time = ncid.createDimension("time", timeNum) + with nc.Dataset(fname) as src: + for name, var in src.variables.items(): + if name == "latitude" or name == "longitude": + nvar = ncid.createVariable(name, var.datatype, ("lat", "lon")) + ncid[name].setncatts(src[name].__dict__) + ncid[name][:] = src[name][:] + if name == "time": + nvar = ncid.createVariable(name, var.datatype, ("time")) + ncid[name].setncatts(src[name].__dict__) + ncid[name][: len(src[name][:])] = src[name][:] + start = 90000000 + length = len(ncid[name][len(src[name][:]) :]) + ncid[name][len(src[name][:]) :] = [start + i for i in range(length)] + + PRECTmms_ID = ncid.createVariable("PRECTmms", "f8", ("time", "lat", "lon")) + TBOT_ID = ncid.createVariable("TBOT", "f8", ("time", "lat", "lon")) + FSDS_ID = ncid.createVariable("FSDS", "f8", ("time", "lat", "lon")) + FLDS_ID = ncid.createVariable("FLDS", "f8", ("time", "lat", "lon")) + + # %% Loops over time + PREC = np.zeros((dim_lon, dim_lat, int(timeNum))) + TBOT = PREC.copy() + FSDS = PREC.copy() + FLDS = PREC.copy() + t_start = 0 + + if i > 0: + PREC[:, :, : int(timeExtra)] = PREC_old[:, :, -int(timeExtra) :] + TBOT[:, :, : int(timeExtra)] = TBOT_old[:, :, -int(timeExtra) :] + FSDS[:, :, : int(timeExtra)] = FSDS_old[:, :, -int(timeExtra) :] + FLDS[:, :, : int(timeExtra)] = FLDS_old[:, :, -int(timeExtra) :] + t_start += timeExtra + + for t in range(int(t_start), int(timeNum)): + + print("Time Step ", t + 1, "/", timeNum) + # %% Random numbers + e_PREC = np.random.randn(len(e)) + e_TBOT = np.random.randn(len(e)) + e_FSDS = np.random.randn(len(e)) + e_FLDS = np.random.randn(len(e)) + + # %% Correlation between variables on small grid + lok = np.array((e_PREC, e_TBOT, e_FSDS, e_FLDS)).T + temp = lok.dot(R_var) + del lok + coeff_corr = np.corrcoef(temp.T) + e_PREC = temp[:, 0] + e_TBOT = temp[:, 1] + e_FSDS = temp[:, 2] + e_FLDS = temp[:, 3] + + # %% Correlate random noise on small grid + e_corrSpat_PREC = R_P.T.dot(e_PREC) + e_corrSpat_TBOT = R_T.T.dot(e_TBOT) + e_corrSpat_FSDS = R_T.T.dot(e_FSDS) + e_corrSpat_FLDS = R_T.T.dot(e_FLDS) + + del e_PREC, e_TBOT, e_FSDS, e_FLDS + + # %% Extend to large grid + e_corrSpat_PREC = griddata((e, f), e_corrSpat_PREC, (c, d), "linear") + e_add_PREC = 0.2 * np.random.randn(dim_lon * dim_lat, 1) + e_corrSpat_PREC = e_corrSpat_PREC.ravel() + e_add_PREC.ravel() + + e_corrSpat_TBOT = griddata((e, f), e_corrSpat_TBOT, (c, d), "linear") + e_add_TBOT = 0.2 * np.random.randn(dim_lon * dim_lat, 1) + e_corrSpat_TBOT = e_corrSpat_TBOT.ravel() + e_add_TBOT.ravel() + + e_corrSpat_FSDS = griddata((e, f), e_corrSpat_FSDS, (c, d), "linear") + e_add_FSDS = 0.2 * np.random.randn(dim_lon * dim_lat, 1) + e_corrSpat_FSDS = e_corrSpat_FSDS.ravel() + e_add_FSDS.ravel() + + e_corrSpat_FLDS = griddata((e, f), e_corrSpat_FLDS, (c, d), "linear") + e_add_FLDS = 0.2 * np.random.randn(dim_lon * dim_lat, 1) + e_corrSpat_FLDS = e_corrSpat_FLDS.ravel() + e_add_FLDS.ravel() + + del e_add_PREC, e_add_TBOT, e_add_FSDS, e_add_FLDS + + if not force_seed: + rnd_state_serialize() + + # %% Correlate temporally + e_corrSpatTemp_PREC = ( + rho * e_prev_PREC + np.sqrt(1 - rho**2) * e_corrSpat_PREC + ) + e_corrSpatTemp_TBOT = ( + rho * e_prev_TBOT + np.sqrt(1 - rho**2) * e_corrSpat_TBOT + ) + e_corrSpatTemp_FSDS = ( + rho * e_prev_FSDS + np.sqrt(1 - rho**2) * e_corrSpat_FSDS + ) + e_corrSpatTemp_FLDS = ( + rho * e_prev_FLDS + np.sqrt(1 - rho**2) * e_corrSpat_FLDS + ) + del e_corrSpat_PREC, e_corrSpat_TBOT, e_corrSpat_FSDS, e_corrSpat_FLDS + del e_prev_PREC, e_prev_TBOT, e_prev_FSDS, e_prev_FLDS + + # %% new previous fields + e_prev_PREC = e_corrSpatTemp_PREC + e_prev_TBOT = e_corrSpatTemp_TBOT + e_prev_FSDS = e_corrSpatTemp_FSDS + e_prev_FLDS = e_corrSpatTemp_FLDS + + # %% Truncate noise + e_corrSpatTemp_PREC = np.exp(mean[0] + std[0] * e_corrSpatTemp_PREC) + e_corrSpatTemp_PREC[e_corrSpatTemp_PREC < trunc_PREC[0]] = trunc_PREC[0] + e_corrSpatTemp_PREC[e_corrSpatTemp_PREC > trunc_PREC[1]] = trunc_PREC[1] + + e_corrSpatTemp_TBOT = std[1] * e_corrSpatTemp_TBOT + e_corrSpatTemp_TBOT[e_corrSpatTemp_TBOT < trunc_TBOT[0]] = trunc_TBOT[0] + e_corrSpatTemp_TBOT[e_corrSpatTemp_TBOT > trunc_TBOT[1]] = trunc_TBOT[1] + + e_corrSpatTemp_FSDS = np.exp(mean[2] + std[2] * e_corrSpatTemp_FSDS) + e_corrSpatTemp_FSDS[e_corrSpatTemp_FSDS < trunc_FSDS[0]] = trunc_FSDS[0] + e_corrSpatTemp_FSDS[e_corrSpatTemp_FSDS > trunc_FSDS[1]] = trunc_FSDS[1] + + e_corrSpatTemp_FLDS = std[3] * e_corrSpatTemp_FLDS + e_corrSpatTemp_FLDS[e_corrSpatTemp_FLDS < trunc_FLDS[0]] = trunc_FLDS[0] + e_corrSpatTemp_FLDS[e_corrSpatTemp_FLDS > trunc_FLDS[1]] = trunc_FLDS[1] + + # %% Apply to data + PREC[:, :, t] = np.reshape(e_corrSpatTemp_PREC, (dim_lon, dim_lat)) + TBOT[:, :, t] = np.reshape(e_corrSpatTemp_TBOT, (dim_lon, dim_lat)) + FSDS[:, :, t] = np.reshape(e_corrSpatTemp_FSDS, (dim_lon, dim_lat)) + FLDS[:, :, t] = np.reshape(e_corrSpatTemp_FLDS, (dim_lon, dim_lat)) + + del ( + e_corrSpatTemp_PREC, + e_corrSpatTemp_TBOT, + e_corrSpatTemp_FSDS, + e_corrSpatTemp_FLDS, + ) + + # %% Save disturbed variables + # Write PREC variable + + PRECTmms_ID[:] = np.transpose(PREC, (2, 1, 0)) + + # Write TBOT variable + TBOT_ID[:] = np.transpose(TBOT, (2, 1, 0)) + + # Write FSDS variable + FSDS_ID[:] = np.transpose(FSDS, (2, 1, 0)) + + # Write FLDS variable + FLDS_ID[:] = np.transpose(FLDS, (2, 1, 0)) + + # pdb.set_trace() + PREC_old = PREC + TBOT_old = TBOT + FSDS_old = FSDS + FLDS_old = FLDS + + del PREC, TBOT, FSDS, FLDS + + i = i + 1 From bbbdf65b84486fad759faa6a889929259fc25d93 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 18 Nov 2025 14:04:36 +0100 Subject: [PATCH 18/21] docs: smoothened text with AI tool --- docs/users_guide/introduction/perturbation.md | 122 ++++++++++++------ 1 file changed, 81 insertions(+), 41 deletions(-) diff --git a/docs/users_guide/introduction/perturbation.md b/docs/users_guide/introduction/perturbation.md index 2041900c8a..3dcd2f217d 100644 --- a/docs/users_guide/introduction/perturbation.md +++ b/docs/users_guide/introduction/perturbation.md @@ -8,24 +8,40 @@ https://gitlab.jsc.fz-juelich.de/detect/cluster-c/c01/perturbationroutineclm5. ## Description -Temperature, precipitation as well as long and short wave radiation -are perturbed. The variables are correlated to each other to ensure -physical consistency (e.g. positive perturbation of incoming shortwave -radiation is related to a negative perturbation of longwave radiation -and a positive perturbation of temperature). For this, the same cross -correlations were used as in Reichle et al. (2007) and Han et -al. (2014). The perturbations are further correlated in space and in -time to ensure that it is still kind of "smooth" and does not look -like random noise. Spatial correlation between two grid cells was -generated by assuming an isotropic correlation structure depending on -the distance between the two grid points. The correlation distance was -set differently for precipitation than for the other variables as -precipitation is much more local (however, the values can be easily -adjusted). The variables are also correlated temporal using a -first-order auto-regressive model following Evensen (2009). There the -temporal persistence was set so that the de-correlation time is 1 day -(again easily adjustable). The perturbation values are then read in -CLM5/eCLM and the variables are perturbed inside the model. +eCLM implements perturbations for atmospheric forcing variables +(temperature, precipitation, longwave radiation, and shortwave +radiation) to support ensemble-based data assimilation in +eCLM-PDAF. The perturbation scheme maintains physical consistency +through cross-variable correlations and spatiotemporal coherence. + +### Physical Consistency + +Perturbed variables are cross-correlated to preserve physical +relationships between forcing fields. For example, a positive +perturbation of incoming shortwave radiation corresponds to a negative +perturbation of longwave radiation and a positive perturbation of +temperature. The cross-correlation structure follows Reichle et +al. (2007) and Han et al. (2014). + +### Spatial Correlation + +Perturbations are spatially correlated using an isotropic correlation +function based on grid cell separation distance. Precipitation uses a +shorter correlation length scale than other variables to reflect its +more localized spatial structure. Correlation parameters are +configurable. + +### Temporal Correlation + +Temporal correlation is generated using a first-order autoregressive +(AR-1) model following Evensen (2009). The default temporal +persistence corresponds to a decorrelation timescale of one day, which +is adjustable based on application requirements. + +### Implementation + +Perturbation fields are read from preprocessed noise files and applied +to atmospheric forcing variables within the eCLM model during runtime. ## Key Components ## @@ -113,37 +129,57 @@ CESM's multi-instance framework. - `src/csm_share/streams/shr_dmodel_mod.F90` - `src/datm/datm_comp_mod.F90` -## Preparing atmospheric forcing noise +## Preparing Atmospheric Forcing Noise -Information on how to prepare atmospherice forcing noise. +This section describes the workflow for generating and configuring +spatiotemporally correlated noise fields for atmospheric forcing +perturbations. -### Python script generating atmospheric forcing noise in NetCDF files +### Noise Generation Tool -See `tools/correlatedNoise.py`. +The `tools/correlatedNoise.py` script generates monthly NetCDF files +containing spatiotemporally correlated noise fields for all ensemble +members. The script implements the correlation structures described in +the [Description](#description) section. -Information on `correlatedNoise.py`: -- generates one noise file per month in which all ensemble members are considered, the number has to be specified -- accounts for correlations in space and time -- paths have to be adjusted -- dt has to be adjusted (3 for 3 hourly forcing data, 1 for 1 hourly forcing data,...) -- ds has to be adjusted (12.5 for 12.5 km grid, 3 for 3 km grid) -- nn could also be adjusted, major limiation: compute power (number of points which are used to generate noise, from those, the values are upscaled to the grid) -- rho could be adjusted depending on time resolution of forcings -- covariance fields and correlations can be adjusted +**Key Configuration Parameters:** +- **Ensemble size**: Number of ensemble members to generate noise for (must be specified) +- **Temporal resolution** (`dt`): Time step in hours (e.g., `3` for 3-hourly data, `1` for hourly data) +- **Spatial resolution** (`ds`): Grid spacing in km (e.g., `12.5` for 12.5 km grid, `3` for 3 km grid) +- **Spatial sampling** (`nn`): Number of points used for noise generation before upscaling (limited by computational resources) +- **Temporal persistence** (`rho`): AR-1 correlation coefficient (adjust based on forcing time resolution) +- **Covariance matrix**: Cross-variable covariance structure and correlation coefficients (configurable) +- **File paths**: Input/output directories (must be adjusted for local system) -### Example stream files for atmospheric forcing noise ### +**Output**: One NetCDF file per month containing noise fields for all ensemble members and all perturbed variables. -Three example streamfiles for noise are found in the following subsections. It is -important to specify the variables at the bottom +### Stream File Configuration -- timeInformation: which time resolution do the forcing files have for which the noise files are generated --> e.g. 3 for 3 hours -- caseId: caseId for the ensemble member that the noise is for --> e.g. 0 for the first ensemble member -- numEns: the number of ensemble members for which the perturbation file was generated -- variables for perturbation: tbot_noise, lwdn_noise, precn_noise, swdn_noise (temperaure, longwave radiation, precipitation, shortwave radiation) +Stream files configure how eCLM reads and applies noise fields to +atmospheric forcing data. Three separate stream files are required +(precipitation, solar radiation, and other variables). Each stream +file must specify the following metadata fields: +**Required Metadata:** -#### user_datm.streams.precip_noise.stream_0000.txt #### +- `timeInformation`: Temporal resolution of forcing data in hours (e.g., `3` for 3-hourly data) +- `caseId`: Ensemble member index (e.g., `0` for first member, `1` for second member) +- `numEns`: Total ensemble size used when generating the noise files +- **Perturbation variables**: `tbot_noise` (temperature), `lwdn_noise` (longwave radiation), `precn_noise` (precipitation), `swdn_noise` (shortwave radiation) + + +### Stream File Examples + +The following examples demonstrate proper stream file configuration +for precipitation, solar radiation, and other atmospheric +variables. Each example must be adapted with correct file paths, +temporal parameters, and ensemble member IDs for your specific +application. + +#### Precipitation Noise Stream File + +**File**: `user_datm.streams.precip_noise.stream_0000.txt` ```xml @@ -200,7 +236,9 @@ important to specify the variables at the bottom ``` -#### user_datm.streams.solar_noise.stream_0000.txt #### +#### Solar Radiation Noise Stream File + +**File**: `user_datm.streams.solar_noise.stream_0000.txt` ```xml @@ -257,7 +295,9 @@ important to specify the variables at the bottom ``` -#### user_datm.streams.other_noise.stream_0000.txt #### +#### Temperature and Longwave Radiation Noise Stream File + +**File**: `user_datm.streams.other_noise.stream_0000.txt` ```xml From 555c1d115dad8689a0600bb0ed44b9c16116cf60 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 18 Nov 2025 14:10:03 +0100 Subject: [PATCH 19/21] docs: remove link --- docs/users_guide/introduction/perturbation.md | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/docs/users_guide/introduction/perturbation.md b/docs/users_guide/introduction/perturbation.md index 3dcd2f217d..42621d6c63 100644 --- a/docs/users_guide/introduction/perturbation.md +++ b/docs/users_guide/introduction/perturbation.md @@ -139,8 +139,7 @@ perturbations. The `tools/correlatedNoise.py` script generates monthly NetCDF files containing spatiotemporally correlated noise fields for all ensemble -members. The script implements the correlation structures described in -the [Description](#description) section. +members. **Key Configuration Parameters:** From e5cbad765a8bf952e4b474577dcc3ba7ed9dfef3 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 28 Nov 2025 18:12:13 +0100 Subject: [PATCH 20/21] remove tools,this is in atm forc gen --- tools/correlatedNoise.py | 433 --------------------------------------- 1 file changed, 433 deletions(-) delete mode 100644 tools/correlatedNoise.py diff --git a/tools/correlatedNoise.py b/tools/correlatedNoise.py deleted file mode 100644 index 8d580c8047..0000000000 --- a/tools/correlatedNoise.py +++ /dev/null @@ -1,433 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Tue Feb 28 14:17:32 2023 - -@author: Yorck Ewerdwalbesloh -""" -# %% Description -# This script generates spatially and temporally correlated noise and correlation of the variables -# Four variables are disturbed: -# - Surface precipitation, log-normal multiplicative -# - temperature at lowest atmospheric level, additive -# - solar radiation, log-normal multilicative -# - long wave radiation, additive - - -# %% Modules -import numpy as np -import scipy.io -import netCDF4 as nc -from scipy.interpolate import griddata -import os -import datetime -from scipy.spatial import distance -import calendar -import json -import pdb - - -def copy_attr_dim(src, dst, usr=None): - # copy attributes - for name in src.ncattrs(): - dst.setncattr("original_attribute_" + name, src.getncattr(name)) - # copy dimensions - for name, dimension in src.dimensions.items(): - dst.createDimension(name, len(dimension)) - # Additional attribute - if usr is None: - usr = os.getlogin() - dst.setncattr("perturbed_by", usr) - dst.setncattr("perturbed_on_date", datetime.datetime.today().strftime("%d.%m.%y")) - - -def compute_distances(e, f): - # Create a 2D array of (x, y) coordinates - points = np.column_stack((e, f)) - dist = distance.cdist(points, points, "euclidean") - - return dist - - -def rnd_state_serialize(): - tmp_state = np.random.get_state() - save_state = () - for i in tmp_state: - if type(i) is np.ndarray: - save_state = save_state + (i.tolist(),) - else: - save_state = save_state + (i,) - json.dump(save_state, open("rnd_state.json", "w")) - - -def rnd_state_deserialize(): - tmp_state = json.load(open("rnd_state.json", "r")) - load_state = () - for i in tmp_state: - if type(i) is list: - load_state = load_state + (np.array(i),) - else: - load_state = load_state + (i,) - np.random.set_state(load_state) - - -# Log normal to normal and vice versa -# for standard deviation and mean (formula from wikipedia) -def ln_to_n(sd_ln, mean_ln): - term = 1.0 + sd_ln * sd_ln / mean_ln / mean_ln - return (np.sqrt(np.log(term)), np.log(mean_ln / np.sqrt(term))) - - -def n_to_ln(sd_n, mean_n): - return ( - (np.exp(sd_n * sd_n) - 1.0) * np.exp(2.0 * mean_n + sd_n * sd_n), - np.exp(mean_n + sd_n * sd_n / 2.0), - ) - - -# %% Parameters -outputpath = "/path/to/noise/" -if outputpath == "/path/to/noise/": - raise ValueError("ERROR: Please set 'outputpath' to your own output directory path!") - -dt = 3 # time interval --> 3 hours for ERA5 data --> if different, this has to be adjusted - -n_ens = 64 # number of ensemble members - -# Correlation Length -L_P = 80 # km for precipitation -L_T = 250 # km for temperature, solar radiation and long wave radiation - -# Persistance parameter -rho = ( - 7 / 8 -) # decorrelation time 24h. time step of forcings 3h --> decorrelation after one days which means 3h are 1/8 --> 7/8 decorrelation time step - -# Covariance matrix for fields -C_var = np.array( - [[1, 0, -0.8, 0.5], [0, 1, 0.4, 0.4], [-0.8, 0.4, 1, -0.5], [0.5, 0.4, -0.5, 1]] -) - -# standard deviations and mean, for precipitation and short wave radiations, log normal distributions are there. - -std = np.array([ln_to_n(0.3, 1.0)[0], 2, ln_to_n(0.3, 1.0)[0], 30]) -mean = np.array([ln_to_n(0.3, 1.0)[1], 0, ln_to_n(0.3, 1.0)[1], 0]) - - -# min and max for disturbance -trunc_PREC = np.array([0.3, 1.7]) -trunc_TBOT = np.array([-5, 5]) -trunc_FSDS = np.array([0.3, 1.7]) -trunc_FLDS = np.array([-70, 70]) - -# Distance between grid points -ds = 12.5 - -# number of points for generating noise -nn = 100 - -# Specify the range of years (e.g., 2003 to 2020) -start_year = 2002 -end_year = 2022 - -# Create the year and month arrays -years = [year for year in range(start_year, end_year + 1) for month in range(1, 13)] -months = [month for year in range(start_year, end_year + 1) for month in range(1, 13)] - - -rnd_state_file = "/path/to/random/state/rnd_state.json" -if rnd_state_file == "/path/to/random/state/rnd_state.json": - raise ValueError("ERROR: Please set 'rnd_state_file' to your own random state file path!") - -force_seed = True -# Either seed random number generator or continue with existing state -if not os.path.isfile(rnd_state_file) or force_seed: - np.random.seed(42) -else: - rnd_state_deserialize() - -print("Starting") - - -# %% Cholesky decomposition of covariance matrix for variable correlation -R_var = np.linalg.cholesky(C_var).T - -# %% Grid of atmospheric data -fname = "/path/to/forcings/era5/2011-01.nc" -if fname == "/path/to/forcings/era5/2011-01.nc": - raise ValueError("ERROR: Please set 'fname' to your own forcing data path!") - -with nc.Dataset(fname) as src: - dim_time = src.dimensions["time"].size - dim_lat = src.dimensions["lat"].size - dim_lon = src.dimensions["lon"].size - - -# %% grid of input data (not real longitude and latitude but only grid on which the noise is simulated) -X = np.arange(0, dim_lon * ds, ds) -Y = np.arange(0, dim_lat * ds, ds) - -c1, d1 = np.meshgrid(X, Y) -c = np.transpose(c1).ravel() -d = np.transpose(d1).ravel() - - -# %% grid for generating noise -X = np.linspace(np.min(c), np.max(c), nn) # grid points -Y = np.linspace(np.min(d), np.max(d), nn) - -e1, f1 = np.meshgrid(X, Y) # raster -e = np.transpose(e1).ravel() -f = np.transpose(f1).ravel() - -# %% generate Matrix for spatial correlation, distances between all grid points have to be computed - -dist = compute_distances(e, f) - -# precipitation --> shorter scales (weather) - -C_P = np.exp(-(dist / L_P)) -R_P = np.linalg.cholesky(C_P).T - -# temperature -C_T = np.exp(-(dist / L_T)) -R_T = np.linalg.cholesky(C_T).T - -# %% Generate random noise for each variable -e_PREC = np.random.randn(len(e)) -e_TBOT = np.random.randn(len(e)) -e_FSDS = np.random.randn(len(e)) -e_FLDS = np.random.randn(len(e)) - -# %% Correlations between variables -temp = np.array((e_PREC, e_TBOT, e_FSDS, e_FLDS)).T.dot(R_var) -e_PREC = temp[:, 0] -e_TBOT = temp[:, 1] -e_FSDS = temp[:, 2] -e_FLDS = temp[:, 3] -del temp - -# %% Correlate random noise -e_corrSpat_PREC = R_P.T.dot(e_PREC) -e_corrSpat_TBOT = R_T.T.dot(e_TBOT) -e_corrSpat_FSDS = R_T.T.dot(e_FSDS) -e_corrSpat_FLDS = R_T.T.dot(e_FLDS) - - -# %% Extend to large grid -e_corrSpat_PREC = griddata((e, f), e_corrSpat_PREC, (c, d), "linear") -e_add_PREC = 0.2 * np.random.randn(dim_lon * dim_lat, 1) -e_corrSpat_PREC = e_corrSpat_PREC.ravel() + e_add_PREC.ravel() - -e_corrSpat_TBOT = griddata((e, f), e_corrSpat_TBOT, (c, d), "linear") -e_add_TBOT = 0.2 * np.random.randn(dim_lon * dim_lat, 1) -e_corrSpat_TBOT = e_corrSpat_TBOT.ravel() + e_add_TBOT.ravel() - -e_corrSpat_FSDS = griddata((e, f), e_corrSpat_FSDS, (c, d), "linear") -e_add_FSDS = 0.2 * np.random.randn(dim_lon * dim_lat, 1) -e_corrSpat_FSDS = e_corrSpat_FSDS.ravel() + e_add_FSDS.ravel() - -e_corrSpat_FLDS = griddata((e, f), e_corrSpat_FLDS, (c, d), "linear") -e_add_FLDS = 0.2 * np.random.randn(dim_lon * dim_lat, 1) -e_corrSpat_FLDS = e_corrSpat_FLDS.ravel() + e_add_FLDS.ravel() - -if not force_seed: - rnd_state_serialize() - -e_prev_PREC = e_corrSpat_PREC -e_prev_TBOT = e_corrSpat_TBOT -e_prev_FSDS = e_corrSpat_FSDS -e_prev_FLDS = e_corrSpat_FLDS - -del e_corrSpat_PREC, e_corrSpat_TBOT, e_corrSpat_FSDS, e_corrSpat_FLDS -del e_add_PREC, e_add_TBOT, e_add_FSDS, e_add_FLDS - -PREC_old = 0 -TBOT_old = 0 -FSDS_old = 0 -FLDS_old = 0 -print("Looping over years") - -i = 0 -for year, month in zip(years, months * len(years)): - print(str(year) + "-" + str(month).zfill(2)) - filename = outputpath + str(year) + "-" + str(month).zfill(2) + ".nc" - - os.makedirs(os.path.dirname(filename), exist_ok=True) - fname = ( - "/path/to/forcings/era5/" - + str(year) - + "-" - + str(month).zfill(2) - + ".nc" - ) - if fname.startswith("/path/to/forcings/era5/"): - raise ValueError("ERROR: Please update the forcing data path in the loop to your own directory!") - - # number of extra entries due to ensemble run - timeExtra = (n_ens - 1) * (24 / dt) - - # number of time steps for files - timeNum = calendar.monthrange(year, month)[1] * (24 / dt) + timeExtra - # pdb.set_trace() - with nc.Dataset(filename, "w") as ncid: - dimid_lon = ncid.createDimension("lon", dim_lon) - dimid_lat = ncid.createDimension("lat", dim_lat) - dimid_time = ncid.createDimension("time", timeNum) - with nc.Dataset(fname) as src: - for name, var in src.variables.items(): - if name == "latitude" or name == "longitude": - nvar = ncid.createVariable(name, var.datatype, ("lat", "lon")) - ncid[name].setncatts(src[name].__dict__) - ncid[name][:] = src[name][:] - if name == "time": - nvar = ncid.createVariable(name, var.datatype, ("time")) - ncid[name].setncatts(src[name].__dict__) - ncid[name][: len(src[name][:])] = src[name][:] - start = 90000000 - length = len(ncid[name][len(src[name][:]) :]) - ncid[name][len(src[name][:]) :] = [start + i for i in range(length)] - - PRECTmms_ID = ncid.createVariable("PRECTmms", "f8", ("time", "lat", "lon")) - TBOT_ID = ncid.createVariable("TBOT", "f8", ("time", "lat", "lon")) - FSDS_ID = ncid.createVariable("FSDS", "f8", ("time", "lat", "lon")) - FLDS_ID = ncid.createVariable("FLDS", "f8", ("time", "lat", "lon")) - - # %% Loops over time - PREC = np.zeros((dim_lon, dim_lat, int(timeNum))) - TBOT = PREC.copy() - FSDS = PREC.copy() - FLDS = PREC.copy() - t_start = 0 - - if i > 0: - PREC[:, :, : int(timeExtra)] = PREC_old[:, :, -int(timeExtra) :] - TBOT[:, :, : int(timeExtra)] = TBOT_old[:, :, -int(timeExtra) :] - FSDS[:, :, : int(timeExtra)] = FSDS_old[:, :, -int(timeExtra) :] - FLDS[:, :, : int(timeExtra)] = FLDS_old[:, :, -int(timeExtra) :] - t_start += timeExtra - - for t in range(int(t_start), int(timeNum)): - - print("Time Step ", t + 1, "/", timeNum) - # %% Random numbers - e_PREC = np.random.randn(len(e)) - e_TBOT = np.random.randn(len(e)) - e_FSDS = np.random.randn(len(e)) - e_FLDS = np.random.randn(len(e)) - - # %% Correlation between variables on small grid - lok = np.array((e_PREC, e_TBOT, e_FSDS, e_FLDS)).T - temp = lok.dot(R_var) - del lok - coeff_corr = np.corrcoef(temp.T) - e_PREC = temp[:, 0] - e_TBOT = temp[:, 1] - e_FSDS = temp[:, 2] - e_FLDS = temp[:, 3] - - # %% Correlate random noise on small grid - e_corrSpat_PREC = R_P.T.dot(e_PREC) - e_corrSpat_TBOT = R_T.T.dot(e_TBOT) - e_corrSpat_FSDS = R_T.T.dot(e_FSDS) - e_corrSpat_FLDS = R_T.T.dot(e_FLDS) - - del e_PREC, e_TBOT, e_FSDS, e_FLDS - - # %% Extend to large grid - e_corrSpat_PREC = griddata((e, f), e_corrSpat_PREC, (c, d), "linear") - e_add_PREC = 0.2 * np.random.randn(dim_lon * dim_lat, 1) - e_corrSpat_PREC = e_corrSpat_PREC.ravel() + e_add_PREC.ravel() - - e_corrSpat_TBOT = griddata((e, f), e_corrSpat_TBOT, (c, d), "linear") - e_add_TBOT = 0.2 * np.random.randn(dim_lon * dim_lat, 1) - e_corrSpat_TBOT = e_corrSpat_TBOT.ravel() + e_add_TBOT.ravel() - - e_corrSpat_FSDS = griddata((e, f), e_corrSpat_FSDS, (c, d), "linear") - e_add_FSDS = 0.2 * np.random.randn(dim_lon * dim_lat, 1) - e_corrSpat_FSDS = e_corrSpat_FSDS.ravel() + e_add_FSDS.ravel() - - e_corrSpat_FLDS = griddata((e, f), e_corrSpat_FLDS, (c, d), "linear") - e_add_FLDS = 0.2 * np.random.randn(dim_lon * dim_lat, 1) - e_corrSpat_FLDS = e_corrSpat_FLDS.ravel() + e_add_FLDS.ravel() - - del e_add_PREC, e_add_TBOT, e_add_FSDS, e_add_FLDS - - if not force_seed: - rnd_state_serialize() - - # %% Correlate temporally - e_corrSpatTemp_PREC = ( - rho * e_prev_PREC + np.sqrt(1 - rho**2) * e_corrSpat_PREC - ) - e_corrSpatTemp_TBOT = ( - rho * e_prev_TBOT + np.sqrt(1 - rho**2) * e_corrSpat_TBOT - ) - e_corrSpatTemp_FSDS = ( - rho * e_prev_FSDS + np.sqrt(1 - rho**2) * e_corrSpat_FSDS - ) - e_corrSpatTemp_FLDS = ( - rho * e_prev_FLDS + np.sqrt(1 - rho**2) * e_corrSpat_FLDS - ) - del e_corrSpat_PREC, e_corrSpat_TBOT, e_corrSpat_FSDS, e_corrSpat_FLDS - del e_prev_PREC, e_prev_TBOT, e_prev_FSDS, e_prev_FLDS - - # %% new previous fields - e_prev_PREC = e_corrSpatTemp_PREC - e_prev_TBOT = e_corrSpatTemp_TBOT - e_prev_FSDS = e_corrSpatTemp_FSDS - e_prev_FLDS = e_corrSpatTemp_FLDS - - # %% Truncate noise - e_corrSpatTemp_PREC = np.exp(mean[0] + std[0] * e_corrSpatTemp_PREC) - e_corrSpatTemp_PREC[e_corrSpatTemp_PREC < trunc_PREC[0]] = trunc_PREC[0] - e_corrSpatTemp_PREC[e_corrSpatTemp_PREC > trunc_PREC[1]] = trunc_PREC[1] - - e_corrSpatTemp_TBOT = std[1] * e_corrSpatTemp_TBOT - e_corrSpatTemp_TBOT[e_corrSpatTemp_TBOT < trunc_TBOT[0]] = trunc_TBOT[0] - e_corrSpatTemp_TBOT[e_corrSpatTemp_TBOT > trunc_TBOT[1]] = trunc_TBOT[1] - - e_corrSpatTemp_FSDS = np.exp(mean[2] + std[2] * e_corrSpatTemp_FSDS) - e_corrSpatTemp_FSDS[e_corrSpatTemp_FSDS < trunc_FSDS[0]] = trunc_FSDS[0] - e_corrSpatTemp_FSDS[e_corrSpatTemp_FSDS > trunc_FSDS[1]] = trunc_FSDS[1] - - e_corrSpatTemp_FLDS = std[3] * e_corrSpatTemp_FLDS - e_corrSpatTemp_FLDS[e_corrSpatTemp_FLDS < trunc_FLDS[0]] = trunc_FLDS[0] - e_corrSpatTemp_FLDS[e_corrSpatTemp_FLDS > trunc_FLDS[1]] = trunc_FLDS[1] - - # %% Apply to data - PREC[:, :, t] = np.reshape(e_corrSpatTemp_PREC, (dim_lon, dim_lat)) - TBOT[:, :, t] = np.reshape(e_corrSpatTemp_TBOT, (dim_lon, dim_lat)) - FSDS[:, :, t] = np.reshape(e_corrSpatTemp_FSDS, (dim_lon, dim_lat)) - FLDS[:, :, t] = np.reshape(e_corrSpatTemp_FLDS, (dim_lon, dim_lat)) - - del ( - e_corrSpatTemp_PREC, - e_corrSpatTemp_TBOT, - e_corrSpatTemp_FSDS, - e_corrSpatTemp_FLDS, - ) - - # %% Save disturbed variables - # Write PREC variable - - PRECTmms_ID[:] = np.transpose(PREC, (2, 1, 0)) - - # Write TBOT variable - TBOT_ID[:] = np.transpose(TBOT, (2, 1, 0)) - - # Write FSDS variable - FSDS_ID[:] = np.transpose(FSDS, (2, 1, 0)) - - # Write FLDS variable - FLDS_ID[:] = np.transpose(FLDS, (2, 1, 0)) - - # pdb.set_trace() - PREC_old = PREC - TBOT_old = TBOT - FSDS_old = FSDS - FLDS_old = FLDS - - del PREC, TBOT, FSDS, FLDS - - i = i + 1 From fbbc8b8cdd9f09f8b289bdeab366805b3cea06ea Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 2 Dec 2025 12:54:45 +0100 Subject: [PATCH 21/21] remove unneeded code --- src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 index 421b0bd07d..935e26c76a 100644 --- a/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/clm5/biogeophys/SoilStateInitTimeConstMod.F90 @@ -180,9 +180,6 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) integer :: begp, endp integer :: begc, endc integer :: begg, endg -#ifdef USE_PDAF - logical :: parameters_in_file, parameters_in_file_adj -#endif !----------------------------------------------------------------------- begp = bounds%begp; endp= bounds%endp