From 0f0413cbfbc1780a242b00869944fc32159ba109 Mon Sep 17 00:00:00 2001 From: Chris Chapman Date: Fri, 7 Mar 2025 01:33:00 +1100 Subject: [PATCH 1/3] modification to modsurface to read spatially varying surface temperature --- src/modstartup.f90 | 27 ++++++++++++++++++++++----- src/modsurface.f90 | 40 +++++++++++++++++++++++++++++++++------- 2 files changed, 55 insertions(+), 12 deletions(-) diff --git a/src/modstartup.f90 b/src/modstartup.f90 index 981875c3f..21a274d61 100644 --- a/src/modstartup.f90 +++ b/src/modstartup.f90 @@ -364,6 +364,9 @@ subroutine startup(path) call initradiation call initchem call initsurface + PRINT *, thls + PRINT *, '=============' + call initdatetime call initemission call initlsm @@ -382,6 +385,7 @@ subroutine startup(path) call readinitfiles ! moved to obtain the correct btime for the timedependent forcings in case of a warmstart call inittimedep !depends on modglobal,modfields, modmpi, modsurf, modradiation call initpois ! hypre solver needs grid and baseprofiles + if(lopenbc) then ! Correct boundaries and initial field for divergence ! Create 1/int(rho) - must be after rhobf has been initialized allocate(rhointi(k1)) @@ -411,18 +415,21 @@ subroutine startup(path) if(myid==0) print *, 'Finished divergence correction initial field' endif + call inittstep call checkinitvalues call timer_toc('modstartup/startup') + PRINT *,'END start up' + PRINT *, 'thls: ', thls end subroutine startup !> Checks whether crucial parameters are set correctly subroutine checkinitvalues - use modsurfdata, only: wtsurf, wqsurf, ustin, thls, isurf, ps, lhetero + use modsurfdata, only: wtsurf, wqsurf, ustin, thls, isurf, ps, lhetero,ltskininp use modglobal, only: itot, jtot, ysize, xsize, dtmax, runtime, & startfile, lwarmstart, eps1, imax, jmax, ih, jh use modmpi, only: myid, nprocx, nprocy, mpierr, MPI_FINALIZE @@ -813,6 +820,7 @@ subroutine readinitfiles Wlm(:,:) = Wl(:,:) case(2) tskin(:,:) = thls + case(3,4) thls = thlprof(1) qts = qtprof(1) @@ -849,7 +857,10 @@ subroutine readinitfiles u0av(1) = uprof(1) thl0av(1) = thlprof(1) svs = svprof(1,:) - + + PRINT *, 'call baseprofs' + PRINT *, 'thls', thls + PRINT *, '=============' call baseprofs ! call baseprofs before thermodynamics #if defined(_OPENACC) @@ -1558,7 +1569,7 @@ subroutine baseprofs ! They are nevertheless calculated and printed to the stdin/baseprof files for user convenience use modfields, only : rhobf,rhobh,drhobdzf,drhobdzh use modglobal, only : k1,kmax,zf,zh,dzf,dzh,rv,rd,grav,cp,pref0,lwarmstart,ibas_prf,cexpnr,ifinput,ifoutput - use modsurfdata, only : thls,ps,qts + use modsurfdata, only : thls,ps,qts,ltskininp use modmpi, only : myid,comm3d,mpierr,D_MPI_BCAST implicit none @@ -1573,6 +1584,11 @@ subroutine baseprofs real,dimension(4) :: pmat real,dimension(4) :: tmat + PRINT *, '=====================' + PRINT *, 'base profile' + PRINT *, ibas_prf + PRINT *, '=====================' + allocate (height(k1),pb(k1),tb(k1)) if(myid==0)then @@ -1581,8 +1597,8 @@ subroutine baseprofs ibas_prf = 5 print *, 'WARNING: warm start requires input files for density. ibas_prf defaulted to 5' endif - if(ibas_prf <= 3 .and. thls < 0) then - STOP 'thls has not been initialized but is needed for setting up the base profiles.' + if(ibas_prf <= 3 .and. thls < 0) then + STOP 'thls has not been initialized but is needed for setting up the base profiles.' end if if(ibas_prf==1) then !thv constant and hydrostatic balance @@ -1616,6 +1632,7 @@ subroutine baseprofs enddo close(ifoutput) elseif(ibas_prf==3) then! use standard atmospheric lapse rate with surface temperature offset + tsurf=thls*(ps/pref0)**(rd/cp) pmat(1)=exp((log(ps)*lapserate(1)*rd+log(tsurf+zsurf*lapserate(1))*grav-& log(tsurf+zmat(1)*lapserate(1))*grav)/(lapserate(1)*rd)) diff --git a/src/modsurface.f90 b/src/modsurface.f90 index 49ffd374e..5225962fc 100644 --- a/src/modsurface.f90 +++ b/src/modsurface.f90 @@ -71,10 +71,10 @@ module modsurface contains !> Reads the namelists and initialises the soil. subroutine initsurface - use modglobal, only : i1, j1, i2, j2, itot, jtot,imax,jmax, nsv, ifnamopt, fname_options, ifinput, cexpnr, checknamelisterror, handle_err - use modraddata, only : iradiation,rad_shortw,irad_par,irad_user,irad_rrtmg,irad_rte_rrtmgp - use modmpi, only : myid, myidx, myidy, comm3d, mpierr, D_MPI_BCAST - use modtracers, only : tracer_prop + use modglobal, only : i1, j1, i2, j2, itot, jtot,imax,jmax, nsv, ifnamopt, fname_options, ifinput, cexpnr, checknamelisterror, handle_err + use modraddata, only : iradiation,rad_shortw,irad_par,irad_user,irad_rrtmg,irad_rte_rrtmgp + use modmpi, only : myid, myidx, myidy, comm3d, mpierr, D_MPI_BCAST + use modtracers, only : tracer_prop use netcdf implicit none @@ -728,6 +728,7 @@ subroutine initsurface endif if(ltskininp) then ! Use tskin.inp.iexpnr.nc to define surface skin temperature + PRINT *, "Reading skin temperature from file" !--- open tskin.inp.xxx.nc --- STATUS = NF90_OPEN('tskin.inp.'//cexpnr//'.nc', nf90_nowrite, NCID) if (STATUS .ne. nf90_noerr) call handle_err(STATUS) @@ -743,7 +744,7 @@ subroutine initsurface STATUS = NF90_GET_VAR (NCID, VARID, ttskin, start=(/1/), count=(/nttskin/) ) if (STATUS .ne. nf90_noerr) call handle_err(STATUS) !--- read tskin input - allocate(tskininp(imax,jmax,nttskin)) + allocate(tskininp(2:i1,2:j1,nttskin)) STATUS = NF90_INQ_VARID(NCID,'tskin', VARID) if (STATUS .ne. nf90_noerr) call handle_err(STATUS) STATUS = NF90_GET_VAR (NCID, VARID, tskininp, start=(/myidx*imax+1,myidy*jmax+1,1/), & @@ -751,8 +752,10 @@ subroutine initsurface if (STATUS .ne. nf90_noerr) call handle_err(STATUS) STATUS = NF90_CLOSE(NCID) if (STATUS .ne. nf90_noerr) call handle_err(STATUS) - endif + + endif + PRINT *, 'skin temp read in' dqtdz = 0 ! need to initialize, otherwise undefined in the first call to thermodynamics, before call surface (cold start) ustar = 0 ! need to initialize, otherwise undefined values in the corners in the first exchange @@ -760,6 +763,10 @@ subroutine initsurface !$acc& ustar, dudz, dvdz, thlflux, qtflux, & !$acc& dqtdz, dthldz, svflux, svs, horv, ra, rs, wsvsurf) + if (ltskininp.and.thls<0) then + + thls = sum(tskininp)/size(tskininp) + endif call timer_toc('modsurface/initsurface') end subroutine initsurface @@ -872,7 +879,8 @@ end subroutine calc_aerodynamic_resistance !> Prescribes the skin temperature subroutine presc_skin_temperature - use modglobal, only: i1, j1 + use modglobal, only: i1, j1,imax,jmax + use modmpi, only: myidx, myidy implicit none integer :: i, j @@ -883,6 +891,19 @@ subroutine presc_skin_temperature tskin(i,j) = thls_patch(patchxnr(i), patchynr(j)) end do end do + + else if (ltskininp) then + !! read tskin from input + !$acc parallel loop collapse(2) default(present) + + do j =2, j1 + do i=2, i1 + !!!!tskin(i,j) = tskininp(myidx*imax + i, myidy*jmax + j,1) !!!!! NOTE: the 1 is temporary. Add temporally evolving later + tskin(i,j) = tskininp(i,j,1) !!!!! NOTE: the 1 is temporary. Add temporally evolving later + end do + end do + + else !$acc parallel loop collapse(2) default(present) do j = 2, j1 @@ -910,6 +931,7 @@ subroutine calc_surface_scalars ! TODO: check if splitting these loops speeds things up on the GPU (async) !$acc parallel loop collapse(2) default(present) + do j = 2, j1 do i = 2, i1 tskin(i,j) = min(max(thlflux(i,j) / (Cs(i,j) * horv(i,j)), -10.), 10.) + thl0(i,j,1) @@ -1372,6 +1394,8 @@ subroutine getobl ! L is capped at 1e6 below, so use the same cap here L = 1e6 write(*,*) 'Obukhov length: Rib = 0 -> setting L=1e6' + PRINT *, 'stopping on line 1400' + STOP else iter = 0 L = obl(i,j) @@ -1513,6 +1537,8 @@ subroutine getobl ! L is capped at 1e6 below, so use the same cap here L = 1e6 write(*,*) 'Obukhov length: Rib = 0 -> setting L=1e6 (2nd point)' + PRINT *, 'Stopping on line 1543' + STOP else iter = 0 L = oblav From 34e834fdcb894812578474a80a33960782596a0c Mon Sep 17 00:00:00 2001 From: Chris Chapman Date: Fri, 28 Mar 2025 23:53:40 +1100 Subject: [PATCH 2/3] Added helper functions to modnetcdf --- src/modnetcdf.f90 | 99 +++++++++++++++++++++++++++++++++++++++++++++- src/modsurface.f90 | 70 ++++++++++++++++++++------------ 2 files changed, 142 insertions(+), 27 deletions(-) diff --git a/src/modnetcdf.f90 b/src/modnetcdf.f90 index 09f7a667e..d229b6e1b 100644 --- a/src/modnetcdf.f90 +++ b/src/modnetcdf.f90 @@ -25,9 +25,8 @@ module modnetcdf implicit none - public - public :: check + public :: check, readfile_dim_len, readfile1d, readfile2d contains @@ -52,4 +51,100 @@ subroutine check(status, file, line) end subroutine check + + subroutine readfile1d(file,varname,out1d) + + use modglobal, only : handle_err + character(*), intent(in) :: file + character(*), intent(in) :: varname + real, intent(inout) :: out1d(:) + + + integer :: NCID, STATUS, varID, var_len + var_len = size(out1d) + + + STATUS = NF90_OPEN(file, nf90_nowrite, NCID) + if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + + STATUS = NF90_INQ_DIMID(NCID, varname, varID) + if (STATUS .ne. nf90_noerr) call handle_err(status) + + STATUS = NF90_GET_VAR (NCID, varID, out1d, start=(/1/), count=(/var_len/) ) + if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + + + STATUS = NF90_CLOSE(NCID) + if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + end subroutine readfile1d + + + subroutine readfile_dim_len(file,dimname,dim_len) + + + use modglobal, only : handle_err + character(*), intent(in) :: file + character(*), intent(in) :: dimname + integer, intent(out) :: dim_len + + + integer :: NCID, STATUS, dimID + character(len = nf90_max_name) :: RecordDimName + + STATUS = NF90_OPEN(file, nf90_nowrite, NCID) + if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + + STATUS = NF90_INQ_DIMID(NCID, dimname, dimID) + if (STATUS .ne. nf90_noerr) call handle_err(status) + + STATUS = nf90_INQUIRE_DIMENSION(NCID, dimID, len=dim_len) + if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + + STATUS = NF90_CLOSE(NCID) + if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + + end subroutine readfile_dim_len + + subroutine readfile2d(file,varname,out2d) + + use modmpi, only : myidx, myidy + use modglobal, only : handle_err, i1, j1,imax,jmax + + character(*), intent(in) :: file + character(*), intent(in) :: varname + real, intent(inout) :: out2d(:,:,:) + + integer :: NCID, STATUS, timeID, nt, varID + character(len = nf90_max_name) :: RecordDimName + real, allocatable :: time(:) + + + STATUS = NF90_OPEN(file, nf90_nowrite, NCID) + if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + + STATUS = NF90_INQ_DIMID(NCID, "time", timeID) + if (STATUS .ne. nf90_noerr) call handle_err(status) + + STATUS = nf90_INQUIRE_DIMENSION(NCID, timeID, len=nt, name=RecordDimName) + if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + + + allocate(time(nt)) + STATUS = NF90_INQ_VARID(NCID, 'time', timeID) + if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + STATUS = NF90_GET_VAR (NCID, timeID, time, start=(/1/), count=(/nt/) ) + if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + PRINT *, time + + STATUS = NF90_INQ_VARID(NCID,'tskin', VARID) + if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + STATUS = NF90_GET_VAR (NCID, VARID, out2d, start=(/myidx*imax+1,myidy*jmax+1,1/), & + & count=(/imax,jmax,nt/)) + if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + STATUS = NF90_CLOSE(NCID) + if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + + + end subroutine readfile2d + end module modnetcdf diff --git a/src/modsurface.f90 b/src/modsurface.f90 index 5225962fc..e874ab225 100644 --- a/src/modsurface.f90 +++ b/src/modsurface.f90 @@ -75,6 +75,7 @@ subroutine initsurface use modraddata, only : iradiation,rad_shortw,irad_par,irad_user,irad_rrtmg,irad_rte_rrtmgp use modmpi, only : myid, myidx, myidy, comm3d, mpierr, D_MPI_BCAST use modtracers, only : tracer_prop + use modnetcdf, only : readfile2d, readfile_dim_len,readfile1d use netcdf implicit none @@ -728,31 +729,46 @@ subroutine initsurface endif if(ltskininp) then ! Use tskin.inp.iexpnr.nc to define surface skin temperature - PRINT *, "Reading skin temperature from file" + + + !PRINT *, "Reading skin temperature from file" !--- open tskin.inp.xxx.nc --- - STATUS = NF90_OPEN('tskin.inp.'//cexpnr//'.nc', nf90_nowrite, NCID) - if (STATUS .ne. nf90_noerr) call handle_err(STATUS) - !--- get time dimensions - STATUS = NF90_INQ_DIMID(NCID, "time", timeID) - if (STATUS /= nf90_noerr) call handle_err(status) - STATUS = nf90_INQUIRE_DIMENSION(NCID, timeID, len=nttskin, name=RecordDimName) - if (STATUS .ne. nf90_noerr) call handle_err(STATUS) - !--- read time + !STATUS = NF90_OPEN('tskin.inp.'//cexpnr//'.nc', nf90_nowrite, NCID) + !if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + !!--- get time dimensions + !STATUS = NF90_INQ_DIMID(NCID, "time", timeID) + !if (STATUS /= nf90_noerr) call handle_err(status) + !STATUS = nf90_INQUIRE_DIMENSION(NCID, timeID, len=nttskin, name=RecordDimName) + !if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + !!--- read time + !allocate(ttskin(nttskin)) + !STATUS = NF90_INQ_VARID(NCID, 'time', VARID) + !if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + !STATUS = NF90_GET_VAR (NCID, VARID, ttskin, start=(/1/), count=(/nttskin/) ) + !if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + !!--- read tskin input + !allocate(tskininp(2:i1,2:j1,nttskin)) + + !STATUS = NF90_INQ_VARID(NCID,'tskin', VARID) + !if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + !STATUS = NF90_GET_VAR (NCID, VARID, tskininp, start=(/myidx*imax+1,myidy*jmax+1,1/), & + ! & count=(/imax,jmax,nttskin/)) + !if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + !STATUS = NF90_CLOSE(NCID) + !if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + + call readfile_dim_len('tskin.inp.'//cexpnr//'.nc','time',nttskin) + PRINT *, nttskin + PRINT *, '-----------' + PRINT *, ttskin + allocate(ttskin(nttskin)) - STATUS = NF90_INQ_VARID(NCID, 'time', VARID) - if (STATUS .ne. nf90_noerr) call handle_err(STATUS) - STATUS = NF90_GET_VAR (NCID, VARID, ttskin, start=(/1/), count=(/nttskin/) ) - if (STATUS .ne. nf90_noerr) call handle_err(STATUS) - !--- read tskin input + call readfile1d('tskin.inp.'//cexpnr//'.nc','time',ttskin) + + PRINT *, ttskin allocate(tskininp(2:i1,2:j1,nttskin)) - STATUS = NF90_INQ_VARID(NCID,'tskin', VARID) - if (STATUS .ne. nf90_noerr) call handle_err(STATUS) - STATUS = NF90_GET_VAR (NCID, VARID, tskininp, start=(/myidx*imax+1,myidy*jmax+1,1/), & - & count=(/imax,jmax,nttskin/)) - if (STATUS .ne. nf90_noerr) call handle_err(STATUS) - STATUS = NF90_CLOSE(NCID) - if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + call readfile2d('tskin.inp.'//cexpnr//'.nc','tskin',tskininp) endif PRINT *, 'skin temp read in' @@ -763,9 +779,13 @@ subroutine initsurface !$acc& ustar, dudz, dvdz, thlflux, qtflux, & !$acc& dqtdz, dthldz, svflux, svs, horv, ra, rs, wsvsurf) - if (ltskininp.and.thls<0) then + if (ltskininp) then - thls = sum(tskininp)/size(tskininp) + thls = sum(tskininp(:,:,1))/size(tskininp(:,:,1)) + !PRINT *, sum(tskininp_test(:,:,1))/size(tskininp_test(:,:,1)) + PRINT *, 'thls: ', thls + !STOP + !SLEEP(5) endif call timer_toc('modsurface/initsurface') end subroutine initsurface @@ -899,7 +919,8 @@ subroutine presc_skin_temperature do j =2, j1 do i=2, i1 !!!!tskin(i,j) = tskininp(myidx*imax + i, myidy*jmax + j,1) !!!!! NOTE: the 1 is temporary. Add temporally evolving later - tskin(i,j) = tskininp(i,j,1) !!!!! NOTE: the 1 is temporary. Add temporally evolving later + tskin(i,j) = tskininp(i,j,1) !!!!! NOTE: the 1 is temporary. Add temporally evolving later + !PRINT *, tskin(i,j) end do end do @@ -2477,5 +2498,4 @@ subroutine do_lsm call qtsurf end subroutine do_lsm - end module modsurface From ea45666a300d5d39eb78be95e9db0e286279220f Mon Sep 17 00:00:00 2001 From: Chris Chapman Date: Thu, 3 Apr 2025 22:13:02 +1100 Subject: [PATCH 3/3] Updated pres_tskin to interpolate temporally --- src/modsurface.f90 | 28 ++++++++++++++++++++-------- 1 file changed, 20 insertions(+), 8 deletions(-) diff --git a/src/modsurface.f90 b/src/modsurface.f90 index e874ab225..6f0ad21e6 100644 --- a/src/modsurface.f90 +++ b/src/modsurface.f90 @@ -899,12 +899,12 @@ end subroutine calc_aerodynamic_resistance !> Prescribes the skin temperature subroutine presc_skin_temperature - use modglobal, only: i1, j1,imax,jmax + use modglobal, only: i1, j1,imax,jmax,rtimee use modmpi, only: myidx, myidy implicit none - integer :: i, j - + integer :: i, j, t + real :: fac if (lhetero) then do j = 2, j1 do i = 2, i1 @@ -913,13 +913,25 @@ subroutine presc_skin_temperature end do else if (ltskininp) then - !! read tskin from input - !$acc parallel loop collapse(2) default(present) - + !! skin temperature read from file + !$acc parallel loop collapse(2) default(present) + t=1 + do while(rtimee>ttskin(t+1)) + t=t+1 + end do + PRINT *, 'surface forcing timestep' + PRINT *, rtimee + PRINT *, ttskin(t) + PRINT *, ttskin(t+1) + PRINT *,'-------------------------' + + fac = ( rtimee-ttskin(t) ) / ( ttskin(t+1)-ttskin(t)) do j =2, j1 do i=2, i1 - !!!!tskin(i,j) = tskininp(myidx*imax + i, myidy*jmax + j,1) !!!!! NOTE: the 1 is temporary. Add temporally evolving later - tskin(i,j) = tskininp(i,j,1) !!!!! NOTE: the 1 is temporary. Add temporally evolving later + + ! wqsurf = wqsurft(t) + fac * ( wqsurft(t+1) - wqsurft(t) ) + + tskin(i,j) = tskininp(i,j,t) + fac * (tskininp(i,j,t+1) - tskininp(i,j,t)) !PRINT *, tskin(i,j) end do end do