diff --git a/config/cheops_2018.cmake b/config/cheops_2018.cmake new file mode 100644 index 00000000..5397387b --- /dev/null +++ b/config/cheops_2018.cmake @@ -0,0 +1,26 @@ +# CHEOPS GCC - v. 2018 +# use with: +# module load gnu/4.8.2 netcdf/4.1.3-gcc-4.8.2 intelmpi/2018 +# module load hdf5/1.8.11 szlib +# module load cmake + +set(CMAKE_Fortran_COMPILER "gfortran" ) # "ifort") #"gfortran") +set(Fortran_COMPILER_WRAPPER mpif90 ) # mpiifort) + +set(USER_Fortran_FLAGS "-fbacktrace -finit-real=nan -fdefault-real-8 -fno-f2c -ffree-line-length-none") +set(USER_Fortran_FLAGS_RELEASE "-funroll-all-loops -O3") +set(USER_Fortran_FLAGS_DEBUG "-W -Wall -Wuninitialized -fcheck=all -fbacktrace -O0 -g -ffpe-trap=invalid,zero,overflow") + +set(NETCDF_INCLUDE_DIR "/opt/rrzk/lib/netcdf/4.1.3-gcc-4.8.2/include") # set(NETCDF_INCLUDE_DIR "/opt/rrzk/lib/netcdf/4.1.3/include") +set(NETCDF_LIB_1 "/opt/rrzk/lib/netcdf/4.1.3-gcc-4.8.2/lib/libnetcdff.so") # set(NETCDF_LIB_1 "/opt/rrzk/lib/netcdf/4.1.3/lib/libnetcdff.so") +set(NETCDF_LIB_2 "/opt/rrzk/lib/netcdf/4.1.3-gcc-4.8.2/lib/libnetcdf.so") # set(NETCDF_LIB_2 "/opt/rrzk/lib/netcdf/4.1.3/lib/libnetcdf.so") +# set(NETCDF_LIB_3 "/opt/rrzk/lib/netcdf/netcdf-4.1.3/f90/netcdf.mod") +set(HDF5_LIB_1 "/opt/rrzk/lib/hdf5/1.8.11/lib/libhdf5_hl.so") # set(HDF5_LIB_1 "/opt/rrzk/lib/hdf5/1.8.11/lib/libhdf5_hl.so") +set(HDF5_LIB_2 "/opt/rrzk/lib/hdf5/1.8.11/lib/libhdf5.so") # set(HDF5_LIB_2 "/opt/rrzk/lib/hdf5/1.8.11/lib/libhdf5.so") +set(SZIP_LIB "/opt/rrzk/lib/szip/szip-2.1/lib/libsz.so") +set(LIBS ${NETCDF_LIB_1} ${NETCDF_LIB_2} ${HDF5_LIB_1} ${HDF5_LIB_2} ${SZIP_LIB} m z curl) +# set(LIBS ${NETCDF_LIB_1} ${NETCDF_LIB_2} ${NETCDF_LIB_3} ${HDF5_LIB_1} ${HDF5_LIB_2} ${SZIP_LIB} m z curl) +# + + + diff --git a/config/cheops_CMakeLists.txt b/config/cheops_CMakeLists.txt new file mode 100644 index 00000000..64073d71 --- /dev/null +++ b/config/cheops_CMakeLists.txt @@ -0,0 +1,131 @@ +### Choose CMAKE Type +if(NOT CMAKE_BUILD_TYPE) + set (CMAKE_BUILD_TYPE RELEASE CACHE STRING + "Choose the type of build, options are: None Debug Release." + FORCE) +endif() + +### Set compiler flags +if("$ENV{SYST}" STREQUAL "HUYGENS") + set(CMAKE_Fortran_COMPILER "mpfort") + set(CMAKE_Fortran_FLAGS "-qfree=F90 -qrealsize=8 -qwarn64 -qflttrap=en:ov:zero:inv:imp -qflag=w:e" CACHE STRING "") + set(CMAKE_Fortran_FLAGS_RELEASE "-O4 -qnoipa -qstrict=none:exceptions" CACHE STRING "") + set(CMAKE_Fortran_FLAGS_DEBUG "-O2 -g -qfullpath -C -qflttrp=enable:nanq:overflow:zerodivide -qsigtrap -qinitauto=ff" CACHE STRING "") +elseif("$ENV{SYST}" STREQUAL "CARTESIUS") + set(CMAKE_Fortran_COMPILER "mpiifort") + set(CMAKE_Fortran_FLAGS "-r8 -ftz -extend_source" CACHE STRING "") + set(CMAKE_Fortran_FLAGS_RELEASE "-O3" CACHE STRING "") + set(CMAKE_Fortran_FLAGS_DEBUG "-traceback -fpe1 -O0 -g -check all" CACHE STRING "") +elseif("$ENV{SYST}" STREQUAL "localpc_ifort") + set(CMAKE_Fortran_COMPILER "mpif90") + set(CMAKE_Fortran_FLAGS "-r8 -ftz -extend_source" CACHE STRING "") + set(CMAKE_Fortran_FLAGS_RELEASE "-O3" CACHE STRING "") + set(CMAKE_Fortran_FLAGS_DEBUG "-traceback -fpe1 -O0 -g -check all" CACHE STRING "") +elseif("$ENV{SYST}" STREQUAL "HYDRA") + set(CMAKE_Fortran_COMPILER "mpiifort") + set(CMAKE_Fortran_FLAGS "-r8 -ftz -extend_source" CACHE STRING "") + set(CMAKE_Fortran_FLAGS_RELEASE "-O3" CACHE STRING "") + set(CMAKE_Fortran_FLAGS_DEBUG "-traceback -fpe1 -O0 -g -check all" CACHE STRING "") +else() + set(CMAKE_Fortran_COMPILER "mpif90") + set(CMAKE_Fortran_FLAGS "-finit-real=nan -W -Wall -fdefault-real-8 -ffree-line-length-none" CACHE STRING "") + set (CMAKE_Fortran_FLAGS_RELEASE "-funroll-all-loops -fno-f2c -O3" CACHE STRING "") + set (CMAKE_Fortran_FLAGS_DEBUG "-fbounds-check -fbacktrace -fno-f2c -O0 -g -ffpe-trap=invalid,zero,overflow" CACHE STRING "") +endif() + +## Project parameters +PROJECT(DALES Fortran) +cmake_minimum_required(VERSION 2.6) +set(VERSION_MAJOR "4") +set(VERSION_MINOR "1") +set(VERSION_PATCH "0") + +### If necessary, resort to BASH-methods to find netcdf-directory +EXEC_PROGRAM(${CMAKE_CURRENT_SOURCE_DIR}/findnetcdf OUTPUT_VARIABLE ADDMODULEPATH) + +### Find NetCDF files +FIND_PATH(NETCDF_INCLUDE_DIR netcdf.mod NETCDF.mod + PATHS + $ENV{SARA_NETCDF_INCLUDE} + $ENV{SURFSARA_NETCDF_INCLUDE} + $ENV{NETCDF_INCLUDE} + ${ADDMODULEPATH}/include + /usr/include + $ENV{HOME}/include + /usr/lib64/gfortran/modules + /opt/rrzk/lib/netcdf/4.1.3-gcc-4.8.2/f90/ + DOC "NetCDF include directory (must contain netcdf.mod)" +) + +FIND_LIBRARY(NETCDF_C_LIB netcdf + PATHS + $ENV{SARA_NETCDF_LIB} + $ENV{SURFSARA_NETCDF_LIB} + $ENV{NETCDF_LIB} + ${ADDMODULEPATH}/lib + ${ADDMODULEPATH}/lib64 + /usr/lib + /usr/lib64 + $ENV{HOME}/lib + $ENV{HOME}/lib64 + /opt/rrzk/lib/netcdf/4.1.3-gcc-4.8.2/lib/ + DOC "NetCDF C library" +) + +FIND_LIBRARY(NETCDF_FORTRAN_LIB netcdff + PATHS + $ENV{SARA_NETCDF_LIB} + $ENV{SURFSARA_NETCDF_LIB} + $ENV{NETCDF_LIB} + ${ADDMODULEPATH}/lib + ${ADDMODULEPATH}/lib64 + /usr/lib + /usr/lib64 + $ENV{HOME}/lib + $ENV{HOME}/lib64 + /opt/rrzk/lib/netcdf/4.1.3-gcc-4.8.2/lib/ + DOC "NetCDF Fortran library" +) + +if(NETCDF_INCLUDE_DIR) + include_directories(${NETCDF_INCLUDE_DIR}) +else(NETCDF_INCLUDE_DIR) + MESSAGE(STATUS "WARNING: No NETCDF bindings are found.") +endif(NETCDF_INCLUDE_DIR) + +if(NETCDF_C_LIB) + set(NETCDF_LIBS ${NETCDF_C_LIB}) +else(NETCDF_C_LIB) + MESSAGE(STATUS "WARNING: No NETCDF bindings are found.") +endif(NETCDF_C_LIB) + +if(NETCDF_FORTRAN_LIB) + set(NETCDF_LIBS ${NETCDF_LIBS} ${NETCDF_FORTRAN_LIB}) +else(NETCDF_FORTRAN_LIB) + MESSAGE(STATUS "WARNING: No Fortran NETCDF bindings are found.") +endif(NETCDF_FORTRAN_LIB) + +### Documentation +INCLUDE(FindDoxygen) +if(DOXYGEN) + ADD_SUBDIRECTORY(utils/doc) +else() + MESSAGE(STATUS "WARNING: Doxygen not found - Reference manual will not be created") +endif() + +### Set case +if(NOT CASE) + set (CASE standard CACHE STRING + "Set the case." + FORCE) +endif() + +### Add case specific file +FILE(GLOB usrfile "${CMAKE_CURRENT_SOURCE_DIR}/cases/${CASE}/moduser.f90") +if(usrfile STREQUAL "") + set(usrfile "${CMAKE_CURRENT_SOURCE_DIR}/cases/standard/moduser.f90") +endif() +execute_process(COMMAND ${CMAKE_COMMAND} -E copy_if_different ${usrfile} ${CMAKE_CURRENT_SOURCE_DIR}/src/moduser.f90) +MESSAGE(STATUS "Case " ${CASE} " uses " ${usrfile}) + +ADD_SUBDIRECTORY(src) diff --git a/config/modversion.f90 b/config/modversion.f90 new file mode 100644 index 00000000..47239b35 --- /dev/null +++ b/config/modversion.f90 @@ -0,0 +1,12 @@ +! This file is used on Cheops and Juwels +! to supplement missing Git/Cmake functionalities +! Variables in this file are filled in by CMake at build time, +! used to record the git version and hash. + +module modversion + + implicit none + + character(80) :: git_version='2.17' ! "@GIT_VERSION@" + +end module modversion diff --git a/rundales/runDALES_parallel_16_cheops_master.sh b/rundales/runDALES_parallel_16_cheops_master.sh new file mode 100755 index 00000000..501ccde5 --- /dev/null +++ b/rundales/runDALES_parallel_16_cheops_master.sh @@ -0,0 +1,27 @@ +#!/bin/bash -l +#SBATCH --ntasks=16 +#SBATCH --mem-per-cpu=20gb +#SBATCH --time=120:00:00 +#SBATCH --account=UniKoeln + +# number of nodes in $SLURM_NNODES (default: 1) +# number of tasks in $SLURM_NTASKS (default: 1) +# number of tasks per node in $SLURM_NTASKS_PER_NODE (default: 1) +# number of threads per task in $SLURM_CPUS_PER_TASK (default: 1) + + +#module avail +#module load intelmpi/4.1.0 netcdf/4.1.3 + module load gnu/4.8.2 netcdf/4.1.3-gcc-4.8.2 intelmpi/2018 + module load hdf5/1.8.11 szlib +# module load cmake + +#echo $PATH + +echo "# of nodes, tasks:" $SLURM_NNODES $SLURM_NTASKS + +#LD_LIBRARY_PATH="$LD_LIBRARY_PATH:/home/rneggers/bin/netcdf-4.3.0_ifort/lib" +#LD_LIBRARY_PATH="$LD_LIBRARY_PATH:/opt/rrzk/lib/netcdf/4.1.3/lib" +LD_LIBRARY_PATH="$LD_LIBRARY_PATH:/opt/rrzk/lib/netcdf/4.1.3-gcc-4.8.2/lib" + +srun -n $SLURM_NTASKS ./dales4 diff --git a/src/modbulkmicro3.f90 b/src/modbulkmicro3.f90 index 94fb221b..b75b9794 100644 --- a/src/modbulkmicro3.f90 +++ b/src/modbulkmicro3.f90 @@ -77,7 +77,7 @@ subroutine initbulkmicro3 ,l_sb_all_or, l_sb_dbg & ,l_setclouds, l_setccn & ! flag whether to set cloud ,l_corr_neg_qt & ! flag whether to adjust qt and thlp in hydrometeor corrections - ,l_sb_lim_aggr, l_sb_stickyice & ! ice and snow aggregation flags + ,l_snow, l_sb_lim_aggr, l_sb_stickyice & ! ice and snow aggregation flags ,l_sb_conv_par & ! conversion flags ,l_c_ccn,l_sb_sat_max & ! flags for cloud nucleation ,l_sb_nuc_sat,l_sb_nuc_expl,l_sb_nuc_diff & ! flags for cloud nucleation @@ -156,6 +156,7 @@ subroutine initbulkmicro3 call MPI_BCAST(l_ccn_init, 1, MPI_LOGICAL ,0,comm3d,ierr) call MPI_BCAST(l_corr_neg_qt, 1, MPI_LOGICAL ,0,comm3d,ierr) ! + call MPI_BCAST(l_snow, 1, MPI_LOGICAL ,0,comm3d,ierr) call MPI_BCAST(l_sb_lim_aggr, 1, MPI_LOGICAL ,0,comm3d,ierr) call MPI_BCAST(l_sb_stickyice, 1, MPI_LOGICAL ,0,comm3d,ierr) call MPI_BCAST(l_sb_conv_par , 1, MPI_LOGICAL ,0,comm3d,ierr) @@ -169,7 +170,7 @@ subroutine initbulkmicro3 call MPI_BCAST(l_sb_inuc_expl , 1, MPI_LOGICAL ,0,comm3d,ierr) call MPI_BCAST(l_sb_reisner , 1, MPI_LOGICAL ,0,comm3d,ierr) ! - call MPI_BCAST(l_sb_reisner, 1, MPI_LOGICAL ,0,comm3d,ierr) + ! call MPI_BCAST(l_sb_reisner, 1, MPI_LOGICAL ,0,comm3d,ierr) ! call MPI_BCAST(N_inuc_R , 1, MY_REAL ,0,comm3d,ierr) call MPI_BCAST(c_inuc_R , 1, MY_REAL ,0,comm3d,ierr) @@ -930,7 +931,7 @@ subroutine bulkmicro3 !********************************************************************* ! remove negative values for hydrometereors and clouds !********************************************************************* - if (l_rain) then + if (l_rain) then ! old: ! if (sum(qr, qr<0.) > 0.000001*sum(qr)) then ! write(*,*)'amount of neg. qr and Nr thrown away is too high ',timee,' sec' @@ -1138,48 +1139,31 @@ subroutine bulkmicro3 ! if (any((qt0(2:i1,2:j1,1:k1)-qvsl(2:i1,2:j1,1:k1)-svm(2:i1,2:j1,1:k1,iq_cl)-svm(2:i1,2:j1,1:k1,iq_ci)).lt. 0.)) then ! write(6,*) 'undersaturated clouds in gridpoints ',count((qt0(2:i1,2:j1,1:k1)-qvsl(2:i1,2:j1,1:k1)-svm(2:i1,2:j1,1:k1,iq_cl)-svm(2:i1,2:j1,1:k1,iq_ci)).lt. 0.) ! endif - if (any((qt0(2:i1,2:j1,1:k1)-qvsl(2:i1,2:j1,1:k1)-q_cl(2:i1,2:j1,1:k1)-q_ci(2:i1,2:j1,1:k1)).lt. 0.)) then - write(6,*) 'undersaturated clouds in gridpoints ',count((qt0(2:i1,2:j1,1:k1)-qvsl(2:i1,2:j1,1:k1)-q_cl(2:i1,2:j1,1:k1)-q_ci(2:i1,2:j1,1:k1)).lt. 0.) - endif + if (any((qt0(2:i1,2:j1,1:k1)-qvsl(2:i1,2:j1,1:k1)-q_cl(2:i1,2:j1,1:k1)-q_ci(2:i1,2:j1,1:k1)).lt. 0.)) then + write(6,*) 'undersaturated clouds in gridpoints ',count((qt0(2:i1,2:j1,1:k1)-qvsl(2:i1,2:j1,1:k1)-q_cl(2:i1,2:j1,1:k1)-q_ci(2:i1,2:j1,1:k1)).lt. 0.) + endif - do k=1,k1 - do j=2,j1 - do i=2,i1 - nrtest=n_cl(i,j,k) ! considering n_clp only - ! in clouds only - if (nrtest.le.0.0) then + do k=1,k1 + do j=2,j1 + do i=2,i1 + nrtest=n_cl(i,j,k) ! considering n_clp only + ! in clouds only + if (nrtest.le.0.0) then xtest(i,j,k) =0.0 - else + else ! calculating new droplet size xtest(i,j,k) = q_cl(i,j,k) / ( n_cl(i,j,k)+eps0) ! o: rhof(k) * (svm(i,j,k,iqsv) + delt*q_clp(i,j,k)) / (svm(i,j,k,insv)+delt*n_clp(i,j,k)) - endif - enddo - enddo - enddo - ! warning if droplets too large - if(any( xtest(2:i1,2:j1,1:k1) .gt. xc_bmax)) then + endif + enddo + enddo + enddo + ! warning if droplets too large + if(any( xtest(2:i1,2:j1,1:k1) .gt. xc_bmax)) then write(6,*) 'warning: cloud droplets at the start of the step too large' write(6,*) ' in gridpoints ', count(xtest(2:i1,2:j1,1:k1).gt.xc_bmax), ' larger than ',xc_bmax write(6,*) ' max value ', maxval(xtest), ' in ', maxloc(xtest) - ! #dbg START and list it -! do k=1,k1 -! do j=2,j1 -! do i=2,i1 -! if( xc(i,j,k) .gt. xc_bmax) then -! write(6,*) 'i,j,k =', i, j, k, ' xc1= ', xc(i,j,k), ' ql0= ', ql0(i,j,k) -! write(6,*) ' n_cl-1 = ', svm(i,j,k,insv), ' n_cl0 = ', sv0(i,j,k,insv), ' n_cl1 = ', svm(i,j,k,insv)+delt*n_clp(i,j,k) ! svp(i,j,k,insv) -! write(6,*) ' q_cl-1 = ', svm(i,j,k,iqsv), ' q_cl0 = ', sv0(i,j,k,iqsv), ' q_cl1 = ', svm(i,j,k,iqsv)+delt*q_clp(i,j,k) ! svp(i,j,k,iqsv) -! write(6,*) ' n_ci-1 = ', svm(i,j,k,in_ci), ' n_ci0 = ', sv0(i,j,k,in_ci), ' n_ci1 = ', svm(i,j,k,in_ci)+delt*n_cip(i,j,k) ! svp(i,j,k,insv) -! write(6,*) ' q_ci-1 = ', svm(i,j,k,iq_ci), ' q_ci0 = ', sv0(i,j,k,iq_ci), ' q_ci1 = ', svm(i,j,k,iq_ci)+delt*q_cip(i,j,k) ! svp(i,j,k,iqsv) -! write(6,*) ' delta q_cl1_mphys= ', delt*q_clp(i,j,k),' delta q_cl1_other= ', delt*svp(i,j,k,iqsv) -! ! write(6,*) ' sat adjustment ', (1.0/delt)*(ql0(i,j,k)-q_cl(i,j,k)-q_ci(i,j,k))-max(0.0,(q_clp(i,j,k)+q_cip(i,j,k))-qtpmcr(i,j,k)) -! ! i.e. listing water already in clouds vs total amount of water -! endif ! #dbg END -! enddo -! enddo -! enddo ! #dbg END - endif + endif endif ! l_sb_dbg_extra @@ -1232,83 +1216,59 @@ subroutine bulkmicro3 ! -------------------------------------------------------------- ! cloud forming processes ! -------------------------------------------------------------- - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'before.')! #d call nucleation3 ! cloud nucleation - if(l_sb_dbg_extra)call check_nan(flag_do_dbg,flag_nan,'nuclea.')! #d - call icenucle3 ! ice nucleation ! #Bb ! #iceout - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'ice nuc')! #d + call icenucle3 ! ice nucleation ! #Bb ! #iceout ! not needed -treated automatically in nucleation_3 !call cor_nucl3 ! limits nucleation so there are still ccn left - ! if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'cor.nuc')! #d ! -------------------------------------------------------------- ! freezing of water droplets ! -------------------------------------------------------------- - call homfreez3 ! homogeneous freezing of cloud droplets ! #iceout - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'hom.fre')! #d + call homfreez3 ! homogeneous freezing of cloud droplets ! #iceout call hetfreez3 ! heterogeneous freezing! #iceout - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'het.fre')! #d - ! -------------------------------------------------------------- ! deposition processes ! -------------------------------------------------------------- call deposit_ice3 ! deposition of vapour to cloud ice - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'ice dep')! #d call deposit_snow3 ! deposition of vapour to snow - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'sno.dep')! #d call deposit_graupel3 ! deposition of vapour to graupel - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'gra.dep')! #d call cor_deposit3 ! correction for deposition - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'cor.con')! #d ! -------------------------------------------------------------- ! snow aggregation and self-collection ! -------------------------------------------------------------- - call ice_aggr3 ! ice selfcollection !#b5 ! #a4d1 ! #Ba !#b2t2 !#b2t1 - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'ice agg')! ice aggregation - tendencies only + if(l_snow) call ice_aggr3 ! ice selfcollection !#b5 ! #a4d1 ! #Ba !#b2t2 !#b2t1 call snow_self3 ! snow selfcollection - tendency only - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'snow sc')! #d ! -------------------------------------------------------------- ! collision processes for snow ! -------------------------------------------------------------- - call coll_sis3 ! snow selfcollection !#t1 - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'s+i->s ')! #d - - call coll_gsg3 ! snow selfcollection !#t1 - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'g+s->g ')! #d + call coll_sis3 ! snow selfcollection !#t1 + call coll_gsg3 ! snow selfcollection !#t1 ! -------------------------------------------------------------- ! - rimings ! -------------------------------------------------------------- - call coll_ici3 ! riming i+c -> i ! #t1 - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'i+c->i ')! #d + call coll_ici3 ! riming i+c -> i ! #t1 call coll_scs3 ! riming s+c -> s ! #t1 - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'s+c->s ')! #d call coll_gcg3 ! riming g+c -> g !#t2 !#t4 - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'g+c->g ')! #d !o call rime_srs3 ! riming g+r -> g !#t2 !#t4 !o if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'s+r->s ')! #d call coll_grg3 ! riming g+r -> g !#t2 !#t4 - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'g+r->g ')! #d ! -------------------------------------------------------------- ! - raindrop freezing ! -------------------------------------------------------------- - call rainhetfreez3 ! heterogeneous freezing! #iceout - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'fre.rai')! #d - + call rainhetfreez3 ! heterogeneous freezing! #iceout ! -------------------------------------------------------------- ! - collision with conversion ! -------------------------------------------------------------- call coll_rig3 ! riming r+i -> g !#t2 !#t4 - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'r+i->g ')! #d - call coll_rsg3 ! riming r+i -> g !#t7 - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'r+s->g ')! #d + call coll_rsg3 ! riming r+i -> g !#t7 !-------------------------------------------------------------- ! conversions @@ -1316,47 +1276,35 @@ subroutine bulkmicro3 ! ice multiplication of Hallet and Mossop call ice_multi3 - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'ice.mul')! #d ! partial conversion - if (l_sb_conv_par) then - call conv_partial3 !#t3 !#b2t2 - endif - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'par.con')! #d + if (l_sb_conv_par) call conv_partial3 !#t3 !#b2t2 ! -------------------------------------------------------------- ! - melting and evaporation of ice particles ! -------------------------------------------------------------- ! melting of ice particles - call evapmelting3 - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'melting')! #d - + call evapmelting3 ! -------------------------------------------------------------- ! - basic warm processes ! -------------------------------------------------------------- call autoconversion3 ! call bulkmicrotend - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'autocon')! #d call cloud_self3 ! call bulkmicrotend - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'cl self')! #d call accretion3 ! call bulkmicrotend - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'accret.')! #d - ! rain evaporation call evap_rain3 - ! call bulkmicrotend - if(l_sb_dbg_extra) call check_nan(flag_do_dbg,flag_nan,'evap_r.')! #d + ! call bulkmicrotend ! ============================= ! saturation adjustment ! ============================= call satadj3 - !================================== ! sedimentation part !================================== @@ -1374,25 +1322,7 @@ subroutine bulkmicro3 !******************************************* ! statistics microphysics !******************************************* - - - - ! - !#d2 call cor_loss_cl3 - ! #d4 call check_update(flag_do_dbg,flag_warn_update,'cor.lcl')! #d - ! #d1 call check_sizes(flag_do_dbg,flag_warn_size,'cor.lcl') ! #d - ! call check_nan(flag_do_dbg,flag_nan,'cor.lcl') ! #d - ! : call bulkmicrotend - ! #t_evap: testing without evaporation of rain - - ! : call bulkmicrotend - ! - ! - and finally the saturation adjustment - ! - actually within this code - ! #sb3 END - ! endif ! <- it should happen even with rain off - - ! saving water number and mixing ratios into + ! : call bulkmicrotend ! #sb3 ! ! svp(2:i1,2:j1,1:k1,in_cl)=n_clp(2:i1,2:j1,1:k1) @@ -1450,9 +1380,7 @@ subroutine bulkmicro3 end if enddo enddo - enddo - ! check update - ! #d4 call check_test(flag_do_dbg,flag_warn_test,'test_r.') ! #d + enddo ! #sb3 START ! testing negative values also for other variables ! @@ -1481,8 +1409,6 @@ subroutine bulkmicro3 enddo enddo enddo - ! - ! #d4 call check_test(flag_do_dbg,flag_warn_test,'test_s.') ! #d ! ! == graupel == iqsv = iq_hg @@ -1510,8 +1436,6 @@ subroutine bulkmicro3 enddo enddo ! - ! #d4 call check_test(flag_do_dbg,flag_warn_test,'test_g.')! #d - ! ! == cloud ice == iqsv = iq_ci insv = in_ci @@ -1544,7 +1468,6 @@ subroutine bulkmicro3 enddo enddo ! - ! #d4 call check_test(flag_do_dbg,flag_warn_test,'test_ci')! #d ! ! == cloud liquid water == iqsv = iq_cl @@ -1615,8 +1538,6 @@ subroutine bulkmicro3 enddo enddo enddo - ! check update - ! #d4 call check_test(flag_do_dbg,flag_warn_test,'test_r.') ! #d ! #sb3 START ! testing negative values also for other variables ! @@ -1645,8 +1566,6 @@ subroutine bulkmicro3 enddo enddo enddo - ! - ! #d4 call check_test(flag_do_dbg,flag_warn_test,'test_s.') ! #d ! ! == graupel == iqsv = iq_hg @@ -1674,8 +1593,6 @@ subroutine bulkmicro3 enddo enddo ! - ! #d4 call check_test(flag_do_dbg,flag_warn_test,'test_g.')! #d - ! ! == cloud ice == iqsv = iq_ci insv = in_ci @@ -1708,8 +1625,6 @@ subroutine bulkmicro3 enddo enddo ! - ! #d4 call check_test(flag_do_dbg,flag_warn_test,'test_ci')! #d - ! ! == cloud liquid water == iqsv = iq_cl insv = in_cl @@ -1737,8 +1652,6 @@ subroutine bulkmicro3 enddo endif ! l_corr__neg_qt ! - ! #d4 call check_test(flag_do_dbg,flag_warn_test,'test_cl')! #d - ! ! == CCN checking == if(.not.l_c_ccn) then @@ -1759,7 +1672,6 @@ subroutine bulkmicro3 endif ! not l_c_ccn ! #sb3 END - ! call check_ssat(flag_do_dbg) ! #sb3 ! == update ===== ! and updating main prognostic variables by contribution from mphys processes @@ -1771,55 +1683,7 @@ subroutine bulkmicro3 enddo enddo enddo - - ! #t START - if(l_sb_dbg_extra) then - ! estimating new xc - insv = in_cl - iqsv = iq_cl - do k=1,k1 - do j=2,j1 - do i=2,i1 - nrtest=svm(i,j,k,insv)+delt*svp(i,j,k,insv) ! since n_clp already included - ! in clouds only - if (nrtest<0.0 ) then - xtest(i,j,k) =0.0 - else - ! calculating new droplet size - xtest(i,j,k) = (svm(i,j,k,iqsv) + delt*svp(i,j,k,iqsv)) / (svm(i,j,k,insv)+delt*svp(i,j,k,insv)+eps0) - ! rhof(k) * (svm(i,j,k,iqsv) + delt*svp(i,j,k,iqsv)) / (svm(i,j,k,insv)+delt*svp(i,j,k,insv)) - endif - enddo - enddo - enddo - ! warning if droplets too large - if(any( xtest(2:i1,2:j1,1:k1) .gt. xc_bmax)) then - write(6,*) 'warning: cloud droplets at the end of the step too large' - write(6,*) ' in ', count(xtest(2:i1,2:j1,1:k1).gt.xc_bmax), ' gridpoints larger than ',xc_bmax - write(6,*) ' max value ', maxval(xtest), ' in ', maxloc(xtest) - ! #dbg START and list it -! do k=1,k1 -! do j=2,j1 -! do i=2,i1 -! if( xc(i,j,k) .gt. xc_bmax) then -! write(6,*) 'i,j,k =', i, j, k, ' xc1= ', xc(i,j,k), ' ql0= ', ql0(i,j,k) -! write(6,*) ' n_cl-1 = ', svm(i,j,k,insv), ' n_cl0 = ', sv0(i,j,k,insv), ' n_cl1 = ', svm(i,j,k,insv)+delt*svp(i,j,k,insv) -! write(6,*) ' q_cl-1 = ', svm(i,j,k,iqsv), ' q_cl0 = ', sv0(i,j,k,iqsv), ' q_cl1 = ', svm(i,j,k,iqsv)+delt*svp(i,j,k,iqsv) -! write(6,*) ' n_ci-1 = ', svm(i,j,k,in_ci), ' n_ci0 = ', sv0(i,j,k,in_ci), ' n_ci1 = ', svm(i,j,k,in_ci)+delt*svp(i,j,k,in_ci) -! write(6,*) ' q_ci-1 = ', svm(i,j,k,iq_ci), ' q_ci0 = ', sv0(i,j,k,iq_ci), ' q_ci1 = ', svm(i,j,k,iq_ci)+delt*svp(i,j,k,iq_ci) -! write(6,*) ' delta q_cl1_mphys= ', delt*q_clp(i,j,k),' delta q_cl1_all= ', delt*svp(i,j,k,iqsv) -! write(6,*) ' qt0 = ', qt0(i,j,k),' delta q_t = ' , delt*qtp(i,j,k) -! write(6,*) ' by qtpmcr = ', delt*qtp(i,j,k), ' by qtpmcr = ', delt*qtpmcr(i,j,k) -! ! write(6,*) ' old sat adjustment = ', (ql0(i,j,k)-q_cl(i,j,k)-q_ci(i,j,k))-delt*max(0.0,(q_clp(i,j,k)+q_cip(i,j,k))-qtpmcr(i,j,k)) -! ! write(6,*) ' new sat adjustment = ', -! ! i.e. listing water already in clouds vs total amount of water -! endif -! enddo -! enddo -! enddo ! #dbg END - endif - endif ! l_sb_dbg - ! #t END + ! # statistics @@ -1828,12 +1692,12 @@ subroutine bulkmicro3 call bulkmicrostat3 ! #t5 - if(l_sb_dbg_extra) call check_sizes(flag_do_dbg,flag_warn_size,'all_ups') - call check_allupdates(flag_do_dbg_up) ! if(l_sb_dbg) call check_allupdates(flag_do_dbg) + if(l_sb_dbg) call check_sizes(flag_do_dbg,flag_warn_size,'all_ups') + if(l_sb_dbg) call check_allupdates(flag_do_dbg_up) ! if(l_sb_dbg) call check_allupdates(flag_do_dbg) ! ending line - if(l_sb_dbg_extra) then - write(6,*) ' ==========================' - endif + !if(l_sb_dbg) then + ! write(6,*) ' ==========================' + !endif ! deallocating process variables deallocate( xtest) ! #d @@ -9648,190 +9512,7 @@ subroutine check_allupdates(flag) endif ! flag endif - end subroutine check_allupdates - - - ! checking for strange of incorrect resuls - subroutine check_test(flag_dbg,flag, proc_name) - use modglobal, only : ih,i1,jh,j1,k1 - use modfields, only : sv0,svm,svp,qtp,thlp,ql0,exnf,rhof,qvsl,qvsi,qt0, thl0 - implicit none - logical, intent (inout) ::flag_dbg,flag - character (len=7), intent (in) :: proc_name - real :: lim_thlp, lim_qtp, ssat, ssat_ice, uncond - integer:: i,j,k - - ! inner setting - lim_thlp = 4.0 - lim_qtp = 0.0005 - - if(flag_dbg) then - - ! checking for large updates - if(any(thlp(2:i1,2:j1,1:k1).gt.lim_thlp).or.any(thlp(2:i1,2:j1,1:k1).lt.(-lim_thlp))) then - write(6,*) 'WARNING: thl problem in: ', proc_name - if (flag) then - write(6,*) 'problems in: ',count(thlp(2:i1,2:j1,1:k1).lt.(-lim_qtp)),'and',count(thlp(2:i1,2:j1,1:k1).gt.lim_qtp) - else ! flag - flag = .true. - do k=1,k1 - do j=2,j1 - do i=2,i1 - if (thlp(i,j,k).gt.lim_thlp) then - ! calculate - uncond=ql0(i,j,k)-q_cl(i,j,k) - ssat = 100.0*((qt0(i,j,k)-q_cl(i,j,k))/qvsl(i,j,k)-1.0) - ssat_ice = 100.0*((qt0(i,j,k)-q_cl(i,j,k))/qvsi(i,j,k)-1.0) - !#iceout uncond=ql0(i,j,k)-q_cl(i,j,k)-q_ci(i,j,k) - !#iceout ssat = 100.0*((qt0(i,j,k)-q_cl(i,j,k)-q_ci(i,j,k))/qvsl(i,j,k)-1.0) - !#iceout ssat_ice = 100.0*((qt0(i,j,k)-q_cl(i,j,k)-q_ci(i,j,k))/qvsi(i,j,k)-1.0) - ! outputs - write(6,*) ' (i,j,k) = ',i,j,k,' thlp = ', thlp(i,j,k) - write(6,*) ' uncondensed = ', uncond,' ssat = ',ssat,'% ',' ssat_i = ',ssat_ice,'% ' - write(6,*) ' qt = ',qt0(i,j,k),' thl = ', thl0(i,j,k), ' n_ccn = ', n_cc(i,j,k) - write(6,*) ' x_cl =', x_cl(i,j,k), ' n_cl =', n_cl(i,j,k),' q_cl =', q_cl (i,j,k) - write(6,*) ' x_ci =', x_ci(i,j,k), ' n_ci =', n_ci(i,j,k),' q_ci =', q_ci (i,j,k) - write(6,*) ' x_hr =', x_hr(i,j,k), ' n_hr =', n_hr(i,j,k),' q_hr =', q_hr (i,j,k) - write(6,*) ' x_hs =', x_hs(i,j,k), ' n_hs =', n_hs(i,j,k),' q_hs =', q_hs (i,j,k) - write(6,*) ' x_hg =', x_hg(i,j,k), ' n_hg =', n_hg(i,j,k),' q_hg =', q_hg (i,j,k) - write(6,*) ' q_clp =',q_clp(i,j,k),' n_clp =', n_clp(i,j,k),' n_cl_m =', svm(i,j,k,in_cl),' q_cl_m =',svm(i,j,k,iq_cl) - write(6,*) ' q_cip =',q_cip(i,j,k),' n_cip =', n_cip(i,j,k),' n_ci_m =', svm(i,j,k,in_ci),' q_ci_m =',svm(i,j,k,iq_ci) - write(6,*) ' q_hrp =',q_hrp(i,j,k),' n_hrp =', n_hrp(i,j,k),' n_hr_m =', svm(i,j,k,in_hr),' q_hr_m =',svm(i,j,k,iq_hr) - write(6,*) ' q_hsp =',q_hsp(i,j,k),' n_hsp =', n_hsp(i,j,k),' n_hs_m =', svm(i,j,k,in_hs),' q_hs_m =',svm(i,j,k,iq_hs) - write(6,*) ' q_hgp =',q_hgp(i,j,k),' n_hgp =', n_hgp(i,j,k),' n_hg_m =', svm(i,j,k,in_hg),' q_hg_m =',svm(i,j,k,iq_hg) - endif - if (thlp(i,j,k).lt.(-lim_thlp)) then - ! calculate - uncond=ql0(i,j,k)-q_cl(i,j,k) - ssat = 100.0*((qt0(i,j,k)-q_cl(i,j,k))/qvsl(i,j,k)-1.0) - ssat_ice = 100.0*((qt0(i,j,k)-q_cl(i,j,k))/qvsi(i,j,k)-1.0) - !#iceout uncond=ql0(i,j,k)-q_cl(i,j,k)-q_ci(i,j,k) - !#iceout ssat = 100.0*((qt0(i,j,k)-q_cl(i,j,k)-q_ci(i,j,k))/qvsl(i,j,k)-1.0) - !#iceout ssat_ice = 100.0*((qt0(i,j,k)-q_cl(i,j,k)-q_ci(i,j,k))/qvsi(i,j,k)-1.0) - ! outputs - write(6,*) ' (i,j,k) = ',i,j,k,' thlp = ', thlp(i,j,k) - write(6,*) ' uncondensed = ', uncond,' ssat = ',ssat,'% ',' ssat_i = ',ssat_ice,'% ' - write(6,*) ' qt = ',qt0(i,j,k),' thl = ', thl0(i,j,k), ' n_ccn = ', n_cc(i,j,k) - write(6,*) ' x_cl =', x_cl(i,j,k), ' n_cl =', n_cl(i,j,k),' q_cl =', q_cl (i,j,k) - write(6,*) ' x_ci =', x_ci(i,j,k), ' n_ci =', n_ci(i,j,k),' q_ci =', q_ci (i,j,k) - write(6,*) ' x_hr =', x_hr(i,j,k), ' n_hr =', n_hr(i,j,k),' q_hr =', q_hr (i,j,k) - write(6,*) ' x_hs =', x_hs(i,j,k), ' n_hs =', n_hs(i,j,k),' q_hs =', q_hs (i,j,k) - write(6,*) ' x_hg =', x_hg(i,j,k), ' n_hg =', n_hg(i,j,k),' q_hg =', q_hg (i,j,k) - write(6,*) ' q_clp =',q_clp(i,j,k),' n_clp =', n_clp(i,j,k),' n_cl_m =', svm(i,j,k,in_cl),' q_cl_m =',svm(i,j,k,iq_cl) - write(6,*) ' q_cip =',q_cip(i,j,k),' n_cip =', n_cip(i,j,k),' n_ci_m =', svm(i,j,k,in_ci),' q_ci_m =',svm(i,j,k,iq_ci) - write(6,*) ' q_hrp =',q_hrp(i,j,k),' n_hrp =', n_hrp(i,j,k),' n_hr_m =', svm(i,j,k,in_hr),' q_hr_m =',svm(i,j,k,iq_hr) - write(6,*) ' q_hsp =',q_hsp(i,j,k),' n_hsp =', n_hsp(i,j,k),' n_hs_m =', svm(i,j,k,in_hs),' q_hs_m =',svm(i,j,k,iq_hs) - write(6,*) ' q_hgp =',q_hgp(i,j,k),' n_hgp =', n_hgp(i,j,k),' n_hg_m =', svm(i,j,k,in_hg),' q_hg_m =',svm(i,j,k,iq_hg) - endif - enddo - enddo - enddo - endif ! flag - endif - - if(any(qtp(2:i1,2:j1,1:k1).lt.(-lim_qtp)).or.any(qtp(2:i1,2:j1,1:k1).gt.lim_qtp)) then - write(6,*) 'WARNING: qtp problem in: ', proc_name - if (flag) then - write(6,*) 'problems in: ',count(qtp(2:i1,2:j1,1:k1).lt.(-lim_qtp)),'and',count(qtp(2:i1,2:j1,1:k1).gt.lim_qtp) - else ! flag - flag = .true. - do k=1,k1 - do j=2,j1 - do i=2,i1 - if (qtp(i,j,k).gt.lim_qtp) then - ! calculate - uncond=ql0(i,j,k)-q_cl(i,j,k) - ssat = 100.0*((qt0(i,j,k)-q_cl(i,j,k))/qvsl(i,j,k)-1.0) - ssat_ice = 100.0*((qt0(i,j,k)-q_cl(i,j,k))/qvsi(i,j,k)-1.0) - !#iceout uncond=ql0(i,j,k)-q_cl(i,j,k)-q_ci(i,j,k) - !#iceout ssat = 100.0*((qt0(i,j,k)-q_cl(i,j,k)-q_ci(i,j,k))/qvsl(i,j,k)-1.0) - !#iceout ssat_ice = 100.0*((qt0(i,j,k)-q_cl(i,j,k)-q_ci(i,j,k))/qvsi(i,j,k)-1.0) - ! outputs - write(6,*) ' (i,j,k) = ',i,j,k,' qtp = ', qtp (i,j,k) - write(6,*) ' uncondensed = ', uncond,' ssat = ',ssat,'% ',' ssat_i = ',ssat_ice,'% ' - write(6,*) ' qt = ',qt0(i,j,k),' thl = ', thl0(i,j,k), ' n_ccn = ', n_cc(i,j,k) - write(6,*) ' x_cl =', x_cl(i,j,k), ' n_cl =', n_cl(i,j,k),' q_cl =', q_cl (i,j,k) - write(6,*) ' x_ci =', x_ci(i,j,k), ' n_ci =', n_ci(i,j,k),' q_ci =', q_ci (i,j,k) - write(6,*) ' x_hr =', x_hr(i,j,k), ' n_hr =', n_hr(i,j,k),' q_hr =', q_hr (i,j,k) - write(6,*) ' x_hs =', x_hs(i,j,k), ' n_hs =', n_hs(i,j,k),' q_hs =', q_hs (i,j,k) - write(6,*) ' x_hg =', x_hg(i,j,k), ' n_hg =', n_hg(i,j,k),' q_hg =', q_hg (i,j,k) - write(6,*) ' q_clp =',q_clp(i,j,k),' n_clp =', n_clp(i,j,k),' n_cl_m =', svm(i,j,k,in_cl),' q_cl_m =',svm(i,j,k,iq_cl) - write(6,*) ' q_cip =',q_cip(i,j,k),' n_cip =', n_cip(i,j,k),' n_ci_m =', svm(i,j,k,in_ci),' q_ci_m =',svm(i,j,k,iq_ci) - write(6,*) ' q_hrp =',q_hrp(i,j,k),' n_hrp =', n_hrp(i,j,k),' n_hr_m =', svm(i,j,k,in_hr),' q_hr_m =',svm(i,j,k,iq_hr) - write(6,*) ' q_hsp =',q_hsp(i,j,k),' n_hsp =', n_hsp(i,j,k),' n_hs_m =', svm(i,j,k,in_hs),' q_hs_m =',svm(i,j,k,iq_hs) - write(6,*) ' q_hgp =',q_hgp(i,j,k),' n_hgp =', n_hgp(i,j,k),' n_hg_m =', svm(i,j,k,in_hg),' q_hg_m =',svm(i,j,k,iq_hg) - ! - endif - if (qtp(i,j,k).lt.(-lim_qtp)) then - ! calculate - uncond=ql0(i,j,k)-q_cl(i,j,k) - ssat = 100.0*((qt0(i,j,k)-q_cl(i,j,k))/qvsl(i,j,k)-1.0) - ssat_ice = 100.0*((qt0(i,j,k)-q_cl(i,j,k))/qvsi(i,j,k)-1.0) - !#iceout uncond=ql0(i,j,k)-q_cl(i,j,k)-q_ci(i,j,k) - !#iceout ssat = 100.0*((qt0(i,j,k)-q_cl(i,j,k)-q_ci(i,j,k))/qvsl(i,j,k)-1.0) - !#iceout ssat_ice = 100.0*((qt0(i,j,k)-q_cl(i,j,k)-q_ci(i,j,k))/qvsi(i,j,k)-1.0) - ! outputs - write(6,*) ' (i,j,k) = ',i,j,k,' qtp = ', qtp (i,j,k) - write(6,*) ' uncondensed = ', uncond,' ssat = ',ssat,'% ',' ssat_i = ',ssat_ice,'% ' - write(6,*) ' qt = ',qt0(i,j,k),' thl = ', thl0(i,j,k), ' n_ccn = ', n_cc(i,j,k) - write(6,*) ' x_cl =', x_cl(i,j,k), ' n_cl =', n_cl(i,j,k),' q_cl =', q_cl (i,j,k) - write(6,*) ' x_ci =', x_ci(i,j,k), ' n_ci =', n_ci(i,j,k),' q_ci =', q_ci (i,j,k) - write(6,*) ' x_hr =', x_hr(i,j,k), ' n_hr =', n_hr(i,j,k),' q_hr =', q_hr (i,j,k) - write(6,*) ' x_hs =', x_hs(i,j,k), ' n_hs =', n_hs(i,j,k),' q_hs =', q_hs (i,j,k) - write(6,*) ' x_hg =', x_hg(i,j,k), ' n_hg =', n_hg(i,j,k),' q_hg =', q_hg (i,j,k) - write(6,*) ' q_clp =',q_clp(i,j,k),' n_clp =', n_clp(i,j,k),' n_cl_m =', svm(i,j,k,in_cl),' q_cl_m =',svm(i,j,k,iq_cl) - write(6,*) ' q_cip =',q_cip(i,j,k),' n_cip =', n_cip(i,j,k),' n_ci_m =', svm(i,j,k,in_ci),' q_ci_m =',svm(i,j,k,iq_ci) - write(6,*) ' q_hrp =',q_hrp(i,j,k),' n_hrp =', n_hrp(i,j,k),' n_hr_m =', svm(i,j,k,in_hr),' q_hr_m =',svm(i,j,k,iq_hr) - write(6,*) ' q_hsp =',q_hsp(i,j,k),' n_hsp =', n_hsp(i,j,k),' n_hs_m =', svm(i,j,k,in_hs),' q_hs_m =',svm(i,j,k,iq_hs) - write(6,*) ' q_hgp =',q_hgp(i,j,k),' n_hgp =', n_hgp(i,j,k),' n_hg_m =', svm(i,j,k,in_hg),' q_hg_m =',svm(i,j,k,iq_hg) - ! - endif - enddo - enddo - enddo - endif - endif - endif ! flag_dbg - - end subroutine check_test - - - !! point out- gives output in a point - subroutine point_out(i,j,k, proc_name) - use modglobal, only : ih,i1,jh,j1,k1 - use modfields, only : sv0,svm,svp,qtp,thlp,ql0,exnf,rhof,qvsl,qvsi,qt0, thl0 - implicit none - ! logical, intent (inout) ::flag - character (len=7), intent (in) :: proc_name - real :: lim_thlp, lim_qtp, ssat, ssat_ice, uncond - integer, intent (in) :: i,j,k - - ! calculate - uncond=ql0(i,j,k)-q_cl(i,j,k) - ssat = 100.0*((qt0(i,j,k)-q_cl(i,j,k))/qvsl(i,j,k)-1.0) - ssat_ice = 100.0*((qt0(i,j,k)-q_cl(i,j,k))/qvsi(i,j,k)-1.0) - !#iceout uncond=ql0(i,j,k)-q_cl(i,j,k)-q_ci(i,j,k) - !#iceout ssat = 100.0*((qt0(i,j,k)-q_cl(i,j,k)-q_ci(i,j,k))/qvsl(i,j,k)-1.0) - !#iceout ssat_ice = 100.0*((qt0(i,j,k)-q_cl(i,j,k)-q_ci(i,j,k))/qvsi(i,j,k)-1.0) - ! ssat(i,j,k) = (100./qvsl(i,j,k))*max( qt0(i,j,k)-qvsl(i,j,k),0.0) - ! outputs - write(6,*) ' (i,j,k) = ',i,j,k,' update in: ', proc_name - write(6,*) ' uncondensed = ', uncond,' ssat = ',ssat,'% ',' ssat_i = ',ssat_ice,'% ' - write(6,*) ' qt = ',qt0(i,j,k),' thl = ', thl0(i,j,k), ' n_ccn = ', n_cc(i,j,k) - write(6,*) ' x_cl =', x_cl(i,j,k), ' n_cl =', n_cl(i,j,k),' q_cl =', q_cl (i,j,k) - write(6,*) ' x_ci =', x_ci(i,j,k), ' n_ci =', n_ci(i,j,k),' q_ci =', q_ci (i,j,k) - write(6,*) ' x_hr =', x_hr(i,j,k), ' n_hr =', n_hr(i,j,k),' q_hr =', q_hr (i,j,k) - write(6,*) ' x_hs =', x_hs(i,j,k), ' n_hs =', n_hs(i,j,k),' q_hs =', q_hs (i,j,k) - write(6,*) ' x_hg =', x_hg(i,j,k), ' n_hg =', n_hg(i,j,k),' q_hg =', q_hg (i,j,k) - write(6,*) ' q_clp =',q_clp(i,j,k),' n_clp =', n_clp(i,j,k),' n_cl_m =', svm(i,j,k,in_cl),' q_cl_m =',svm(i,j,k,iq_cl) - write(6,*) ' q_cip =',q_cip(i,j,k),' n_cip =', n_cip(i,j,k),' n_ci_m =', svm(i,j,k,in_ci),' q_ci_m =',svm(i,j,k,iq_ci) - write(6,*) ' q_hrp =',q_hrp(i,j,k),' n_hrp =', n_hrp(i,j,k),' n_hr_m =', svm(i,j,k,in_hr),' q_hr_m =',svm(i,j,k,iq_hr) - write(6,*) ' q_hsp =',q_hsp(i,j,k),' n_hsp =', n_hsp(i,j,k),' n_hs_m =', svm(i,j,k,in_hs),' q_hs_m =',svm(i,j,k,iq_hs) - write(6,*) ' q_hgp =',q_hgp(i,j,k),' n_hgp =', n_hgp(i,j,k),' n_hg_m =', svm(i,j,k,in_hg),' q_hg_m =',svm(i,j,k,iq_hg) - - - end subroutine point_out - + end subroutine check_allupdates subroutine check_ssat(flag_dbg) diff --git a/src/modlsmcrosssection.f90 b/src/modlsmcrosssection.f90 index 4c0ba834..9c6e5dab 100644 --- a/src/modlsmcrosssection.f90 +++ b/src/modlsmcrosssection.f90 @@ -74,7 +74,7 @@ subroutine initlsmcrosssection implicit none integer :: ierr - + namelist/NAMLSMCROSSSECTION/ & lcross, dtav, crossheight, crossplane @@ -166,6 +166,7 @@ end subroutine initlsmcrosssection subroutine lsmcrosssection use modglobal, only : rk3step,timee,dt_lim use modstat_nc, only : writestat_nc + use modsurfdata, only : isurf !#cgn implicit none @@ -177,9 +178,13 @@ subroutine lsmcrosssection end if tnext = tnext+idtav dt_lim = minval((/dt_lim,tnext-timee/)) - call wrtvert - call wrthorz - call wrtsurf + if (isurf==1) then ! #cgn - only if the soil model is used + call wrtvert + call wrthorz + call wrtsurf + else !#cgn - for other setups + call wrtsurf_basic + endif !#cgn end subroutine lsmcrosssection @@ -338,6 +343,98 @@ subroutine wrtsurf end subroutine wrtsurf + + !> Do the xy lsmcrosssections and dump them to file + subroutine wrtsurf_basic + use modglobal, only : imax,jmax,i1,j1,cexpnr,ifoutput,rtimee + use modsurfdata, only : Qnet, H, LE, G0, rs, ra, tskin, tendskin, & + cliq,rsveg,rssoil,Wl + use modstat_nc, only : lnetcdf, writestat_nc + implicit none + + + ! LOCAL + integer i,j + real, allocatable :: vars(:,:,:) + + + + + ! open(ifoutput,file='movh_qnet.'//cexpnr,position='append',action='write') + ! write(ifoutput,'(es12.5)') ((qnet(i,j),i=2,i1),j=2,j1) + ! close(ifoutput) + ! + ! open(ifoutput,file='movh_h.'//cexpnr,position='append',action='write') + ! write(ifoutput,'(es12.5)') ((h(i,j),i=2,i1),j=2,j1) + ! close(ifoutput) + ! + ! open(ifoutput,file='movh_le.'//cexpnr,position='append',action='write') + ! write(ifoutput,'(es12.5)') ((le(i,j),i=2,i1),j=2,j1) + ! close(ifoutput) + ! + ! open(ifoutput,file='movh_g0.'//cexpnr,position='append',action='write') + ! write(ifoutput,'(es12.5)') ((g0(i,j),i=2,i1),j=2,j1) + ! close(ifoutput) + ! + ! open(ifoutput,file='movh_rs.'//cexpnr,position='append',action='write') + ! write(ifoutput,'(es12.5)') ((rs(i,j),i=2,i1),j=2,j1) + ! close(ifoutput) + ! + ! open(ifoutput,file='movh_ra.'//cexpnr,position='append',action='write') + !write(ifoutput,'(es12.5)') ((ra(i,j),i=2,i1),j=2,j1) + ! close(ifoutput) + ! + + open(ifoutput,file='movh_tskin.'//cexpnr,position='append',action='write') + write(ifoutput,'(es12.5)') ((tskin(i,j),i=2,i1),j=2,j1) + close(ifoutput) + + ! + !open(ifoutput,file='movh_tendskin.'//cexpnr,position='append',action='write') + !write(ifoutput,'(es12.5)') ((tendskin(i,j),i=2,i1),j=2,j1) + !close(ifoutput) + ! + ! open(ifoutput,file='movh_cliq.'//cexpnr,position='append',action='write') + ! write(ifoutput,'(es12.5)') ((cliq(i,j),i=2,i1),j=2,j1) + ! close(ifoutput) + ! + ! open(ifoutput,file='movh_rsveg.'//cexpnr,position='append',action='write') + ! write(ifoutput,'(es12.5)') ((rsveg(i,j),i=2,i1),j=2,j1) + ! close(ifoutput) + ! + ! open(ifoutput,file='movh_rssoil.'//cexpnr,position='append',action='write') + ! write(ifoutput,'(es12.5)') ((rssoil(i,j),i=2,i1),j=2,j1) + ! close(ifoutput) + ! + ! open(ifoutput,file='movh_wl.'//cexpnr,position='append',action='write') + ! write(ifoutput,'(es12.5)') ((wl(i,j),i=2,i1),j=2,j1) + ! close(ifoutput) + + + if (lnetcdf) then + allocate(vars(1:imax,1:jmax,nvar3)) + vars = 0.0 + ! vars(:,:,1) = 0.0 ! qnet(2:i1,2:j1) + ! vars(:,:,2) = 0.0 ! h(2:i1,2:j1) + ! vars(:,:,3) = 0.0 ! le(2:i1,2:j1) + ! vars(:,:,4) = 0.0 ! g0(2:i1,2:j1) + vars(:,:,5) = tskin(2:i1,2:j1) + ! vars(:,:,6) = 0.0 ! tendskin(2:i1,2:j1) + ! vars(:,:,7) = 0.0 ! rs(2:i1,2:j1) + ! vars(:,:,8) = 0.0 ! ra(2:i1,2:j1) + ! vars(:,:,9) = 0.0 ! cliq(2:i1,2:j1) + ! vars(:,:,10) =0.0 ! Wl(2:i1,2:j1) + ! vars(:,:,11) = 0.0 !rssoil(2:i1,2:j1) + ! vars(:,:,12) = 0.0 ! rsveg(2:i1,2:j1) + call writestat_nc(ncid3,1,tncname3,(/rtimee/),nrec3,.true.) + call writestat_nc(ncid3,nvar3,ncname3(1:nvar3,:),vars,nrec3,imax,jmax) + deallocate(vars) + end if + + + end subroutine wrtsurf_basic + + !> Clean up when leaving the run subroutine exitlsmcrosssection use modstat_nc, only : exitstat_nc,lnetcdf diff --git a/src/modmicrodata3.f90 b/src/modmicrodata3.f90 index b11b9994..0a8dce34 100644 --- a/src/modmicrodata3.f90 +++ b/src/modmicrodata3.f90 @@ -46,6 +46,7 @@ module modmicrodata3 ! flags for different parts of microphysical processes logical :: l_sb_lim_aggr = .true. & !< whether to start with snow aggregation only after specific size + ,l_snow = .true. & !< whether to use aggregation of ice to snow ,l_sb_stickyice = .false. & !< whether the ice is sticky or the stickyness limited as in Seifert ,l_sb_clnuc_first = .true. & !< whether to first correct for nucleated cloud water, then other ,l_sb_conv_par = .true. & !< whether to use partial conversion to graupel diff --git a/src/modsurface.f90 b/src/modsurface.f90 index e0272294..67afcf68 100644 --- a/src/modsurface.f90 +++ b/src/modsurface.f90 @@ -87,7 +87,9 @@ subroutine initsurface ! Jarvis-Steward related variables rsminav, rssoilminav, LAIav, gDav, & ! Prescribed values for isurf 2, 3, 4 - z0, thls, ps, ustin, wtsurf, wqsurf, wsvsurf, & + z0, thls, ps, ustin, wtsurf, wqsurf, wsvsurf, & + ! ice surface + l_surfice, & ! Heterogeneous variables lhetero, xpatches, ypatches, land_use, loldtable, & ! AGS variables @@ -146,6 +148,7 @@ subroutine initsurface call MPI_BCAST(ps ,1,MY_REAL ,0,comm3d,mpierr) call MPI_BCAST(thls ,1,MY_REAL ,0,comm3d,mpierr) + call MPI_BCAST(l_surfice , 1, MPI_LOGICAL, 0, comm3d, mpierr) call MPI_BCAST(lhetero , 1, MPI_LOGICAL, 0, comm3d, mpierr) call MPI_BCAST(loldtable , 1, MPI_LOGICAL, 0, comm3d, mpierr) call MPI_BCAST(lrsAgs , 1, MPI_LOGICAL, 0, comm3d, mpierr) @@ -813,7 +816,11 @@ subroutine surface end do end do - call qtsurf + if(l_surfice) then + call qicesurf + else + call qtsurf + end if end if @@ -1100,6 +1107,92 @@ subroutine qtsurf end subroutine qtsurf + +!> Calculate the surface humidity assuming saturation with respect to ice + subroutine qicesurf + use modglobal, only : tmelt,bt,at,rd,rv,cp,es0,pref0,ijtot,i1,j1 & + ,esatitab,ttab + use modfields, only : qt0 + !use modsurfdata, only : rs, ra + use modmpi, only : my_real,mpierr,comm3d,mpi_sum,mpi_integer + + implicit none + real :: exner, tsurf, qsatsurf, surfwet, es, qtsl + integer :: i,j, patchx, patchy + integer :: Npatch(xpatches,ypatches), SNpatch(xpatches,ypatches) + integer :: tlonr + real :: lqts_patch(xpatches,ypatches) + real :: esi1s, tlo, thi + + patchx = 0 + patchy = 0 + + if(isurf <= 2) then + qtsl = 0. + do j = 2, j1 + do i = 2, i1 + exner = (ps / pref0)**(rd/cp) + tsurf = tskin(i,j) * exner + ! es = es0 * exp(at*(tsurf-tmelt) / (tsurf-bt)) + ! instead use tabulated values as in modthermodynamics + tlonr=int((tsurf-150.)*5.) ! tlonr=int((ttry-150.)*5.) + tlo=ttab(tlonr) + thi=ttab(tlonr+1) + esi1s =(thi-tsurf)*5.*esatitab(tlonr)+(tsurf-tlo)*5.*esatitab(tlonr+1) ! esi1=(thi-Tnr)*5.*esatitab(tlonr)+(Tnr-tlo)*5.*esatitab(thinr) + qsatsurf = rd / rv * esi1s / ps + surfwet = ra(i,j) / (ra(i,j) + rs(i,j)) + qskin(i,j) = surfwet * qsatsurf + (1. - surfwet) * qt0(i,j,1) + qtsl = qtsl + qskin(i,j) + end do + end do + + call MPI_ALLREDUCE(qtsl, qts, 1, MY_REAL, & + MPI_SUM, comm3d,mpierr) + qts = qts / ijtot + thvs = thls * (1. + (rv/rd - 1.) * qts) + + if (lhetero) then + lqts_patch = 0. + Npatch = 0 + do j = 2, j1 + do i = 2, i1 + patchx = patchxnr(i) + patchy = patchynr(j) + exner = (ps_patch(patchx,patchy) / pref0)**(rd/cp) + tsurf = tskin(i,j) * exner + ! esi1s = es0 * exp(at*(tsurf-tmelt) / (tsurf-bt)) + ! instead use tabulated values as in modthermodynamics + tlonr=int((tsurf-150.)*5.) ! tlonr=int((ttry-150.)*5.) + tlo=ttab(tlonr) + thi=ttab(tlonr+1) + ! #tb_heterosurf : here add land switch type + ! if (land_use(i,j).eq. type_ice) then + esi1s =(thi-tsurf)*5.*esatitab(tlonr)+(tsurf-tlo)*5.*esatitab(tlonr+1) + ! : else + ! : esi1s = es0 * exp(at*(tsurf-tmelt) / (tsurf-bt)) + ! #tb_heterosurf + qsatsurf = rd / rv * esi1s / ps_patch(patchx,patchy) + surfwet = ra(i,j) / (ra(i,j) + rs(i,j)) + qskin(i,j) = surfwet * qsatsurf + (1. - surfwet) * qt0(i,j,1) + + lqts_patch(patchx,patchy) = lqts_patch(patchx,patchy) + qskin(i,j) + Npatch(patchx,patchy) = Npatch(patchx,patchy) + 1 + enddo + enddo + call MPI_ALLREDUCE(lqts_patch(1:xpatches,1:ypatches), qts_patch(1:xpatches,1:ypatches),& + xpatches*ypatches, MY_REAL,MPI_SUM, comm3d,mpierr) + call MPI_ALLREDUCE(Npatch(1:xpatches,1:ypatches) , SNpatch(1:xpatches,1:ypatches) ,& + xpatches*ypatches,MPI_INTEGER ,MPI_SUM, comm3d,mpierr) + qts_patch = qts_patch / SNpatch + thvs_patch = thls_patch * (1. + (rv/rd - 1.) * qts_patch) + endif + end if + + return + + end subroutine qicesurf + + !> Calculates the Obukhov length iteratively. subroutine getobl use modglobal, only : zf, rv, rd, grav, i1, j1, i2, j2, cu, cv @@ -2207,7 +2300,11 @@ subroutine do_lsm thls_patch = thls_patch / SNpatch endif - call qtsurf + if(l_surfice) then + call qicesurf + else + call qtsurf + end if end subroutine do_lsm diff --git a/src/modsurfdata.f90 b/src/modsurfdata.f90 index 277eda05..cc4e8a87 100644 --- a/src/modsurfdata.f90 +++ b/src/modsurfdata.f90 @@ -210,6 +210,9 @@ module modsurfdata real, allocatable :: gD (:,:) !< Response factor vegetation to vapor pressure deficit [-] real :: gDav + ! Ice curface flag + logical :: l_surfice = .false. !< Flag for surface saturated with respect to ice + ! Turbulent exchange variables logical :: lmostlocal = .false. !< Switch to apply MOST locally to get local Obukhov length logical :: lsmoothflux = .false. !< Create uniform sensible and latent heat flux over domain diff --git a/src/modtestbed.f90 b/src/modtestbed.f90 index 3095605c..e35f7cac 100644 --- a/src/modtestbed.f90 +++ b/src/modtestbed.f90 @@ -34,7 +34,7 @@ module modtestbed implicit none PRIVATE -PUBLIC :: inittestbed, testbednudge, exittestbed, ltestbed, & +PUBLIC :: inittestbed, testbednudge, exittestbed, ltestbed, & testbed_getinttime, ntnudge, nknudge, & tb_time,tb_ps,tb_qts,tb_thls,tb_wqs,tb_wts, & tb_z0h, tb_z0m, tb_alb, tb_Qnet, & @@ -42,6 +42,7 @@ module modtestbed tb_dqtdxls,tb_dqtdyls, & tb_qtadv,tb_thladv,tb_uadv,tb_vadv, & tb_sv,tb_svadv,tb_svs, & !#sb3 + ltb_sv, & !#sb3 tb_tsoilav,tb_phiwav, & tbrad_p, tbrad_ql, tbrad_qv, tbrad_t, tbrad_o3 SAVE @@ -54,9 +55,11 @@ module modtestbed real, dimension(:,:), allocatable :: tb_svs !#sb3 real, dimension(:) , allocatable :: tb_time, tb_ps, tb_qts, tb_thls, tb_wqs, tb_wts, tb_z0h, tb_z0m, tb_alb, tb_Qnet real :: tb_taunudge = 10800. & - ,tb_zminnudge = 0. ! #tb + ,tb_zminnudge = 0. & !#tb + ,tb_zmidnudge = 0. !#tb logical :: ltestbed = .false., & ltb_nudge = .false., & + ltb_loadsv= .false., & ltb_sv = .false., & !#sb3 ltb_u,ltb_v,ltb_w,ltb_thl,ltb_qt integer :: nknudge,ntnudge @@ -73,19 +76,22 @@ subroutine inittestbed use modnudge, only : lnudge ! #tb implicit none - + + real, dimension(:,:,:), allocatable :: dumsv !#sb3 real, dimension(:,:), allocatable :: dumomega,dumqv,dumql,dumqi,dumt,dumpf, dumo3,& dumheight,dumqt,dumthl,dumu,dumv,dumw, & dumug,dumvg,dumqtadv,dumthladv,dumuadv,dumvadv, & dumqadv,dumladv,dumiadv,dumtadv, & dumtsoilav,dumphiwav,dumswi,& dumlwnet,dumswnet + real, dimension(:,:), allocatable :: dumsvs ! #sb3 real, dimension(:), allocatable :: dumheights real :: dumphifc,dumphiwp INTEGER NCID, STATUS, VARID, timID INTEGER start2(2), count2(2) + INTEGER start_sv(3), count_sv(3),start_svs(2), count_svs(2) character(len = nf90_max_name) :: RecordDimName integer :: ierr,i,k,ik,nknudgep1,nknudges @@ -94,7 +100,9 @@ subroutine inittestbed namelist /NAMTESTBED/ & ltestbed, ltb_nudge, tb_taunudge & ,tb_zminnudge & !#tb - ,ltb_sv !#sb3 + ,tb_zmidnudge & !#tb + ,ltb_loadsv !#sb3 + ! ltb_sv !#sb3 if(myid==0)then @@ -111,12 +119,23 @@ subroutine inittestbed ! write(6,*) " OVERRIDING, setting lnudge=.false. " endif ! #tb END - + ! check for + if (tb_zmidnudge.lt.tb_zminnudge) then + tb_zmidnudge = tb_zminnudge + endif + + if (nsv>0 .and. ltb_loadsv) then + ltb_sv = .true. + endif + + end if call MPI_BCAST(ltestbed , 1,MPI_LOGICAL,0,comm3d,mpierr) call MPI_BCAST(ltb_nudge , 1,MPI_LOGICAL,0,comm3d,mpierr) + call MPI_BCAST(tb_taunudge , 1,MY_REAL ,0,comm3d,mpierr) call MPI_BCAST(tb_zminnudge , 1,MY_REAL ,0,comm3d,mpierr) !#tb + call MPI_BCAST(tb_zmidnudge , 1,MY_REAL ,0,comm3d,mpierr) !#tb call MPI_BCAST(ltb_sv , 1,MPI_LOGICAL,0,comm3d,mpierr) !#sb3 if (.not. ltestbed) return @@ -172,8 +191,9 @@ subroutine inittestbed tb_qtadv (ntnudge,k1), & tb_thladv (ntnudge,k1), & tb_uadv (ntnudge,k1), & - tb_vadv (ntnudge,k1), & - tb_time (ntnudge), & + tb_vadv (ntnudge,k1) ) + + allocate( tb_time (ntnudge), & tb_ps (ntnudge), & tb_qts (ntnudge), & tb_thls (ntnudge), & @@ -182,18 +202,26 @@ subroutine inittestbed tb_z0m (ntnudge), & tb_z0h (ntnudge), & tb_alb (ntnudge), & - tb_Qnet (ntnudge), & - tb_tsoilav(ntnudge,ksoilmax), & - tb_phiwav (ntnudge,ksoilmax), & - tbrad_p (ntnudge, nknudge), & + tb_Qnet (ntnudge) ) + + allocate( tb_tsoilav(ntnudge,ksoilmax), & + tb_phiwav (ntnudge,ksoilmax) ) + + allocate( tbrad_p (ntnudge, nknudge), & tbrad_t (ntnudge, nknudge), & tbrad_qv (ntnudge, nknudge), & tbrad_ql (ntnudge, nknudge), & - tbrad_o3 (ntnudge, nknudge), & - tb_sv (ntnudge,k1,nsv), & ! #sb3 + tbrad_o3 (ntnudge, nknudge) ) + + if (nsv>0) then + allocate( tb_sv (ntnudge,k1,nsv), & ! #sb3 tb_svadv (ntnudge,k1,nsv), & ! #sb3 tb_svs (ntnudge,nsv) & ! #sb3 - ) + ) + tb_sv = 0 + tb_svadv = 0 + tb_svs = 0 + endif tnudge = tb_taunudge !nudging timescale @@ -262,7 +290,11 @@ subroutine inittestbed dumphiwav (nknudges,ntnudge), & dumswi (nknudges,ntnudge) & ) - + if (nsv>0 .and. ltb_loadsv) then + allocate(dumsv (nsv,nknudge,ntnudge) ) ! 3D fields + ! add dumsv_adv + allocate(dumsvs (nsv,ntnudge) ) ! surface values + endif !--- timeseries --- @@ -537,8 +569,38 @@ subroutine inittestbed if (status /= nf90_noerr) call handle_err(status) dumswi = ( dumphiwav - dumphiwp ) / ( dumphifc - dumphiwp ) !soil wetness index, using input values for wilting point and field capacity + + ! loading scalars + if (nsv>0 .and. ltb_loadsv) then + + ! ltb_sv = .true. <- done before + start_svs = (/ 1 , 1 /) + count_svs = (/ nsv , ntnudge /) + start_sv = (/ 1 , 1 , 1 /) + count_sv = (/ nsv , nknudge , ntnudge /) + ! other options: + !2) + ! count_sv = (/ nsv , ntnudge , nknudge /) + !3) + ! count_sv = (/ nknudge , ntnudge , nsv /) + ! + ! surface timeseries + write(6,*) 'Now loading svs variables' + write(6,*) count_svs + STATUS = NF90_INQ_VARID(NCID, 'sv_flx', VARID) + if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + STATUS = NF90_GET_VAR (NCID, VARID, dumsvs, start=start_svs, count=count_svs) + if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + write(6,*) 'Now loading sv variables' + write(6,*) count_sv + ! sv profiles for nsv species + STATUS = NF90_INQ_VARID(NCID, 'sv', VARID) + if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + STATUS = NF90_GET_VAR (NCID, VARID, dumsv, start=start_sv, count=count_sv) + if (STATUS .ne. nf90_noerr) call handle_err(STATUS) + endif !--- close nc file --- STATUS = NF90_CLOSE(NCID) @@ -575,7 +637,18 @@ subroutine inittestbed tb_thladv (i,k) = dumthladv (ik,i) + fac * ( dumthladv (ik-1,i) - dumthladv (ik,i) ) !if (i.eq.1) write(6,*) k, zf(k), " : ", ik, dumheight(ik,i), ik-1, dumheight(ik-1,i) - + ! and for scalars + + ! if (ltb_sv .and. i.eq.1) then + if (ltb_sv) then + !tb_sv(k,:) = dumsv(:,ik0) + fac * ( dumsv(:,ik1) - dumsv(:,ik0) ) + tb_sv(i,k,:) = dumsv(:,ik,i) + fac * ( dumsv(:,ik-1,i) - dumsv(:,ik,i) ) + ! outputs + ! write(0,'(a,i8,f8.2,i8,f15.7,i8,2f15.7,f20.3)') 'inittestbed: sv field loaded from scm_in.nc: ', k, zf(k), & + ! ik0, dumheight(ik0,i), ik1, dumheight(ik1,i), fac, tb_sv(k,6) + ! write(0,'(a,i8,f8.2,i8,f15.7,i8,2f15.7,f20.3)') 'inittestbed: sv field loaded from scm_in.nc: ', k, zf(k), & + ! ik, dumheight(ik,i), ik-1, dumheight(ik-1,i), fac, tb_sv(i,k,6) + endif enddo ! @@ -646,7 +719,11 @@ subroutine inittestbed deallocate(dumtsoilav) deallocate(dumphiwav) deallocate(dumswi) - + + if (ltb_sv) then + deallocate(dumsv) + deallocate(dumsvs) + endif !--- do some output to screen --- ! do i=1,2 @@ -710,8 +787,11 @@ subroutine inittestbed call MPI_BCAST(tbrad_t ,ntnudge*nknudge,MY_REAL ,0,comm3d,mpierr) call MPI_BCAST(tbrad_o3 ,ntnudge*nknudge,MY_REAL ,0,comm3d,mpierr) - call MPI_BCAST(tb_sv ,ntnudge*k1*nsv,MY_REAL ,0,comm3d,mpierr) - + if (nsv>0) then + call MPI_BCAST(tb_sv ,ntnudge*k1*nsv ,MY_REAL ,0,comm3d,mpierr) + call MPI_BCAST(tb_svs ,ntnudge*nsv ,MY_REAL ,0,comm3d,mpierr) + endif + ltb_u = any(abs(tb_u)>1e-8) ltb_v = any(abs(tb_v)>1e-8) ltb_w = any(abs(tb_w)>1e-8) @@ -729,7 +809,7 @@ subroutine testbednudge implicit none integer k,t - real :: dtm,dtp,currtnudge, qttnudge,qtthres + real :: dtm,dtp,currtnudge, qttnudge,qtthres, zfac, zfacu if (.not.(ltestbed .and. ltb_nudge)) return @@ -749,20 +829,33 @@ subroutine testbednudge qtthres = 1e-6 do k=1,kmax - if ( zf(k).gt. tb_zminnudge) then ! tb_zminnudge #tb START + zfac = 0. + if ( zf(k).gt. tb_zminnudge) then + ! set factor in nudging intensity: linear increase from 0 to 1 in heightrange tb_zmidnudge > z > tb_zminnudge + zfac = zf(k) - tb_zminnudge + if (tb_zmidnudge.gt.tb_zminnudge) then + zfac = max( 0., min( 1., zfac / (tb_zmidnudge - tb_zminnudge) ) ) + else + zfac = 0.5 + sign( 0.5, zfac ) + endif + !write(0,*) 'modtestbed: ', k, zf(k), tb_zminnudge, tb_zmidnudge, zfac + endif + + zfacu = zfac + !zfacu = 1. currtnudge = max(rdt,tnudge(t,k)*dtp+tnudge(t+1,k)*dtm) - if (ltb_u) up(2:i1,2:j1,k) = up(2:i1,2:j1,k) - & + if (ltb_u) up(2:i1,2:j1,k) = up(2:i1,2:j1,k) - zfacu * & ( u0av(k) - (tb_u(t,k) *dtp + tb_u(t+1,k) *dtm) ) / currtnudge - if (ltb_v) vp(2:i1,2:j1,k) = vp(2:i1,2:j1,k) - & + if (ltb_v) vp(2:i1,2:j1,k) = vp(2:i1,2:j1,k) - zfacu * & ( v0av(k) - (tb_v(t,k) *dtp + tb_v(t+1,k) *dtm) ) / currtnudge - if (ltb_w) wp(2:i1,2:j1,k) = wp(2:i1,2:j1,k) - & + if (ltb_w) wp(2:i1,2:j1,k) = wp(2:i1,2:j1,k) - zfac * & ( - (tb_w(t,k) *dtp + tb_w(t+1,k) *dtm) ) / currtnudge - if (ltb_thl) thlp(2:i1,2:j1,k) = thlp(2:i1,2:j1,k) - & + if (ltb_thl) thlp(2:i1,2:j1,k) = thlp(2:i1,2:j1,k) -zfac * & ( thl0av(k) - (tb_thl(t,k)*dtp + tb_thl(t+1,k)*dtm) ) / currtnudge if (ltb_qt) then @@ -774,7 +867,7 @@ subroutine testbednudge qtp(2:i1,2:j1,k) = qtp(2:i1,2:j1,k) - & ( qt0av(k) - (tb_qt(t,k) *dtp + tb_qt(t+1,k) *dtm) ) / qttnudge end if - endif ! tb_zminnudge #tb END + !endif ! tb_zminnudge #tb END end do ! if(ltb_sv.and. ltb_nudge) then @@ -816,7 +909,8 @@ end subroutine testbed_getinttime !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine exittestbed - if (allocated(tb_time)) then + use modglobal, only: nsv + if (allocated(tb_time)) then deallocate( tnudge & ,tb_u & ,tb_v & @@ -848,11 +942,14 @@ subroutine exittestbed ,tbrad_qv & ,tbrad_ql & ,tbrad_o3 & - ,tb_sv & + ) + if(nsv>0) then + deallocate(tb_sv & ,tb_svadv & ,tb_svs & - ) - end if + ) + end if + end if end subroutine exittestbed !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!