Read skin temperature from input netcdf file

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

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))

src/modsurface.f90

@@ -71,10 +71,11 @@ 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 modnetcdf,   only : readfile2d, readfile_dim_len,readfile1d
     use netcdf
 
     implicit none
@@ -728,38 +729,64 @@ 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)
-      !--- 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
-      allocate(tskininp(imax,jmax,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)
-    endif
+      call readfile1d('tskin.inp.'//cexpnr//'.nc','time',ttskin)
 
+      PRINT *, ttskin
+      allocate(tskininp(2:i1,2:j1,nttskin))
+
+      call readfile2d('tskin.inp.'//cexpnr//'.nc','tskin',tskininp)
+
+    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
 
     !$acc enter data copyin(z0m, z0h, obl, tskin, qskin, Cm, Cs, &
     !$acc&                  ustar, dudz, dvdz, thlflux, qtflux, &
     !$acc&                  dqtdz, dthldz, svflux, svs, horv, ra, rs, wsvsurf)
 
+    if (ltskininp) then 
+
+        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
 
@@ -872,17 +899,44 @@ 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,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
           tskin(i,j) = thls_patch(patchxnr(i), patchynr(j))
         end do
       end do
+
+    else if (ltskininp) then
+     !! 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
+
+         !    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
+
+
     else
       !$acc parallel loop collapse(2) default(present)
       do j = 2, j1
@@ -910,6 +964,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 +1427,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 +1570,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
@@ -2451,5 +2510,4 @@ subroutine do_lsm
     call qtsurf
 
   end subroutine do_lsm
-
 end module modsurface