diff --git a/cases/botany-dephy/backrad.inp.001.nc b/cases/botany-dephy/backrad.inp.001.nc new file mode 100644 index 00000000..c455b145 Binary files /dev/null and b/cases/botany-dephy/backrad.inp.001.nc differ diff --git a/cases/botany-dephy/convert_input.py b/cases/botany-dephy/convert_input.py new file mode 100644 index 00000000..844622e6 --- /dev/null +++ b/cases/botany-dephy/convert_input.py @@ -0,0 +1,77 @@ +"""Converts prof.inp and lscale.inp input files to new NetCDF-based input.""" +import argparse + +import numpy as np +import netCDF4 as nc + +PROF_NAMES = ["zh", "thetal", "qt", "ua", "va", "tke"] +LSCALE_NAMES = ["zh", "ug", "vg", "wa", "dqtdxls", "dqtdyls", "tnqt_adv", "tnthetal_rad"] + + +def parse_args(): + parser = argparse.ArgumentParser() + + parser.add_argument("-p", "--prof", type=str, default="prof.inp.001") + parser.add_argument("-l", "--lscale", type=str, default="lscale.inp.001") + parser.add_argument("-o", "--output", type=str, default="init.001.nc") + + return parser.parse_args() + + +def main(): + args = parse_args() + + prof = np.genfromtxt(args.prof, skip_header=2) + lscale = np.genfromtxt(args.lscale, skip_header=2) + + with nc.Dataset(args.output, "w") as ds: + height = ds.createDimension("zh", len(prof)) + heights = ds.createVariable("zh", "f", ("zh",)) + heights[:] = prof[:,0] + + for ivar in range(1, len(prof[0])): + nc_var = ds.createVariable(PROF_NAMES[ivar], "f", ("zh",)) + nc_var[:] = prof[:,ivar] + + for ivar in range(1, len(lscale[0])): + nc_var = ds.createVariable(LSCALE_NAMES[ivar], "f", ("zh",)) + nc_var[:] = lscale[:,ivar] + + # Nudging + dim = ds.createDimension("time") + var = ds.createVariable("time", "f", ("time",)) + var[0] = 0 + var[1] = 1.00000000E+07 + + a = 2 + b = 3 + c = 7.4 + nudging_constant = 6 * 3600 + (b * (0.5 * np.pi + np.arctan(a * 0.5 * np.pi * (1 - prof[:,0] / 3000))))**c + + var = ds.createVariable("ua_nud", "f", ("time", "zh")) + var[0,:] = prof[:,3] + var[1,:] = prof[:,3] + + var = ds.createVariable("nudging_constant_ua", "f", ("time", "zh")) + var[0,:] = nudging_constant + var[1,:] = nudging_constant + + var = ds.createVariable("thetal_nud", "f", ("time", "zh")) + var[0,:] = prof[:,1] + var[1,:] = prof[:,1] + + var = ds.createVariable("nudging_constant_thetal", "f", ("time", "zh")) + var[0,:] = nudging_constant + var[1,:] = nudging_constant + + var = ds.createVariable("qt_nud", "f", ("time", "zh")) + var[0,:] = prof[:,2] + var[1,:] = prof[:,2] + + var = ds.createVariable("nudging_constant_qt", "f", ("time", "zh")) + var[0,:] = nudging_constant + var[1,:] = nudging_constant + + +if __name__ == "__main__": + main() diff --git a/cases/botany-dephy/init.001.nc b/cases/botany-dephy/init.001.nc new file mode 100644 index 00000000..2643e285 Binary files /dev/null and b/cases/botany-dephy/init.001.nc differ diff --git a/cases/botany-dephy/namoptions.001 b/cases/botany-dephy/namoptions.001 new file mode 100644 index 00000000..8e69b9e4 --- /dev/null +++ b/cases/botany-dephy/namoptions.001 @@ -0,0 +1,184 @@ +&RUN +iexpnr = 001 +lwarmstart = .false. +startfile = 'initdlatestmx000y000.001' +runtime = 43200 +trestart = 43200 +ladaptive = .true. +irandom = 44 +randthl = 0.1 +randqt = 2.5e-5 +nsv = 0 +courant = .5 +peclet = .1 +loutdirs = .true. +/ + +&DOMAIN +itot = 1536 +jtot = 1536 +kmax = 175 + +xsize = 153600 +ysize = 153600 + +xlat = 13.1 +xlon = 302 +xyear = 2020 ! Currently only used in netCDF time unit +xday = 32 +xtime = 0.0 +/ + +&PHYSICS +z0 = 0.00016 ! 1.6e-4 +ustin = 0.32 +ps = 101604.87 +thls = 298.5 +! thls = 299.135 ! EURECA +! thls = 298.5 ! RICO + +lmoist = .true. +isurf = 2 !Forced surface temperature; fluxes are calculated +timerad = 60 +iradiation = 4 +lcoriol = .true. +lmomsubs = .false. +/ +&NAMSURFACE +thls = 298.5 ! duplicate temperature ! +z0mav = 1.6e-4 +z0hav = 3.2e-5 +albedoav = 0.07 ! From RCEMIP. note albedoav is used if lCnstAlbedo = .true. in NAMRADIATION +/ +&NAMRADIATION +lCnstAlbedo = .true. ! use a albedoav as albedo instead of angle-dependent parameterization +! locean = .true. ! use ocean parameterization of albedo, otherwise land. (only for lCnstAlbedo = .false.) +/ + +&DYNAMICS +llsadv = .false. +lqlnr = .false. ! obsolete +cu = 0. ! -5 in RICO +cv = 0. + +iadv_mom = 62 +iadv_tke = 62 +iadv_thl = 62 +iadv_qt = 52 +iadv_sv = 52 +/ + +&NAMMICROPHYSICS +imicro = 2 +Nc_0 = 70000000.0 +l_sb = .true. +/ + +&NAMBULKMICROSTAT +lmicrostat = .true. +timeav = 300 +dtav = 60 +/ +&NAMCHECKSIM +tcheck = 0 +/ +&NAMSAMPLING +lsampcl = .false. +lsampco = .false. +dtav = 60 +timeav = 21600 +/ +&NAMTIMESTAT +ltimestat = .true. +dtav = 60 +/ +&NAMCROSSSECTION +lcross = .true. +lxz = .false. +lyz = .false. +crossheight = 1, 13, 35, 52 ! 7.5m, 198m, 611m, 999m +dtav = 300 +/ + +&NAMGENSTAT +lstat = .true. +dtav = 300 +timeav = 300 +/ +&NAMCAPE +lcape = .true. +dtav = 300 +/ +&NAMFIELDDUMP +lfielddump = .false. +! lclassic = .false. +dtav = 3600 +khigh = 149 ! 132=4030m 149=5048m +! tmin = 43200 ! start fielddumps at 12h (implemented in to4.4_Fredrik) +! select vartiables to store (implemented in to4.4_Fredrik, not in v4.3) +lu = .true. ! consider omitting u, v +lv = .true. +lw = .true. +lqt = .true. +lql = .true. +lthl = .true. +lbuoy = .false. +lsv = .false., .true. ! not Nr, just qr +/ +&NAMSTATTEND +dtav = 60 +timeav = 300 +ltend = .false. +/ +&NAMRADSTAT +lstat = .true. +dtav = 60 +timeav = 300 +/ + +&NAMRADFIELD +lradfield = .true. +dtav = 3600 +timeav = 3600 +/ +&NAMVARBUDGET +lvarbudget = .true. ! Note requires interactive radiation or it breaks +dtav = 300 +timeav = 300 +/ +&NAMCLOUDFIELD +dtav=1440 +lcloudfield = .false. +/ +&NAMNETCDFSTATS +lnetcdf = .true. +! lsync = .true. +lclassic = .true. ! NetCDF 3. No compression but saves RAM. +/ +&SOLVER +solver_id = 100 ! FFTW +/ + +&NAMNUDGE +! nudge to profiles given in nudge.inp.NNN +lnudge = .true. +tnudgefac = 1. +/ + +&VVUQ_extra +dudz = 0.0022 +u0 = -10 +thl_low = +w0 = 0.0045 +wpamp = -0.00085 +case = 2 +thl_Gamma = 5.0 +qt_lambda = 1850 +qt0 = 0.01425 +z_ml = 500 +thl_tend0 = -5.78e-06 +qt_tend0 = -1.73e-08 +dthl0 = 1.25 +!thl_tend_z_max = 2000 +!qt_tend_z_max = 4000 +/ diff --git a/cases/botany-dephy/rrtmg_lw.nc b/cases/botany-dephy/rrtmg_lw.nc new file mode 100644 index 00000000..57d0d77c Binary files /dev/null and b/cases/botany-dephy/rrtmg_lw.nc differ diff --git a/cases/botany-dephy/rrtmg_sw.nc b/cases/botany-dephy/rrtmg_sw.nc new file mode 100644 index 00000000..5db12c4f Binary files /dev/null and b/cases/botany-dephy/rrtmg_sw.nc differ diff --git a/src/modnudge.f90 b/src/modnudge.f90 index a39f51ef..472ad0af 100644 --- a/src/modnudge.f90 +++ b/src/modnudge.f90 @@ -1,14 +1,5 @@ !> \file modnudge.f90 -!! Nudges theta_l and q_t profiles to the initial profiles on a timescale tnudgeT -!> - -!> -!! Nudges theta_l and q_t profiles to the initial profiles on a timescale tnudgeT -!> -!! \author Thijs Heus,MPI-M -!! \par Revision list -!! \todo Documentation -! This file is part of DALES. +! This file is part of DALES. ! ! DALES is free software; you can redistribute it and/or modify ! it under the terms of the GNU General Public License as published by @@ -23,168 +14,411 @@ ! You should have received a copy of the GNU General Public License ! along with this program. If not, see . ! -! Copyright 1993-2009 Delft University of Technology, Wageningen University, Utrecht University, KNMI +! Copyright 1993-2024 The DALES team. ! +!> Module for nudging prognostic fields to some provided profiles. +module modnudge + use modprecision, only: field_r + use modtimer, only: timer_tic, timer_toc + implicit none + private -module modnudge -use modprecision, only: field_r - -implicit none -PRIVATE -PUBLIC :: initnudge, nudge,exitnudge -SAVE - real(field_r), dimension(:,:), allocatable :: tnudge,unudge,vnudge,wnudge,thlnudge,qtnudge - real, dimension(:) , allocatable :: timenudge - real :: tnudgefac = 1. - logical :: lnudge = .false.,lunudge,lvnudge,lwnudge,lthlnudge,lqtnudge - integer :: ntnudge = 10000 + save + + character(*), parameter :: modname = "modnudge" + + ! Switches for enabling/disabling nudging + logical :: lnudge = .false. + logical :: lunudge = .true. + logical :: lvnudge = .true. + logical :: lwnudge = .true. + logical :: lthlnudge = .true. + logical :: lqtnudge = .true. + + ! Nudging profiles + real(field_r), allocatable :: tnudge(:,:) + real(field_r), allocatable :: unudge(:,:) + real(field_r), allocatable :: vnudge(:,:) + real(field_r), allocatable :: wnudge(:,:) + real(field_r), allocatable :: thlnudge(:,:) + real(field_r), allocatable :: qtnudge(:,:) + ! Nudging constants (only supported with DEPHY input) + real(field_r), allocatable :: tunudge(:,:) + real(field_r), allocatable :: tvnudge(:,:) + real(field_r), allocatable :: twnudge(:,:) + real(field_r), allocatable :: tthlnudge(:,:) + real(field_r), allocatable :: tqtnudge(:,:) + + real(field_r), allocatable :: timenudge(:) + + ! Nudging timescale + real(field_r) :: tnudgefac = 1. + ! Number of nudging time steps + integer :: ntnudge = 10000 + + public :: initnudge + public :: nudge + public :: exitnudge contains + + !> Read nudging profiles and nudging constants from input, and allocate + !! memory. subroutine initnudge - use modmpi, only :myid,mpierr,comm3d,D_MPI_BCAST - use modglobal,only :ifnamopt,fname_options,runtime,cexpnr,ifinput,k1,kmax,checknamelisterror - implicit none + use modmpi, only: myid, mpierr, comm3d, D_MPI_BCAST + use modglobal, only: ifnamopt, fname_options, runtime, cexpnr, ifinput, & + k1, kmax, checknamelisterror, lstart_netcdf + use modstat_nc - integer :: ierr,k,t - real,allocatable,dimension(:) :: height + character(*), parameter :: routine = modname//"::initnudge" + + integer :: ierr, k, t + integer :: ncid, varid, dimid character(1) :: chmess1 - namelist /NAMNUDGE/ & - lnudge,tnudgefac - allocate(tnudge(k1,ntnudge),unudge(k1,ntnudge),vnudge(k1,ntnudge),wnudge(k1,ntnudge),thlnudge(k1,ntnudge),qtnudge(k1,ntnudge)) - allocate(timenudge(0:ntnudge), height(k1)) - tnudge = 0 - unudge=0 - vnudge=0 - wnudge=0 - thlnudge=0 - qtnudge=0 - timenudge=0 - - if(myid==0)then - - open(ifnamopt,file=fname_options,status='old',iostat=ierr) - read (ifnamopt,NAMNUDGE,iostat=ierr) + real, allocatable, dimension(:) :: height + + namelist /NAMNUDGE/ lnudge, lunudge, lvnudge, lwnudge, lthlnudge, & + lqtnudge, tnudgefac + + if (myid == 0) then + open(ifnamopt, file=fname_options, status='old', iostat=ierr) + read(ifnamopt, NAMNUDGE, iostat=ierr) call checknamelisterror(ierr, ifnamopt, 'NAMNUDGE') - write(6 ,NAMNUDGE) + write(6, NAMNUDGE) close(ifnamopt) end if - call D_MPI_BCAST(lnudge , 1,0,comm3d,mpierr) + + call D_MPI_BCAST(lnudge, 1, 0, comm3d, mpierr) + call D_MPI_BCAST(lunudge, 1, 0, comm3d, mpierr) + call D_MPI_BCAST(lvnudge, 1, 0, comm3d, mpierr) + call D_MPI_BCAST(lwnudge, 1, 0, comm3d, mpierr) + call D_MPI_BCAST(lthlnudge, 1, 0, comm3d, mpierr) + call D_MPI_BCAST(lqtnudge, 1, 0, comm3d, mpierr) + call D_MPI_BCAST(tnudgefac, 1, 0, comm3d, mpierr) if (.not. lnudge) return - if(myid==0) then - t = 0 - open (ifinput,file='nudge.inp.'//cexpnr) - do while (timenudge(t) < runtime) - t = t + 1 - if (t > ntnudge) then - write (*,*) "Too many time points in file ", 'nudge.inp.'//cexpnr, ", the limit is ntnudge = ", ntnudge - stop + call timer_tic(routine, 0) + + if (lstart_netcdf) then + if (myid == 0) then + call nchandle_error(nf90_open("init."//cexpnr//".nc", NF90_NOWRITE, & + ncid)) + + ! Get the number of forcing time steps + call nchandle_error(nf90_inq_dimid(ncid, "time", dimid)) + call nchandle_error(nf90_inquire_dimension(ncid, dimid, len=ntnudge)) + end if + + call D_MPI_BCAST(ntnudge, 1, 0, comm3d, mpierr) + allocate(timenudge(0:ntnudge)) + timenudge(:) = 0 + if (myid == 0) then + call nchandle_error(nf90_inq_varid(ncid, "time", varid)) + call nchandle_error(nf90_get_var(ncid, varid, timenudge(1:))) + end if + call D_MPI_BCAST(timenudge, ntnudge + 1, 0, comm3d, mpierr) + + if (lunudge) then + allocate(unudge(k1,ntnudge), tunudge(k1,ntnudge)) + if (myid == 0) then + call nchandle_error(nf90_inq_varid(ncid, "ua_nud", varid)) + call nchandle_error(nf90_get_var(ncid, varid, unudge(1:kmax,:))) + call nchandle_error(nf90_inq_varid(ncid, "nudging_constant_ua", & + varid)) + call nchandle_error(nf90_get_var(ncid, varid, tunudge(1:kmax,:))) + end if + end if + + if (lvnudge) then + allocate(vnudge(k1,ntnudge), tvnudge(k1,ntnudge)) + if (myid == 0) then + call nchandle_error(nf90_inq_varid(ncid, "va_nud", varid)) + call nchandle_error(nf90_get_var(ncid, varid, vnudge(1:kmax,:))) + call nchandle_error(nf90_inq_varid(ncid, "nudging_constant_va", & + varid)) + call nchandle_error(nf90_get_var(ncid, varid, tvnudge(1:kmax,:))) + end if + end if + + if (lwnudge) then + allocate(wnudge(k1,ntnudge), twnudge(k1,ntnudge)) + if (myid == 0) then + call nchandle_error(nf90_inq_varid(ncid, "wa_nud", varid)) + call nchandle_error(nf90_get_var(ncid, varid, wnudge)) + call nchandle_error(nf90_inq_varid(ncid, "nudging_constant_wa", & + varid)) + call nchandle_error(nf90_get_var(ncid, varid, twnudge)) + end if + end if + + if (lthlnudge) then + allocate(thlnudge(k1,ntnudge), tthlnudge(k1,ntnudge)) + if (myid == 0) then + call nchandle_error(nf90_inq_varid(ncid, "thetal_nud", varid)) + call nchandle_error(nf90_get_var(ncid, varid, thlnudge(1:kmax,:))) + call nchandle_error(nf90_inq_varid(ncid, "nudging_constant_thetal", & + varid)) + call nchandle_error(nf90_get_var(ncid, varid, tthlnudge(1:kmax,:))) + end if + end if + + if (lqtnudge) then + allocate(qtnudge(k1,ntnudge), tqtnudge(k1,ntnudge)) + if (myid == 0) then + call nchandle_error(nf90_inq_varid(ncid, "qt_nud", varid)) + call nchandle_error(nf90_get_var(ncid, varid, qtnudge(1:kmax,:))) + call nchandle_error(nf90_inq_varid(ncid, "nudging_constant_qt", & + varid)) + call nchandle_error(nf90_get_var(ncid, varid, tqtnudge(1:kmax,:))) end if + end if + + if (myid == 0) then + call nchandle_error(nf90_close(ncid)) + end if + else + allocate(tnudge(k1,ntnudge), unudge(k1,ntnudge), vnudge(k1,ntnudge), & + wnudge(k1,ntnudge), thlnudge(k1,ntnudge), qtnudge(k1,ntnudge)) + allocate(tunudge(k1,ntnudge), tvnudge(k1,ntnudge), & + twnudge(k1,ntnudge), tthlnudge(k1,ntnudge), & + tqtnudge(k1,ntnudge)) + allocate(timenudge(0:ntnudge), height(k1)) + + tnudge = 0 + unudge = 0 + vnudge = 0 + wnudge = 0 + thlnudge = 0 + qtnudge = 0 + timenudge = 0 + t = 0 + if (myid == 0) then + open(ifinput, file='nudge.inp.'//cexpnr) - chmess1 = "#" - ierr = 1 ! not zero - !search for the next line consisting of "# time", from there onwards the profiles will be read - do while (.not.(chmess1 == "#" .and. ierr ==0)) - read(ifinput,*,iostat=ierr) chmess1, timenudge(t) - if (ierr < 0) then - stop 'STOP: No time dependend nudging data for end of run' + do while (timenudge(t) < runtime) + t = t + 1 + if (t > ntnudge) then + write(*, *) "Too many time points in file ", 'nudge.inp.'//cexpnr, & + &", the limit is ntnudge = ", ntnudge + stop end if - end do - write(6,*) 'time', timenudge(t) - write(6,*) ' height t_nudge u_nudge v_nudge w_nudge thl_nudge qt_nudge' - do k=1,kmax - read (ifinput,*) & - height (k), & - tnudge (k,t), & - unudge (k,t), & - vnudge (k,t), & - wnudge (k,t), & - thlnudge(k,t), & - qtnudge(k,t) - end do + chmess1 = "#" + ierr = 1 ! not zero + !search for the next line consisting of "# time", from there onwards the profiles will be read + do while (.not. (chmess1 == "#" .and. ierr == 0)) + read(ifinput, *, iostat=ierr) chmess1, timenudge(t) + if (ierr < 0) then + stop 'STOP: No time dependend nudging data for end of run' + end if - do k=kmax,1,-1 - write (6,'(f7.1,6e12.4)') & - height (k), & - tnudge (k,t), & - unudge (k,t), & - vnudge (k,t), & - wnudge (k,t), & - thlnudge(k,t), & - qtnudge(k,t) + end do + write(6, *) 'time', timenudge(t) + write(6, *) ' height t_nudge u_nudge v_nudge w_nudge & + &thl_nudge qt_nudge' + do k = 1, kmax + read(ifinput, *) & + height(k), & + tnudge(k,t), & + unudge(k,t), & + vnudge(k,t), & + wnudge(k,t), & + thlnudge(k,t), & + qtnudge(k,t) + end do + + do k = kmax, 1, -1 + write(6, '(f7.1,6e12.4)') & + height(k), & + tnudge(k,t), & + unudge(k,t), & + vnudge(k,t), & + wnudge(k,t), & + thlnudge(k,t), & + qtnudge(k,t) + end do end do - end do - close(ifinput) - tnudge = tnudgefac*tnudge - end if - call D_MPI_BCAST(timenudge ,ntnudge+1 ,0,comm3d,mpierr) - call D_MPI_BCAST(tnudge ,k1*ntnudge,0,comm3d,mpierr) - call D_MPI_BCAST(unudge ,k1*ntnudge,0,comm3d,mpierr) - call D_MPI_BCAST(vnudge ,k1*ntnudge,0,comm3d,mpierr) - call D_MPI_BCAST(wnudge ,k1*ntnudge,0,comm3d,mpierr) - call D_MPI_BCAST(thlnudge ,k1*ntnudge,0,comm3d,mpierr) - call D_MPI_BCAST(qtnudge ,k1*ntnudge,0,comm3d,mpierr) - lunudge = any(abs(unudge)>1e-8) - lvnudge = any(abs(vnudge)>1e-8) - lwnudge = any(abs(wnudge)>1e-8) - lthlnudge = any(abs(thlnudge)>1e-8) - lqtnudge = any(abs(qtnudge)>1e-8) - - !$acc enter data copyin(timenudge, tnudge, unudge, vnudge, wnudge, thlnudge, qtnudge) - end subroutine initnudge + close (ifinput) + end if + + lunudge = any(abs(unudge) > 1e-8) + lvnudge = any(abs(vnudge) > 1e-8) + lwnudge = any(abs(wnudge) > 1e-8) + lthlnudge = any(abs(thlnudge) > 1e-8) + lqtnudge = any(abs(qtnudge) > 1e-8) + + tnudge = tnudgefac * tnudge + + tunudge(:,:) = tnudge(:,:) + tvnudge(:,:) = tnudge(:,:) + twnudge(:,:) = tnudge(:,:) + tthlnudge(:,:) = tnudge(:,:) + tqtnudge(:,:) = tnudge(:,:) + end if + + + call D_MPI_BCAST(timenudge, ntnudge + 1, 0, comm3d, mpierr) + if (lunudge) then + call D_MPI_BCAST(unudge, k1 * ntnudge, 0, comm3d, mpierr) + call D_MPI_BCAST(tunudge, k1 * ntnudge, 0, comm3d, mpierr) + end if + if (lvnudge) then + call D_MPI_BCAST(vnudge, k1 * ntnudge, 0, comm3d, mpierr) + call D_MPI_BCAST(tvnudge, k1 * ntnudge, 0, comm3d, mpierr) + end if + if (lwnudge) then + call D_MPI_BCAST(wnudge, k1 * ntnudge, 0, comm3d, mpierr) + call D_MPI_BCAST(twnudge, k1 * ntnudge, 0, comm3d, mpierr) + end if + if (lthlnudge) then + call D_MPI_BCAST(thlnudge, k1 * ntnudge, 0, comm3d, mpierr) + call D_MPI_BCAST(tthlnudge, k1 * ntnudge, 0, comm3d, mpierr) + end if + if (lqtnudge) then + call D_MPI_BCAST(qtnudge, k1 * ntnudge, 0, comm3d, mpierr) + call D_MPI_BCAST(tqtnudge, k1 * ntnudge, 0, comm3d, mpierr) + end if -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !$acc enter data copyin(timenudge, unudge, vnudge, wnudge, thlnudge, & + !$acc& qtnudge, tunudge, tvnudge, twnudge, tthlnudge, & + !$acc& tqtnudge) + call timer_toc(routine) + end subroutine initnudge + + !> Perform nudging of velocities, temperature and humidity fields. subroutine nudge - use modglobal, only : timee,rtimee,i1,j1,kmax,rdt - use modfields, only : up,vp,wp,thlp, qtp,u0av,v0av,qt0av,thl0av - implicit none + use modglobal, only: timee, rtimee, i1, j1, kmax, rdt + use modfields, only: up, vp, wp, thlp, qtp, u0av, v0av, qt0av, thl0av + + character(*), parameter :: routine = modname//"::nudge" - integer i,j,k,t - real :: dtm,dtp,currtnudge + integer :: i, j, k, t + real(field_r) :: dtm, dtp, currtnudge - if (.not.(lnudge)) return -! if (rk3step/=3) return - if (timee==0) return + if (.not. (lnudge)) return - t=1 - do while(rtimee>timenudge(t)) - t=t+1 + if (timee == 0) return + + call timer_tic(routine, 0) + + t = 1 + do while (rtimee > timenudge(t)) + t = t + 1 end do - if (rtimee>timenudge(1)) then - t=t-1 - end if - - dtm = ( rtimee-timenudge(t) ) / ( timenudge(t+1)-timenudge(t) ) - dtp = ( timenudge(t+1)-rtimee)/ ( timenudge(t+1)-timenudge(t) ) - - !$acc parallel loop collapse(3) default(present) - do k=1,kmax - do j = 2,j1 - do i = 2,i1 - currtnudge = max(1.0*rdt,tnudge(k,t)*dtp+tnudge(k,t+1)*dtm) ! 1.0: convert to double for nvfortran - if(lunudge ) up (i,j,k)=up (i,j,k)-& - (u0av (k)-(unudge (k,t)*dtp+unudge (k,t+1)*dtm))/currtnudge - if(lvnudge ) vp (i,j,k)=vp (i,j,k)-& - (v0av (k)-(vnudge (k,t)*dtp+vnudge (k,t+1)*dtm))/currtnudge - if(lwnudge ) wp (i,j,k)=wp (i,j,k)-& - ( -(wnudge (k,t)*dtp+wnudge (k,t+1)*dtm))/currtnudge - if(lthlnudge) thlp(i,j,k)=thlp(i,j,k)-& - (thl0av(k)-(thlnudge(k,t)*dtp+thlnudge(k,t+1)*dtm))/currtnudge - if(lqtnudge ) qtp (i,j,k)=qtp (i,j,k)-& - (qt0av (k)-(qtnudge (k,t)*dtp+qtnudge (k,t+1)*dtm))/currtnudge + if (rtimee > timenudge(1)) then + t = t - 1 + end if + + dtm = (rtimee - timenudge(t)) / (timenudge(t + 1) - timenudge(t)) + dtp = (timenudge(t + 1) - rtimee) / (timenudge(t + 1) - timenudge(t)) + + if (lunudge) then + !$acc parallel loop collapse(3) private(currtnudge) default(present) async + do k = 1, kmax + do j = 2, j1 + do i = 2, i1 + currtnudge = max(1.0_field_r * rdt, & + tunudge(k,t) * dtp + tunudge(k,t + 1) * dtm) + up(i,j,k) = up(i,j,k) - (u0av(k) - (unudge(k,t) * dtp + & + unudge(k,t + 1) * dtm)) / currtnudge end do - end do - end do + end do + end do + end if + + if (lvnudge) then + !$acc parallel loop collapse(3) default(present) private(currtnudge) async + do k = 1, kmax + do j = 2, j1 + do i = 2, i1 + currtnudge = max(1.0_field_r * rdt, & + tvnudge(k,t) * dtp + tvnudge(k,t + 1) * dtm) + vp(i,j,k) = vp(i,j,k) - (v0av(k) - (vnudge(k,t) * dtp + & + vnudge(k,t + 1) * dtm)) / currtnudge + end do + end do + end do + end if + + if (lwnudge) then + !$acc parallel loop collapse(3) default(present) private(currtnudge) async + do k = 1, kmax + do j = 2, j1 + do i = 2, i1 + currtnudge = max(1.0_field_r * rdt, & + twnudge(k,t) * dtp + twnudge(k,t + 1) * dtm) + wp(i,j,k) = wp(i,j,k) - ((wnudge(k,t) * dtp + wnudge(k,t + 1) & + * dtm)) / currtnudge + end do + end do + end do + end if + + if (lthlnudge) then + !$acc parallel loop collapse(3) default(present) private(currtnudge) async + do k = 1, kmax + do j = 2, j1 + do i = 2, i1 + currtnudge = max(1.0_field_r * rdt, & + tthlnudge(k,t) * dtp + tthlnudge(k,t + 1) * dtm) + thlp(i,j,k) = thlp(i,j,k) - (thl0av(k) - (thlnudge(k,t) * dtp + & + thlnudge(k,t + 1) * dtm)) / currtnudge + end do + end do + end do + end if + + if (lqtnudge) then + !$acc parallel loop collapse(3) default(present) private(currtnudge) async + do k = 1, kmax + do j = 2, j1 + do i = 2, i1 + currtnudge = max(1.0_field_r * rdt, & + tqtnudge(k,t) * dtp + tqtnudge(k,t + 1) * dtm) + qtp(i,j,k) = qtp(i,j,k) - (qt0av(k) - (qtnudge(k,t) * dtp + & + qtnudge(k,t + 1) * dtm)) / currtnudge + end do + end do + end do + end if + + !$acc wait + + call timer_toc(routine) end subroutine nudge + !> Deallocates memory. subroutine exitnudge - deallocate(timenudge) + if (.not. lnudge) return + + deallocate(timenudge) + + if (allocated(tnudge)) deallocate(tnudge) + + if (lunudge) then + deallocate(unudge, tunudge) + end if + + if (lvnudge) then + deallocate(vnudge, tvnudge) + end if + + if (lwnudge) then + deallocate(wnudge, twnudge) + end if + + if (lthlnudge) then + deallocate(thlnudge, tthlnudge) + end if + + if (lqtnudge) then + deallocate(qtnudge, tqtnudge) + end if end subroutine exitnudge end module