!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! GNU General Public License !! !! !! !! This file is part of the Flexible Modeling System (FMS). !! !! !! !! FMS is free software; you can redistribute it and/or modify !! !! it and are expected to follow the terms of the GNU General Public !! !! License as published by the Free Software Foundation. !! !! !! !! FMS is distributed in the hope that it will be useful, !! !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! !! GNU General Public License for more details. !! !! !! !! You should have received a copy of the GNU General Public License !! !! along with FMS; if not, write to: !! !! Free Software Foundation, Inc. !! !! 59 Temple Place, Suite 330 !! !! Boston, MA 02111-1307 USA !! !! or see: !! !! http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! MODULE MY25_TURB_MOD !======================================================================= ! MELLOR-YAMADA LEVEL 2.5 TURBULENCE CLOSURE SCHEME - GFDL VERSION ! !======================================================================= use Fms_Mod, ONLY: FILE_EXIST, OPEN_NAMELIST_FILE, ERROR_MESG, FATAL, & CLOSE_FILE, NOTE, READ_DATA, & check_nml_error, mpp_pe, mpp_root_pe, & write_version_number, stdlog, open_restart_file use fms_io_mod, only: register_restart_field, restart_file_type use fms_io_mod, only: save_restart, restore_state use Tridiagonal_Mod, ONLY: TRI_INVERT, CLOSE_TRIDIAGONAL use constants_mod, only: grav, vonkarm use monin_obukhov_mod, only : mo_diff !--------------------------------------------------------------------- implicit none private !--------------------------------------------------------------------- public :: MY25_TURB, MY25_TURB_INIT, MY25_TURB_END, TKE_SURF, get_tke public :: my25_turb_restart !--------------------------------------------------------------------- ! --- GLOBAL STORAGE !--------------------------------------------------------------------- real, allocatable, dimension(:,:,:) :: TKE !--------------------------------------------------------------------- character(len=128) :: version = '$Id: my25_turb.F90,v 17.0.2.1 2009/09/24 15:28:34 z1l Exp $' character(len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $' logical :: module_is_initialized = .false. logical :: init_tke integer :: num_total_pts, pts_done ! for netcdf restart file. type(restart_file_type), save :: Tur_restart !--------------------------------------------------------------------- ! --- CONSTANTS !--------------------------------------------------------------------- real :: aa1, aa2, bb1, bb2, ccc real :: ckm1, ckm2, ckm3, ckm4, ckm5, ckm6, ckm7, ckm8 real :: ckh1, ckh2, ckh3, ckh4 real :: cvfq1, cvfq2, bcq real, parameter :: aa1_old = 0.78 real, parameter :: aa2_old = 0.79 real, parameter :: bb1_old = 15.0 real, parameter :: bb2_old = 8.0 real, parameter :: ccc_old = 0.056 real, parameter :: aa1_new = 0.92 real, parameter :: aa2_new = 0.74 real, parameter :: bb1_new = 16.0 real, parameter :: bb2_new = 10.0 real, parameter :: ccc_new = 0.08 real, parameter :: cc1 = 0.27 real, parameter :: t00 = 2.7248e2 real, parameter :: small = 1.0e-10 !--------------------------------------------------------------------- ! --- NAMELIST !--------------------------------------------------------------------- real :: TKEmax = 5.0 real :: TKEmin = 0.0 real :: el0max = 1.0e6 real :: el0min = 0.0 real :: alpha_land = 0.10 real :: alpha_sea = 0.10 real :: akmax = 1.0e4 real :: akmin_land = 5.0 real :: akmin_sea = 0.0 integer :: nk_lim = 2 integer :: init_iters = 20 logical :: do_thv_stab = .true. logical :: use_old_cons = .false. real :: kcrit = 0.01 NAMELIST / my25_turb_nml / & TKEmax, TKEmin, init_iters, & akmax, akmin_land, akmin_sea, nk_lim, & el0max, el0min, alpha_land, alpha_sea, & do_thv_stab, use_old_cons, & kcrit !--------------------------------------------------------------------- contains !####################################################################### SUBROUTINE MY25_TURB( is, js, delt, fracland, phalf, pfull, theta, & um, vm, zhalf, zfull, z0, & el0, el, akm, akh, & mask, kbot, ustar, bstar, h ) !======================================================================= !--------------------------------------------------------------------- ! Arguments (Intent in) ! delt - Time step in seconds ! fracland - Fractional amount of land beneath a grid box ! phalf - Pressure at half levels ! pfull - Pressure at full levels ! theta - Potential temperature ! um, vm - Wind components ! zhalf - Height at half levels ! zfull - Height at full levels ! z0 - Roughness length ! mask - OPTIONAL; floating point mask (0. or 1.) designating ! where data is present ! kbot - OPTIONAL;lowest model level index (integer); ! at levels > kbot, mask = 0. ! ustar - OPTIONAL:friction velocity (m/sec) ! bstar - OPTIONAL:buoyancy scale (m/sec**2) !--------------------------------------------------------------------- integer, intent(in) :: is, js real, intent(in) :: delt real, intent(in), dimension(:,:) :: fracland, z0 real, intent(in), dimension(:,:,:) :: phalf, pfull, zhalf, zfull real, intent(in), dimension(:,:,:) :: um, vm, theta integer, intent(in), OPTIONAL, dimension(:,:) :: kbot real, intent(in), OPTIONAL, dimension(:,:,:) :: mask real, intent(in), OPTIONAL, dimension(:,:) :: ustar, bstar !--------------------------------------------------------------------- ! Arguments (Intent out) ! el0 - characteristic length scale ! el - master length scale ! akm - mixing coefficient for momentum ! akh - mixing coefficient for heat and moisture ! h - OPTIONAL, diagnosed depth of planetary boundary ! layer (m) !--------------------------------------------------------------------- real, intent(out), dimension(:,:) :: el0 real, intent(out), dimension(:,:,:) :: akm, akh, el real, intent(out), OPTIONAL, dimension(:,:) :: h !--------------------------------------------------------------------- ! (Intent local) !--------------------------------------------------------------------- integer :: ix, jx, kx, i, j, k, ism, jsm integer :: kxp, kxm, klim, it, itermax real :: cvfqdt, dvfqdt real, dimension(SIZE(um,1),SIZE(um,2)) :: zsfc, x1, x2, akmin real, dimension(SIZE(um,1),SIZE(um,2),SIZE(um,3)-1) :: & dsdzh, shear, buoync, qm2, qm3, qm4, el2, & aaa, bbb, ccc, ddd, & xxm1, xxm2, xxm3, xxm4, xxm5 real, dimension(SIZE(um,1),SIZE(um,2),SIZE(um,3)) :: & dsdz, qm, xx1, xx2 !==================================================================== ! --- Check to see if MY25_TURB has been initialized if( .not. module_is_initialized ) CALL ERROR_MESG( ' MY25_TURB', & ' MY25_TURB_INIT has not been called',& FATAL ) ! --- Set dimensions etc ix = SIZE( um, 1 ) jx = SIZE( um, 2 ) kx = SIZE( um, 3 ) kxp = kx + 1 kxm = kx - 1 ism = is - 1 jsm = js - 1 !==================================================================== ! --- SURFACE HEIGHT !==================================================================== if( PRESENT( kbot ) ) then do j = 1,jx do i = 1,ix k = kbot(i,j) + 1 zsfc(i,j) = zhalf(i,j,k) end do end do else zsfc(:,:) = zhalf(:,:,kxp) endif !==================================================================== ! --- D( )/DZ OPERATORS: AT FULL LEVELS & AT HALF LEVELS !==================================================================== dsdz(:,:,1:kx) = 1.0 / ( zhalf(:,:,2:kxp) - zhalf(:,:,1:kx) ) dsdzh(:,:,1:kxm) = 1.0 / ( zfull(:,:,2:kx) - zfull(:,:,1:kxm) ) !==================================================================== ! --- WIND SHEAR !==================================================================== xxm1(:,:,1:kxm) = dsdzh(:,:,1:kxm)*( um(:,:,2:kx) - um(:,:,1:kxm) ) xxm2(:,:,1:kxm) = dsdzh(:,:,1:kxm)*( vm(:,:,2:kx) - vm(:,:,1:kxm) ) shear = xxm1 * xxm1 + xxm2 * xxm2 !==================================================================== ! --- BUOYANCY !==================================================================== xxm1(:,:,1:kxm) = theta(:,:,2:kx) - theta(:,:,1:kxm) if( do_thv_stab ) then xxm2(:,:,1:kxm) = 0.5*( theta(:,:,2:kx) + theta(:,:,1:kxm) ) else xxm2(:,:,1:kxm) = t00 end if buoync = grav * dsdzh * xxm1 / xxm2 !==================================================================== ! --- MASK OUT UNDERGROUND VALUES FOR ETA COORDINATE !==================================================================== if( PRESENT( mask ) ) then do k=2,kx do j=1,jx do i=1,ix if(mask(i,j,k) < 0.1) then TKE(ism+i,jsm+j,k+1) = 0.0 dsdz(i,j,k ) = 0.0 dsdzh(i,j,k-1) = 0.0 shear(i,j,k-1) = 0.0 buoync(i,j,k-1) = 0.0 endif enddo enddo enddo endif !==================================================================== ! --- SET ITERATION LOOP IF INITALIZING TKE !==================================================================== itermax = 1 if (init_tke) then itermax = init_iters pts_done = pts_done + ix*jx if (pts_done >= num_total_pts) init_tke = .false. endif ! $$$$$$$$$$$$$$$$$ do it = 1,itermax ! $$$$$$$$$$$$$$$$$ !==================================================================== ! --- SOME TKE STUFF !==================================================================== do k=1,kx do j=1,jx do i=1,ix xx1(i,j,k) = 2*TKE(ism+i,jsm+j,k+1) if(xx1(i,j,k) > 0.0) then qm(i,j,k) = sqrt(xx1(i,j,k)) else qm(i,j,k) = 0.0 endif enddo enddo enddo qm2(:,:,1:kxm) = xx1(:,:,1:kxm) qm3(:,:,1:kxm) = qm(:,:,1:kxm) * qm2(:,:,1:kxm) qm4(:,:,1:kxm) = qm2(:,:,1:kxm) * qm2(:,:,1:kxm) !==================================================================== ! --- CHARACTERISTIC LENGTH SCALE !==================================================================== xx1(:,:,1:kxm) = qm(:,:,1:kxm)*( pfull(:,:,2:kx) - pfull(:,:,1:kxm) ) do k = 1, kxm xx2(:,:,k) = xx1(:,:,k) * ( zhalf(:,:,k+1) - zsfc(:,:) ) end do if( PRESENT( kbot ) ) then xx1(:,:,kx) = 0.0 xx2(:,:,kx) = 0.0 do j = 1,jx do i = 1,ix k = kbot(i,j) xx1(i,j,k) = qm(i,j,k) * ( phalf(i,j,k+1) - pfull(i,j,k) ) xx2(i,j,k) = xx1(i,j,k) * z0(i,j) end do end do else xx1(:,:,kx) = qm(:,:,kx) * ( phalf(:,:,kxp) - pfull(:,:,kx) ) xx2(:,:,kx) = xx1(:,:,kx) * z0(:,:) endif if (PRESENT(mask)) then x1 = SUM( xx1, 3, mask=mask.gt.0.1 ) x2 = SUM( xx2, 3, mask=mask.gt.0.1 ) else x1 = SUM( xx1, 3 ) x2 = SUM( xx2, 3 ) endif !---- should never be equal to zero ---- if (count(x1 <= 0.0) > 0) CALL ERROR_MESG( ' MY25_TURB', & 'divid by zero, x1 <= 0.0', FATAL) el0 = x2 / x1 el0 = el0 * (alpha_land*fracland + alpha_sea*(1.-fracland)) el0 = MIN( el0, el0max ) el0 = MAX( el0, el0min ) !==================================================================== ! --- MASTER LENGTH SCALE !==================================================================== do k = 1, kxm xx1(:,:,k) = vonkarm * ( zhalf(:,:,k+1) - zsfc(:,:) ) end do x1(:,:) = vonkarm * z0(:,:) if( PRESENT( kbot ) ) then do j = 1,jx do i = 1,ix do k = kbot(i,j), kx xx1(i,j,k) = x1(i,j) end do end do end do else xx1(:,:,kx) = x1(:,:) endif do k = 1,kx el(:,:,k+1) = xx1(:,:,k) / ( 1.0 + xx1(:,:,k) / el0(:,:) ) end do el(:,:,1) = el0(:,:) el2(:,:,1:kxm) = el(:,:,2:kx) * el(:,:,2:kx) !==================================================================== ! --- MIXING COEFFICIENTS !==================================================================== xxm3(:,:,1:kxm) = el2(:,:,1:kxm)*buoync(:,:,1:kxm) xxm4(:,:,1:kxm) = el2(:,:,1:kxm)* shear(:,:,1:kxm) xxm5(:,:,1:kxm) = el(:,:,2:kx )* qm3(:,:,1:kxm) !------------------------------------------------------------------- ! --- MOMENTUM !------------------------------------------------------------------- xxm1 = xxm5*( ckm1*qm2 + ckm2*xxm3 ) xxm2 = qm4 + ckm5*qm2*xxm4 + xxm3*( ckm6*xxm4 + ckm7*qm2 + ckm8*xxm3 ) xxm2 = MAX( xxm2, 0.2*qm4 ) xxm2 = MAX( xxm2, small ) akm(:,:,1) = 0.0 akm(:,:,2:kx) = xxm1(:,:,1:kxm) / xxm2(:,:,1:kxm) akm = MAX( akm, 0.0 ) !------------------------------------------------------------------- ! --- HEAT AND MOISTURE !------------------------------------------------------------------- xxm1(:,:,1:kxm) = ckh1*xxm5(:,:,1:kxm) - ckh2*xxm4(:,:,1:kxm)*akm(:,:,2:kx) xxm2(:,:,1:kxm) = qm2(:,:,1:kxm) + ckh3*xxm3(:,:,1:kxm) xxm1 = MAX( xxm1, ckh4*xxm5 ) xxm2 = MAX( xxm2, 0.4*qm2 ) xxm2 = MAX( xxm2, small ) akh(:,:,1) = 0.0 akh(:,:,2:kx) = xxm1(:,:,1:kxm) / xxm2(:,:,1:kxm) !------------------------------------------------------------------- ! --- BOUNDS !------------------------------------------------------------------- ! --- UPPER BOUND akm = MIN( akm, akmax ) akh = MIN( akh, akmax ) ! --- LOWER BOUND ! where( akm(:,:,1:klim) < small ) akm(:,:,1:klim) = 0.0 ! where( akh(:,:,1:klim) < small ) akh(:,:,1:klim) = 0.0 ! --- LOWER BOUND NEAR SURFACE akmin = akmin_land*fracland + akmin_sea*(1.-fracland) if( PRESENT( kbot ) ) then do j = 1,jx do i = 1,ix klim = kbot(i,j) - nk_lim + 1 do k = klim,kbot(i,j) akm(i,j,k) = MAX( akm(i,j,k), akmin(i,j) ) akh(i,j,k) = MAX( akh(i,j,k), akmin(i,j) ) end do end do end do else klim = kx - nk_lim + 1 do k = klim,kx akm(:,:,k) = MAX( akm(:,:,k), akmin(:,:) ) akh(:,:,k) = MAX( akh(:,:,k), akmin(:,:) ) end do endif !------------------------------------------------------------------- ! --- MASK OUT UNDERGROUND VALUES FOR ETA COORDINATE !------------------------------------------------------------------- if( PRESENT( mask ) ) then akm(:,:,1:kx) = akm(:,:,1:kx) * mask(:,:,1:kx) akh(:,:,1:kx) = akh(:,:,1:kx) * mask(:,:,1:kx) endif !==================================================================== ! --- PROGNOSTICATE TURBULENT KE !==================================================================== cvfqdt = cvfq1 * delt dvfqdt = cvfq2 * delt * 2.0 !------------------------------------------------------------------- ! --- PART OF LINEARIZED ENERGY DISIIPATION TERM !------------------------------------------------------------------- xxm1(:,:,1:kxm) = dvfqdt * qm(:,:,1:kxm) / el(:,:,2:kx) !------------------------------------------------------------------- ! --- PART OF LINEARIZED VERTICAL DIFFUSION TERM !------------------------------------------------------------------- xx1(:,:,1:kx) = el(:,:,2:kxp) * qm(:,:,1:kx) xx2(:,:,1) = 0.5* xx1(:,:,1) xx2(:,:,2:kx) = 0.5*( xx1(:,:,2:kx) + xx1(:,:,1:kxm) ) xx1 = xx2 * dsdz !------------------------------------------------------------------- ! --- IMPLICIT TIME DIFFERENCING FOR VERTICAL DIFFUSION ! --- AND ENERGY DISSIPATION TERM !------------------------------------------------------------------- do k=1,kxm do j=1,jx do i=1,ix aaa(i,j,k) = -cvfqdt * xx1(i,j,k+1) * dsdzh(i,j,k) ccc(i,j,k) = -cvfqdt * xx1(i,j,k ) * dsdzh(i,j,k) bbb(i,j,k) = 1.0 - aaa(i,j,k ) - ccc(i,j,k) bbb(i,j,k) = bbb(i,j,k ) + xxm1(i,j,k) ddd(i,j,k) = TKE(ism+i,jsm+j,k+1) enddo enddo enddo ! correction for vertical diffusion of TKE surface boundary condition if (present(kbot)) then do j = 1,jx do i = 1,ix k = kbot(i,j) ddd(i,j,k-1) = ddd(i,j,k-1) - aaa(i,j,k-1) * TKE(ism+i,jsm+j,k+1) enddo enddo else do j = 1,jx do i = 1,ix ddd(i,j,kxm) = ddd(i,j,kxm) - aaa(i,j,kxm) * TKE(ism+i,jsm+j,kxp) enddo enddo endif ! mask out terms below ground if (present(mask)) then where (mask(:,:,2:kx) < 0.1) ddd(:,:,1:kxm) = 0.0 endif CALL TRI_INVERT( xxm1, ddd, aaa, bbb, ccc ) CALL CLOSE_TRIDIAGONAL !------------------------------------------------------------------- ! --- MASK OUT UNDERGROUND VALUES FOR ETA COORDINATE !------------------------------------------------------------------- if( PRESENT( mask ) ) then do k=1,kxm do j=1,jx do i=1,ix if(mask(i,j,k+1) < 0.1) xxm1(i,j,k) = TKE(ism+i,jsm+j,k+1) enddo enddo enddo endif !------------------------------------------------------------------- ! --- SHEAR AND BUOYANCY TERMS !------------------------------------------------------------------- xxm2(:,:,1:kxm) = delt*( akm(:,:,2:kx)* shear(:,:,1:kxm) & - akh(:,:,2:kx)*buoync(:,:,1:kxm) ) !------------------------------------------------------------------- ! --- UPDATE TURBULENT KINETIC ENERGY !------------------------------------------------------------------- do j=1,jx do i=1,ix TKE(ism+i,jsm+j,1) = 0.0 do k=2,kx TKE(ism+i,jsm+j,k) = xxm1(i,j,k-1) + xxm2(i,j,k-1) enddo enddo enddo !==================================================================== ! --- BOUND TURBULENT KINETIC ENERGY !==================================================================== TKE(is:ism+ix,js:jsm+jx,:) = MIN( TKE(is:ism+ix,js:jsm+jx,:), TKEmax ) TKE(is:ism+ix,js:jsm+jx,:) = MAX( TKE(is:ism+ix,js:jsm+jx,:), TKEmin ) if( PRESENT( mask ) ) then do k=1,kx do j=1,jx do i=1,ix if(mask(i,j,k) < 0.1) TKE(ism+i,jsm+j,k+1) = 0.0 enddo enddo enddo endif !==================================================================== ! --- COMPUTE PBL DEPTH IF DESIRED !==================================================================== if (present(h)) then if (.not.present(ustar).or..not.present(bstar)) then CALL ERROR_MESG( ' MY25_TURB', & 'cannot request pbl depth diagnostic if ustar'// & ' and bstar are not also supplied', FATAL ) end if if (present(kbot)) then call k_pbl_depth(ustar,bstar,akm,akh,zsfc,zfull,zhalf,& h,kbot=kbot) else call k_pbl_depth(ustar,bstar,akm,akh,zsfc,zfull,zhalf,h) end if end if !==================================================================== ! $$$$$$$$$$$$$$$$$ end do ! $$$$$$$$$$$$$$$$$ !==================================================================== end SUBROUTINE MY25_TURB !####################################################################### subroutine get_tke(is, ie, js, je, tke_out) integer, intent(in) :: is, ie, js, je real, intent(out), dimension(:,:,:) :: tke_out if( .not. module_is_initialized ) then CALL ERROR_MESG(' MY25_TURB','MY25_TURB_INIT has not been called',FATAL) endif tke_out = TKE(is:ie,js:je,:) end subroutine get_tke !####################################################################### SUBROUTINE MY25_TURB_INIT( ix, jx, kx ) !======================================================================= ! ***** INITIALIZE MELLOR-YAMADA !======================================================================= !--------------------------------------------------------------------- ! Arguments (Intent in) ! ix, jx - Horizontal dimensions for global storage arrays ! kx - Number of vertical levels in model !--------------------------------------------------------------------- integer, intent(in) :: ix, jx, kx !--------------------------------------------------------------------- ! (Intent local) !--------------------------------------------------------------------- integer :: unit, io, ierr, logunit integer :: id_restart !===================================================================== !--------------------------------------------------------------------- ! --- Read namelist !--------------------------------------------------------------------- if( FILE_EXIST( 'input.nml' ) ) then ! ------------------------------------- unit = OPEN_NAMELIST_FILE ( ) ierr = 1 do while( ierr .ne. 0 ) READ ( unit, nml = my25_turb_nml, iostat = io, end = 10 ) ierr = check_nml_error (io, 'my25_turb_nml') end do 10 continue CALL CLOSE_FILE( unit ) ! ------------------------------------- end if !--------------------------------------------------------------------- ! --- Output version !--------------------------------------------------------------------- if ( mpp_pe() == mpp_root_pe() ) then call write_version_number(version, tagname) logunit = stdlog() WRITE( logunit, nml = my25_turb_nml ) endif !--------------------------------------------------------------------- ! --- Initialize constants !--------------------------------------------------------------------- if( use_old_cons ) then aa1 = aa1_old aa2 = aa2_old bb1 = bb1_old bb2 = bb2_old ccc = ccc_old else aa1 = aa1_new aa2 = aa2_new bb1 = bb1_new bb2 = bb2_new ccc = ccc_new end if ckm1 = ( 1.0 - 3.0*ccc )*aa1 ckm3 = 3.0 * aa1*aa2* ( bb2 - 3.0*aa2 ) ckm4 = 9.0 * aa1*aa2*ccc*( bb2 + 4.0*aa1 ) ckm5 = 6.0 * aa1*aa1 ckm6 = 18.0 * aa1*aa1*aa2*( bb2 - 3.0*aa2 ) ckm7 = 3.0 * aa2* ( bb2 + 7.0*aa1 ) ckm8 = 27.0 * aa1*aa2*aa2*( bb2 + 4.0*aa1 ) ckm2 = ckm3 - ckm4 ckh1 = aa2 ckh2 = 6.0 * aa1*aa2 ckh3 = 3.0 * aa2*( bb2 + 4.0*aa1 ) ckh4 = 2.0e-6 * aa2 cvfq1 = 5.0 * cc1 / 3.0 cvfq2 = 1.0 / bb1 bcq = 0.5 * ( bb1**(2.0/3.0) ) !--------------------------------------------------------------------- ! --- Allocate storage for TKE !--------------------------------------------------------------------- if( ALLOCATED( TKE ) ) DEALLOCATE( TKE ) ALLOCATE( TKE(ix,jx,kx+1) ) ; TKE = TKEmin !--------------------------------------------------------------------- ! --- Input TKE !--------------------------------------------------------------------- id_restart = register_restart_field(Tur_restart, 'my25_turb.res', 'TKE', TKE) if (file_exist( 'INPUT/my25_turb.res.nc' )) then if (mpp_pe() == mpp_root_pe() ) then call error_mesg ('my25_turb_mod', 'MY25_TURB_INIT:& &Reading netCDF formatted restart file: & &INPUT/my25_turb.res.nc', NOTE) endif call restore_state(Tur_restart) else if( FILE_EXIST( 'INPUT/my25_turb.res' ) ) then unit = OPEN_restart_FILE ( file = 'INPUT/my25_turb.res', action = 'read' ) call read_data ( unit, TKE ) CALL CLOSE_FILE( unit ) init_tke = .false. else TKE = TKEmin init_tke = .true. num_total_pts = ix*jx pts_done = 0 endif !------------------------------------------------------------------- module_is_initialized = .true. !--------------------------------------------------------------------- !===================================================================== end SUBROUTINE MY25_TURB_INIT !####################################################################### SUBROUTINE MY25_TURB_END !======================================================================= integer :: unit !======================================================================= !-------------------------------------------------------------------- ! local variables: call my25_turb_restart module_is_initialized = .false. !--------------------------------------------------------------------- ! unit = OPEN_RESTART_FILE ( file = 'RESTART/my25_turb.res', action = 'write' ) ! call write_data ( unit, TKE ) ! CALL CLOSE_FILE ( unit ) !===================================================================== end SUBROUTINE MY25_TURB_END !####################################################################### ! ! ! ! write out restart file. ! Arguments: ! timestamp (optional, intent(in)) : A character string that represents the model time, ! used for writing restart. timestamp will append to ! the any restart file name as a prefix. ! ! subroutine my25_turb_restart(timestamp) character(len=*), intent(in), optional :: timestamp if(.not. present(timestamp)) then if (mpp_pe() == mpp_root_pe() ) & call error_mesg ('my25_turb_mod', 'my25_turb_end: & &Writing netCDF formatted restart file as & &requested: RESTART/my25_turb.res.nc', NOTE) endif !---------------------------------------------------------------------- ! write out the restart data which is always needed, regardless of ! when the first donner calculation step is after restart. !---------------------------------------------------------------------- call save_restart(Tur_restart, timestamp) end subroutine my25_turb_restart ! NAME="my25_turb_restart" !####################################################################### SUBROUTINE K_PBL_DEPTH(ustar,bstar,akm,akh,zsfc,zfull,zhalf,h,kbot) !======================================================================= real, intent(in), dimension(:,:) :: ustar, bstar, zsfc real, intent(in), dimension(:,:,:) :: zhalf, zfull real, intent(in), dimension(:,:,:) :: akm, akh real, intent(out),dimension(:,:) :: h integer, intent(in), optional, dimension(:,:) :: kbot !======================================================================= real, dimension(size(zfull,1),size(zfull,2)) :: km_surf, kh_surf real, dimension(size(zfull,1),size(zfull,2)) :: zhalfhalf real, dimension(size(zfull,1),size(zfull,2),size(zfull,3)) :: zfull_ag real, dimension(size(zfull,1),size(zfull,2),size(zfull,3)+1):: zhalf_ag real, dimension(size(zfull,1),size(zfull,2),size(zfull,3)+1):: diff_tm integer, dimension(size(zfull,1),size(zfull,2)) :: ibot integer :: i,j,k,nlev,nlat,nlon nlev = size(zfull,3) nlat = size(zfull,2) nlon = size(zfull,1) !compute height of surface if (present(kbot)) then ibot=kbot else ibot(:,:) = nlev end if !compute density profile, and heights relative to surface do k = 1, nlev zfull_ag(:,:,k) = zfull(:,:,k) - zsfc(:,:) zhalf_ag(:,:,k) = zhalf(:,:,k) - zsfc(:,:) end do zhalf_ag(:,:,nlev+1) = zhalf(:,:,nlev+1) - zsfc(:,:) !compute height half way between surface and lowest model level zhalfhalf=0.5*MINVAL(MAX(zfull_ag,0.0)) !compute k's there by a call to mo_diff call mo_diff(zhalfhalf,ustar,bstar,km_surf,kh_surf) !create combined surface k's and diffusivity matrix diff_tm(:,:,nlev+1) = 0. diff_tm(:,:,1:nlev) = 0.5*(akm+akh) if (present(kbot)) then do j=1,nlat do i=1,nlon diff_tm(i,j,ibot(i,j)+1) = 0.5*(km_surf(i,j)+kh_surf(i,j)) zhalf_ag(i,j,ibot(i,j)+1) = zhalfhalf(i,j) enddo enddo else diff_tm(:,:,nlev+1) = 0.5*(km_surf(:,:)+kh_surf(:,:)) zhalf_ag(:,:,nlev+1) = zhalfhalf(:,:) end if !determine pbl depth as the height above ground where diff_tm !first falls beneath a critical value kcrit. If the value between !ground and level 1 does not exceed kcrit set pbl depth equal to zero. !kcrit is a namelist parameter. do j = 1,nlat do i = 1,nlon if (diff_tm(i,j,ibot(i,j)+1) .gt. kcrit) then k=ibot(i,j)+1 do while (k.gt. 2 .and. & diff_tm(i,j,k-1).gt.kcrit) k=k-1 enddo h(i,j) = zhalf_ag(i,j,k) + & (zhalf_ag(i,j,k-1)-zhalf_ag(i,j,k))* & (diff_tm(i,j,k)-kcrit) / & (diff_tm(i,j,k)-diff_tm(i,j,k-1)) else h(i,j) = 0. end if enddo enddo !----------------------------------------------------------------------- !===================================================================== end SUBROUTINE K_PBL_DEPTH !####################################################################### SUBROUTINE TKE_SURF ( is, js, u_star, kbot ) !======================================================================= !--------------------------------------------------------------------- ! Arguments (Intent in) ! u_star - surface friction velocity (m/s) ! kbot - OPTIONAL;lowest model level index (integer); ! at levels > Kbot, Mask = 0. !--------------------------------------------------------------------- integer, intent(in) :: is, js real, intent(in), dimension(:,:) :: u_star integer, intent(in), OPTIONAL, dimension(:,:) :: kbot !--------------------------------------------------------------------- ! (Intent local) !--------------------------------------------------------------------- real, dimension(SIZE(u_star,1),SIZE(u_star,2)) :: x1 integer :: ix, jx, kxp, i, j, k !======================================================================= ix = SIZE( u_star, 1 ) jx = SIZE( u_star, 2 ) kxp = SIZE( TKE, 3 ) !--------------------------------------------------------------------- x1 = bcq * u_star * u_star if( PRESENT( kbot ) ) then do j = 1,jx do i = 1,ix k = kbot(i,j) + 1 TKE(is+i-1,js+j-1,k) = x1(i,j) end do end do else do j = 1,jx do i = 1,ix TKE(is+i-1,js+j-1,kxp) = x1(i,j) end do end do endif !======================================================================= end SUBROUTINE TKE_SURF !####################################################################### end MODULE MY25_TURB_MOD