!############################# Change Log ################################## ! 5.0.0 ! !########################################################################### ! Copyright (C) 1990, 1995, 1999, 2000, 2003 - All Rights Reserved ! Regional Atmospheric Modeling System - RAMS !########################################################################### subroutine diffuse() ! +-----------------------------------------------------------------+ ! \ this routine is the subdriver to compute tendencies due to \ ! \ subgrid-scale turbulence. \ ! +-----------------------------------------------------------------+ use mem_tend, only: & tend ! %tket, %epst, %ut, %vt, %wt use mem_basic, only: & basic_g ! %up(IN), %vp(IN) use var_tables, only: & num_scalar, & ! INTENT(IN) scalar_tab ! %var_p, %var_t use mem_turb, only: & idiffk, & !INTENT(IN) turb_g, & ! %tkep, %hkm, %vkh xkhkm use mem_grid, only: & if_adap, & ! INTENT(IN) jdim, & ! INTENT(IN) ngrid, & ! INTENT(IN) grid_g, & ! %rtgt INTENT(IN) dzm, & ! INTENT(IN) dzt, & ! INTENT(IN) npatch, & ! INTENT(IN) nstbot, & ! INTENT(IN) nscl, & ! INTENT(IN) naddsc, & ! INTENT(IN) dtlt, & ! INTENT(IN) nnzp use mem_leaf, only: & leaf_g !INTENT(IN) use mem_micro, only: & micro_g !%rcp use mem_scratch, only: & scratch, & ! %vt3da, %vt3db, %vt3dc, %vt3dd, %vt3de, ! %vt3df, %vt3dg, %vt3dh, %vt3di, %vt3dj, ! %vt3dn, %scr2 vctr34 ! use node_mod, only: & mxp, & !INTENT(IN) myp, & !INTENT(IN) mzp, & !INTENT(IN) ia, & !INTENT(IN) iz, & !INTENT(IN) ja, & !INTENT(IN) jz, & !INTENT(IN) ia_1, & !INTENT(IN) ja_1, & !INTENT(IN) iz1, & !INTENT(IN) jz1, & !INTENT(IN) ibcon, & !INTENT(IN) mynum, & !INTENT(IN) nodei0, & !INTENT(IN) nodej0, & !INTENT(IN) nodemyp, & !INTENT(IN) nodemxp, & !INTENT(IN) izu, & !INTENT(IN) jzv, & !INTENT(IN) ia1, & !INTENT(IN) ja1, & !INTENT(IN) iz_1, & !INTENT(IN) jz_1, & !INTENT(IN) mynum use ke_coms, only: & alf_eps, & alf_tke use micphys, only: & level use mem_turb_scalar, only: & turb_s use mem_grell, only: & cuforc_sh_g use mem_cuparm, only : & nnqparm use mem_chem1, only: & chemistry, & nspecies_transported !ml/srf- for stilt/new turb scheme use mem_stilt, only: & stilt_g, & imassflx use ccatt_start, only: & ccatt ! intent(in) implicit none include "i8.h" !local variables: integer(kind=i8) :: mxyzp, ind integer :: n real :: s1,s2,s3 real, pointer :: scalarp(:), scalart(:), vkh_p(:), hkh_p(:) integer :: i,j,k,ksf ! srf - Large Scale Forcing for shallow and deep cumulus !real, pointer :: lsfcupar_p(:,:,:) real, target :: scr3(mxp*myp*mzp) real, target :: scr2(mxp*myp*mzp) !interface ! subroutine PBLforcing(ngrid, m1, m2, m3, ia, iz, ja, jz, & ! vt3df, scp, lsfcupar, nsc) ! ! Arguments: ! integer, INTENT(IN) :: ngrid, m1, m2, m3, ia, iz, ja, jz ! real, INTENT(IN) :: vt3df(m1,m2,m3), scp(m1,m2,m3) ! real, pointer :: lsfcupar(:,:,:) ! integer, intent(IN) :: nsc ! end subroutine PBLforcing !end interface ! CATT ! Nullifing pointer to Large Scale Forcing for GRELL CUPAR - Not used !nullify(lsfcupar_p) ! mxyzp = mxp*myp*mzp scr2 = 0.0 scr3 = 0.0 scratch%vt3dg = 0. !LFR - Avoiding overflow in strain !-srf - opt turb_k com nakanishi !-srf passando rv e rtp em vez de vt3dp e vt3dq !if(level == 0) then ! scratch%vt3dp= 0. ! scratch%vt3dq= 0. !elseif (level > 0) then ! call atob(mxyzp,basic_g(ngrid)%rv,scratch%vt3dp) ! call atob(mxyzp,basic_g(ngrid)%rtp,scratch%vt3dq) !end if if (if_adap==0) then call strain(mzp,mxp,myp,ia,iz,ja,jz & ,ia_1,ja_1,iz1,jz1,jdim & ,basic_g(ngrid)%up (:,:,:) ,basic_g(ngrid)%vp (:,:,:) & ,basic_g(ngrid)%wp (:,:,:) ,scratch%vt3da (:) & ,scratch%vt3db (:) ,scratch%vt3dc (:) & ,scratch%vt3dd (:) ,scratch%vt3de (:) & ,scratch%vt3df (:) ,scratch%vt3dg (:) & ,scratch%vt3dh (:) ,scratch%vt3di (:) & ,scratch%vt3dn (:) ,scr2 (:) & ,idiffk(ngrid)) else call strain_adap(mzp,mxp,myp,ia,iz,ja,jz & ,ia_1,ja_1,iz1,jz1,jdim & ,grid_g(ngrid)%lpu (:,:) ,grid_g(ngrid)%lpv (:,:) & ,grid_g(ngrid)%lpw (:,:) ,basic_g(ngrid)%up (:,:,:) & ,basic_g(ngrid)%vp (:,:,:) ,basic_g(ngrid)%wp (:,:,:) & ,scratch%vt3da (:) ,scratch%vt3db (:) & ,scratch%vt3dc (:) ,scratch%vt3dd (:) & ,scratch%vt3de (:) ,scratch%vt3df (:) & ,scratch%vt3dg (:) ,scratch%vt3dh (:) & ,scratch%vt3di (:) ,scratch%vt3dn (:) & ,scr2 (:) ,idiffk(ngrid) & ,grid_g(ngrid)%dxm (:,:) ,grid_g(ngrid)%dxt (:,:) & ,grid_g(ngrid)%dxu (:,:) ,grid_g(ngrid)%dxv (:,:) & ,grid_g(ngrid)%dym (:,:) ,grid_g(ngrid)%dyt (:,:) & ,grid_g(ngrid)%dyu (:,:) ,grid_g(ngrid)%dyv (:,:) & ,dzm,dzt) endif if (level<=1) & scratch%vt3dp = 0. if (level>=2) & call ae1_l(int(mxyzp,i8), scratch%vt3dp(:), micro_g(ngrid)%rcp(:,:,:)) !-srf 29/12/2008 adapted from OLAM call bruvais_OLAM(mzp, mxp, myp, ia, iz, ja, jz, & basic_g(ngrid)%theta(1,1,1), basic_g(ngrid)%rtp(1,1,1), & basic_g(ngrid)%rv(1,1,1), scratch%vt3dp(1), & basic_g(ngrid)%pp(1,1,1), basic_g(ngrid)%pi0(1,1,1), & scratch%vt3dj(1), grid_g(ngrid)%rtgt(1,1), & grid_g(ngrid)%lpw(1,1)) !ml/srf- for new turn scheme if (idiffk(ngrid) <= 3 .or. idiffk(ngrid) == 7.or. idiffk(ngrid) == 8) then call mxdefm(mzp,mxp,myp,ia,iz,ja,jz,ibcon,jdim & ,scratch%vt3dh (:) ,scratch%vt3di (:) & ,scratch%vt3dj (:) ,scratch%vt3dk (:) & ,scratch%scr1 (:) ,scr2 (:) & ,basic_g(ngrid)%dn0 (:,:,:) ,grid_g(ngrid)%rtgt (:,:) & ,grid_g(ngrid)%dxt (:,:) ,grid_g(ngrid)%dyt (:,:) & ,grid_g(ngrid)%lpw (:,:) ,mynum ) if (CCATT==1 .and. chemistry >= 0) then !srf------ !coef de difusao horizontal diferente para tracers call mxdefm_tracer(mzp,mxp,myp,ia,iz,ja,jz & ,ibcon,jdim,scratch%vt3dh(:),scr3(:) & ,basic_g(ngrid)%dn0(:,:,:),grid_g(ngrid)%dxt(:,:),& grid_g(ngrid)%dyt(:,:),grid_g(ngrid)%lpw(:,:),mynum) !srf------ endif endif !ML -> Nananishi and Niino (2004) scheme based on Mellor-Yamada Level 2.5 if (idiffk(ngrid) == 7) then if(IMASSFLX==1)then call nakanishi(mzp, mxp, myp, npatch, ia, iz, ja, jz, jdim & ,turb_g(ngrid)%tkep ,tend%tket ,scratch%vt3dd & ,scratch%vt3de ,scratch%vt3dh ,scratch%vt3di & ,scratch%vt3dj ,scratch%scr1 ,grid_g(ngrid)%rtgt & !srf ,basic_g(ngrid)%theta ,scratch%vt3dp ,scratch%vt3dq & ,basic_g(ngrid)%theta ,basic_g(ngrid)%rv ,basic_g(ngrid)%rtp & ! ,basic_g(ngrid)%dn0 ,basic_g(ngrid)%up ,basic_g(ngrid)%vp & ,leaf_g(ngrid)%veg_rough ,leaf_g(ngrid)%patch_rough,leaf_g(ngrid)%tstar & ,leaf_g(ngrid)%ustar ,leaf_g(ngrid)%patch_area ,turb_g(ngrid)%sflux_u & ,turb_g(ngrid)%sflux_v ,turb_g(ngrid)%sflux_t ,grid_g(ngrid)%lpu & ,grid_g(ngrid)%lpv ,grid_g(ngrid)%lpw ,turb_g(ngrid)%kpbl & ,stilt_g(ngrid)%pblhgt ,stilt_g(ngrid)%lmo ,stilt_g(ngrid)%ltscale & ,stilt_g(ngrid)%sigw ) !- from Marcos Longo !---------------------------------------------------------------------------------------! ! Integrating average turbulence parameters for mass flux. ! !---------------------------------------------------------------------------------------! call prepare_timeavg_driver(mzp,mxp,myp,ia,iz,ja,jz,dtlt,ngrid,idiffk(ngrid)) !---------------------------------------------------------------------------------------! else call nakanishi_light(mzp, mxp, myp, npatch, ia, iz, ja, jz, jdim & ,turb_g(ngrid)%tkep ,tend%tket ,scratch%vt3dd & ,scratch%vt3de ,scratch%vt3dh ,scratch%vt3di & ,scratch%vt3dj ,scratch%scr1 ,grid_g(ngrid)%rtgt & !srf-opt ,basic_g(ngrid)%theta ,scratch%vt3dp ,scratch%vt3dq & ,basic_g(ngrid)%theta ,basic_g(ngrid)%rv ,basic_g(ngrid)%rtp & ! ,basic_g(ngrid)%dn0 ,basic_g(ngrid)%up ,basic_g(ngrid)%vp & ,leaf_g(ngrid)%veg_rough ,leaf_g(ngrid)%patch_rough,leaf_g(ngrid)%tstar & ,leaf_g(ngrid)%ustar ,leaf_g(ngrid)%patch_area ,turb_g(ngrid)%sflux_u & ,turb_g(ngrid)%sflux_v ,turb_g(ngrid)%sflux_t ,grid_g(ngrid)%lpu & ,grid_g(ngrid)%lpv ,grid_g(ngrid)%lpw & !,turb_g(ngrid)%kpbl & !srf !,stilt_g(ngrid)%pblhgt & !srf ,stilt_g(ngrid)%lmo & !,stilt_g(ngrid)%ltscale & !srf !,stilt_g(ngrid)%sigw & !srf ) endif endif !ML if (idiffk(ngrid)==1) then call tkemy(mzp,mxp,myp,ia,iz,ja,jz,ibcon,jdim,nodei0(mynum,ngrid),nodej0(mynum,ngrid) & ,turb_g(ngrid)%tkep (:,:,:) ,tend%tket (:) & ,scratch%vt3dh (:) ,scratch%vt3di (:) & ,scratch%vt3dj (:) ,scratch%scr1 (:) & ,grid_g(ngrid)%rtgt (:,:) ,basic_g(ngrid)%theta (:,:,:) & ,basic_g(ngrid)%dn0 (:,:,:) ,basic_g(ngrid)%up (:,:,:) & ,basic_g(ngrid)%vp (:,:,:) ,basic_g(ngrid)%wp (:,:,:) & ,turb_g(ngrid)%sflux_u(:,:) ,turb_g(ngrid)%sflux_v(:,:) & ,turb_g(ngrid)%sflux_w(:,:) ,turb_g(ngrid)%sflux_t(:,:),vctr34 & ,grid_g(ngrid)%lpw (:,:) ,grid_g(ngrid)%lpu (:,:) & ,grid_g(ngrid)%lpv (:,:)) endif if (idiffk(ngrid)==4) then call mxtked(mzp,mxp,myp,ia,iz,ja,jz & ,ibcon,jdim & ,turb_g(ngrid)%tkep (:,:,:) ,tend%tket (:) & ,basic_g(ngrid)%up (:,:,:) ,basic_g(ngrid)%vp (:,:,:) & ,basic_g(ngrid)%wp (:,:,:) ,basic_g(ngrid)%rtp (:,:,:) & ,basic_g(ngrid)%rv (:,:,:) ,basic_g(ngrid)%theta (:,:,:) & ,scratch%vt3da (:) ,scratch%vt3dc (:) & ,scratch%vt3dh (:) ,scratch%vt3dj (:) & ,scratch%scr1 (:) ,scr2 (:) & ,turb_g(ngrid)%sflux_u(:,:) ,turb_g(ngrid)%sflux_v(:,:) & ,turb_g(ngrid)%sflux_w(:,:) ,turb_g(ngrid)%sflux_t(:,:) & ,grid_g(ngrid)%dxt (:,:) ,grid_g(ngrid)%rtgt (:,:) & ,grid_g(ngrid)%lpw (:,:) ) endif !_STC............................................................ !_STC Call to subroutine tkescl for E-l closure !_STC (S. Trini Castelli) !_STC............................................................ if (idiffk(ngrid)==5) then call tkescl(mzp,mxp,myp,npatch,ia,iz,ja,jz & ,turb_g(ngrid)%tkep(:,:,:),tend%tket(:) & ,turb_g(ngrid)%epsp(:,:,:),tend%epst(:) & ,scratch%vt3da(:),scratch%vt3dc(:) & ,scratch%vt3dh(:),scratch%vt3di(:) & ,scratch%vt3dj(:),scratch%scr1(:) & ,scr2(:) ,grid_g(ngrid)%rtgt(:,:) & ,scratch%vt3dd(:),scratch%vt3de(:),grid_g(ngrid)%dxt(:,:) & ,leaf_g(ngrid)%ustar(:,:,:),leaf_g(ngrid)%patch_area(:,:,:) & ,grid_g(ngrid)%lpw(:,:),basic_g(ngrid)%dn0(:,:,:) ) endif !_STC............................................................ !_STC Call to subroutine tkeeps for E-eps closure !_STC (S. Trini Castelli) !_STC............................................................ if (idiffk(ngrid)==6) then call tkeeps(mzp,mxp,myp,npatch,ia,iz,ja,jz & ,turb_g(ngrid)%tkep(:,:,:),tend%tket(:) & ,turb_g(ngrid)%epsp(:,:,:),tend%epst(:) & ,scratch%vt3da(:),scratch%vt3dc(:) & ,scratch%vt3dh(:),scratch%vt3di(:) & ,scratch%vt3dj(:),scratch%scr1(:) & ,scr2(:) ,grid_g(ngrid)%rtgt(:,:) & ,leaf_g(ngrid)%ustar(:,:,:),leaf_g(ngrid)%patch_area(:,:,:) & ,grid_g(ngrid)%lpw(:,:),basic_g(ngrid)%dn0(:,:,:) ) endif !_STC.................................................. !_STC Note: from subroutines TKESCL, TKEEPS : !_STC VT3DI=Ke !_STC SCR1=Km !_STC VT3DH = Kh !_STC SCR2 = SCR1 = Km !_STC.................................................. !_STC............................................................ if (idiffk(ngrid) == 8) then !----------------------------------- !Call subrotines HCV_driver (Campos Velho - KZZ) !----------------------------------- call HCV_driver(mzp,mxp,myp,ia,iz,ja,jz,mynum & ,basic_g(ngrid)%up (:,:,:) & ,basic_g(ngrid)%vp (:,:,:) & ,basic_g(ngrid)%theta (:,:,:) & ,basic_g(ngrid)%rv (:,:,:) & ,basic_g(ngrid)%dn0 (:,:,:) & ,turb_g(ngrid)%sflux_t(:,:) & ,turb_g(ngrid)%sflux_r(:,:) & ,turb_g(ngrid)%sflux_u(:,:) & ,turb_g(ngrid)%sflux_v(:,:) & ,grid_g(ngrid)%rtgt (:,:) & ,grid_g(ngrid)%lpw (:,:) & ,scratch%scr1 (:) & ,scratch%vt3dh (:) & ,scratch%vt3di (:) & ,scratch%vt3dj (:) & ,scratch%vt3dk (:) ) endif call klbnd(mzp,mxp,myp,ibcon,jdim & ,scratch%scr1 (:),basic_g(ngrid)%dn0(:,:,:),grid_g(ngrid)%lpw(:,:)) call klbnd(mzp,mxp,myp,ibcon,jdim & ,scr2 (:),basic_g(ngrid)%dn0(:,:,:),grid_g(ngrid)%lpw(:,:)) call klbnd(mzp,mxp,myp,ibcon,jdim & ,scratch%vt3dh(:),basic_g(ngrid)%dn0(:,:,:),grid_g(ngrid)%lpw(:,:)) ! CATT if (CCATT==1 .and. chemistry >= 0) then !srf---- call klbnd(mzp,mxp,myp,ibcon,jdim,scr3(:) & ,basic_g(ngrid)%dn0(:,:,:),grid_g(ngrid)%lpw(:,:)) endif !_STC ....... boundary conditions even on Ke diffusion coefficient if (idiffk(ngrid)==5 .or. idiffk(ngrid)==6) & call klbnd(mzp,mxp,myp,ibcon,jdim & ,scratch%vt3di(:),basic_g(ngrid)%dn0(:,:,:),grid_g(ngrid)%lpw(:,:)) !bob swap new hkm, vkm, and vkh with past time level: lagged K's have !bob internal lateral boundary values from neighboring nodes ind = 0 do j = 1,nodemyp(mynum,ngrid) do i = 1,nodemxp(mynum,ngrid) do k = 1,nnzp(ngrid) ind = ind + 1 s1 = scr2(ind) s2 = scratch%scr1(ind) s3 = scratch%vt3dh(ind) scr2(ind) = turb_g(ngrid)%hkm(k,i,j) scratch%scr1(ind) = turb_g(ngrid)%vkm(k,i,j) scratch%vt3dh(ind) = turb_g(ngrid)%vkh(k,i,j) !! also for vt3di = K(tke) ????? 22 March 02 turb_g(ngrid)%hkm(k,i,j) = s1 turb_g(ngrid)%vkm(k,i,j) = s2 turb_g(ngrid)%vkh(k,i,j) = s3 enddo enddo enddo if (if_adap==0) then !call dumpVarAllLatLonk(basic_g(ngrid)%up, 'tUP' ,415,0,0,1,mxp,1,myp,1,mzp,0.0,0.0) !ok !call dumpVarAllLatLonk(basic_g(ngrid)%Vp, 'tVP' ,415,0,0,1,mxp,1,myp,1,mzp,0.0,0.0) !ok !call dumpVarAllLatLonk(basic_g(ngrid)%Wp, 'tWP' ,415,0,0,1,mxp,1,myp,1,mzp,0.0,0.0) !ok !call dumpVarAllLatLonk(basic_g(ngrid)%dn0, 'DN0' ,415,0,0,1,mxp,1,myp,1,mzp,0.0,0.0) !ok !call dumpVarAllLatLonk(basic_g(ngrid)%dn0v, 'DN0V' ,415,0,0,1,mxp,1,myp,1,mzp,0.0,0.0) !ok !call dumpVarAllLatLonk(basic_g(ngrid)%dn0u, 'DN0U' ,415,0,0,1,mxp,1,myp,1,mzp,0.0,0.0) !ok !call dumpVarAllLatLonk(grid_g(ngrid)%rtgu ,'rtgu',415,0,0,1,mxp,1,myp,1,1,0.0,0.0) !ok !call dumpVarAllLatLonk(grid_g(ngrid)%rtgv ,'rtgv',415,0,0,1,mxp,1,myp,1,1,0.0,0.0) !ok !call dumpVarAllLatLonk(grid_g(ngrid)%rtgt ,'rtgt',415,0,0,1,mxp,1,myp,1,1,0.0,0.0) !ok !call dumpVarAllLatLonk(turb_g(ngrid)%sflux_u,'sflux_u',415,0,0,1,mxp,1,myp,1,1,0.0,0.0) !ok !call dumpVarAllLatLonk(turb_g(ngrid)%sflux_v,'sflux_v',415,0,0,1,mxp,1,myp,1,1,0.0,0.0) !ok !call dumpVarAllLatLonk(turb_g(ngrid)%sflux_w,'sflux_w',415,0,0,1,mxp,1,myp,1,1,0.0,0.0) !***BAD*** call diffvel(mzp,mxp,myp,ia,iz,ja,jz,jdim,ia_1,ja_1 & ,ia1,ja1,iz_1,jz_1,iz1,jz1,izu,jzv,idiffk(ngrid) & ,basic_g(ngrid)%up (:,:,:) ,basic_g(ngrid)%vp (:,:,:) & ,basic_g(ngrid)%wp (:,:,:) ,tend%ut (:) & ,tend%vt (:) ,tend%wt (:) & ,scratch%vt3da (:) ,scratch%vt3db (:) & ,scratch%vt3dc (:) ,scratch%vt3dd (:) & ,scratch%vt3de (:) ,scratch%vt3df (:) & ,scratch%vt3dg (:) ,scratch%vt3dj (:) & ,scratch%vt3dk (:) ,scratch%vt3dl (:) & ,scratch%vt3dm (:) ,scratch%vt3dn (:) & ,scratch%vt3do (:) ,grid_g(ngrid)%rtgu (:,:) & ,grid_g(ngrid)%rtgv (:,:) ,grid_g(ngrid)%rtgt (:,:) & ,turb_g(ngrid)%sflux_u(:,:) ,turb_g(ngrid)%sflux_v(:,:) & ,turb_g(ngrid)%sflux_w(:,:) ,basic_g(ngrid)%dn0 (:,:,:) & ,basic_g(ngrid)%dn0u (:,:,:) ,basic_g(ngrid)%dn0v (:,:,:) & ,scratch%scr1 (:) ,scr2 (:),ibcon,mynum) !call dumpVarAllLatLonk(tend%wt, 'tWT' ,433,0,0,1,mxp,1,myp,1,mzp,0.0,0.0) else call diffvel_adap(mzp,mxp,myp,ia,iz,ja,jz,jdim & ,iz1,jz1,izu,jzv,idiffk(ngrid) & ,basic_g(ngrid)%up (:,:,:) ,basic_g(ngrid)%vp (:,:,:) & ,basic_g(ngrid)%wp (:,:,:) ,tend%ut (:) & ,tend%vt (:) ,tend%wt (:) & ,scratch%vt3da (:) ,scratch%vt3db (:) & ,scratch%vt3dc (:) ,scratch%vt3dd (:) & ,scratch%vt3de (:) ,scratch%vt3df (:) & ,scratch%vt3dg (:) ,scratch%vt3dj (:) & ,scratch%vt3dk (:) ,scratch%vt3dl (:) & ,scratch%vt3dm (:) ,scratch%vt3dn (:) & ,scratch%vt3do (:) ,grid_g(ngrid)%aru (:,:,:) & ,grid_g(ngrid)%arv (:,:,:) ,grid_g(ngrid)%arw (:,:,:) & ,grid_g(ngrid)%volu (:,:,:) ,grid_g(ngrid)%volv (:,:,:) & ,grid_g(ngrid)%volw (:,:,:) ,grid_g(ngrid)%lpu (:,:) & ,grid_g(ngrid)%lpv (:,:) ,grid_g(ngrid)%lpw (:,:) & ,turb_g(ngrid)%sflux_u(:,:) ,turb_g(ngrid)%sflux_v(:,:) & ,turb_g(ngrid)%sflux_w(:,:) ,basic_g(ngrid)%dn0 (:,:,:) & ,basic_g(ngrid)%dn0u (:,:,:) ,basic_g(ngrid)%dn0v (:,:,:) & ,scratch%scr1 (:) ,scr2 (:) & ,grid_g(ngrid)%topma (:,:) ,ibcon,mynum) endif ! Convert momentum K's to scalar K's, if necessary !-ml/srf - for new turb scheme if (idiffk(ngrid) <= 3 .or. idiffk(ngrid) == 7 .or. idiffk(ngrid) == 8) then do ind = 1,mxyzp scr2(ind) = scr2(ind) * xkhkm(ngrid) enddo elseif (idiffk(ngrid) == 4) then do ind = 1,mxyzp scratch%vt3di(ind) = 2. * scratch%scr1(ind) enddo endif !- CATT if (CCATT==1.and. CHEMISTRY >=0) then ind = 0 do j = 1,nodemyp(mynum,ngrid) do i = 1,nodemxp(mynum,ngrid) do k = 1,nnzp(ngrid) ind = ind + 1 s1 = scr3(ind) scr3(ind) = turb_s(ngrid)%hksc(k,i,j) turb_s(ngrid)%hksc(k,i,j) = s1 ! salva o atual coef para o proximo passo no tempo. enddo enddo enddo if (idiffk(ngrid)<=3 .or. idiffk(ngrid) == 7.or. idiffk(ngrid) == 8) then ind = 0 do j = 1,nodemyp(mynum,ngrid) do i = 1,nodemxp(mynum,ngrid) do k = 1,nnzp(ngrid) ind = ind + 1 scr3(ind) = scr3(ind) * xkhkm(ngrid) enddo enddo enddo endif endif do n = 1,num_scalar(ngrid) scalarp => scalar_tab(n,ngrid)%a_var_p scalart => scalar_tab(n,ngrid)%a_var_t scratch%vt2da = 0. if (nstbot==1) then if (scalar_tab(n,ngrid)%name=='THP' .or. & scalar_tab(n,ngrid)%name=='THC') then call atob(mxp*myp, turb_g(ngrid)%sflux_t(:,:), scratch%vt2da(:)) ! Large Scale Forcing for GRELL CUPAR !if (nnqparm(ngrid)>=2) then ! !--------------------------------------- ! !srf- Large Scale Forcing for GRELL CUPAR ! lsfcupar_p => cuforc_sh_g(ngrid)%lsfth !--------------------------------------- ! endif elseif (scalar_tab(n,ngrid)%name=='RTP') then call atob(mxp*myp, turb_g(ngrid)%sflux_r(:,:), scratch%vt2da(:)) ! Large Scale Forcing for GRELL CUPAR not used ! CATT !if (nnqparm(ngrid)>=2) then ! !--------------------------------------- ! !srf- Large Scale Forcing for GRELL CUPAR ! lsfcupar_p => cuforc_sh_g(ngrid)%lsfrt ! !--------------------------------------- !endif endif endif ! 3/10/01 - Define ksf below, the "K scalar flag", to let subroutine diffsclr ! know which vertical K is being passed to it. If diffsclr sees that it's ! a different K from the previous one, diffsclr will re-compute the tridiff ! matrix coefficients. In order to use vertical scalar K's other than ! vt3dh and vt3di, use ksf = 3, ksf = 4, etc. for each different K. !_STC.................................................. !_STC Corrections to account for the new idiffk options !_STC for E-l and E-eps closure. Isotropy hypothesis. !_STC (S. Trini Castelli) !_STC.................................................. if (scalar_tab(n,ngrid)%name=='TKEP') then vkh_p => scratch%vt3di hkh_p => scr2 !-ml/srf - for new turb scheme if (idiffk(ngrid) >= 4 .and. idiffk(ngrid) /= 7 .and. idiffk(ngrid) /= 8) & hkh_p => scratch%vt3di ! if (idiffk(ngrid)>=4) hkh_p => scratch%vt3di ksf = 1 elseif (scalar_tab(n,ngrid)%name=='EPSP') then vkh_p => scratch%vt3di hkh_p => scr2 !-ml/srf - for new turb scheme if (idiffk(ngrid) >= 4 .and. idiffk(ngrid) /= 7.and. idiffk(ngrid) /= 8) & hkh_p => scratch%vt3di ! if (idiffk(ngrid)>=4) hkh_p => scratch%vt3di ksf = 3 ! Convert Ktke to Keps; it will be converted back after use below call ae1t0_l(mxyzp, vkh_p, vkh_p, (ALF_EPS/ALF_TKE)) call ae1t0_l(mxyzp, hkh_p, hkh_p, (ALF_EPS/ALF_TKE)) else vkh_p => scratch%vt3dh hkh_p => scr2 !-ml/srf - for new turb scheme if (idiffk(ngrid)>=4 .and. idiffk(ngrid) /= 7.and. idiffk(ngrid) /= 8) & hkh_p => scratch%vt3dh ! if (idiffk(ngrid)>=4) hkh_p => scratch%vt3di ksf = 2 endif ! CATT if (CCATT==1 .and. CHEMISTRY >=0) then !srf----------------- Hor. Diffusion Coef for tracers if(n > (num_scalar(ngrid) - (NADDSC + NSPECIES_TRANSPORTED))) then if (idiffk(ngrid) < 4 .or. idiffk(ngrid) == 7.or. idiffk(ngrid) == 8) then hkh_p => scr3 endif endif !srf-------------------- endif if (if_adap==0) then ! NEW WAY call diffsclr(mzp, mxp, myp, ia, iz, ja, jz, jdim, & ia_1, ja_1, ia1, ja1, iz_1, jz_1, iz1, jz1, n, ksf, & scalarp(:), scalart(:), & scratch%vt3da(:), scratch%vt3db(:), scratch%vt3df(:), & scratch%vt3dg(:), scratch%vt3dj(:), scratch%vt3dk(:), & scratch%vt3do(:), scratch%vt3dc(:), scratch%vt3dd(:), & scratch%vt3dl(:), scratch%vt3dm(:), scratch%vt2db(:), & grid_g(ngrid)%rtgt(:,:), scratch%vt2da(:), & basic_g(ngrid)%dn0(:,:,:), & vkh_p(:) , hkh_p(:) ) ! IF (nnqparm(ngrid)>=2) then ! SGScale Forcing for GRELL CUPAR if (scalar_tab(n,ngrid)%name=='THP' .or. scalar_tab(n,ngrid)%name=='THC') & call PBLforcing(ngrid, mzp, mxp, myp, ia, iz, ja, jz, & scratch%vt3df, scalarp(:), cuforc_sh_g(ngrid)%lsfth, n) if (scalar_tab(n,ngrid)%name=='RTP') & call PBLforcing(ngrid, mzp, mxp, myp, ia, iz, ja, jz, & scratch%vt3df, scalarp(:), cuforc_sh_g(ngrid)%lsfrt, n) ENDIF else call diffsclr_adap(mzp,mxp,myp,ia,iz,ja,jz,jdim,n,ksf & ,grid_g(ngrid)%lpw(:,:) ,scalarp & ,scalart ,scratch%vt3da (:) & ,scratch%vt3dc (:) ,scratch%vt3df (:) & ,scratch%vt3dg (:) ,scratch%vt3dj (:) & ,scratch%vt3dk (:) ,scratch%vt3dl (:) & ,scratch%vt3dm (:) ,scratch%vt3do (:) & ,scratch%vt2da (:) ,scratch%vt2db (:) & ,basic_g(ngrid)%dn0 (:,:,:) ,vkh_p & ,hkh_p ,grid_g(ngrid)%aru (:,:,:) & ,grid_g(ngrid)%arv (:,:,:) ,grid_g(ngrid)%arw (:,:,:) & ,grid_g(ngrid)%volt (:,:,:) ,scratch%vt3db (:) & ,grid_g(ngrid)%dxu (:,:) ,grid_g(ngrid)%dyv (:,:) & ,grid_g(ngrid)%topma(:,:) ) endif if (scalar_tab(n,ngrid)%name == 'EPSP') then call ae1t0_l(mxyzp, vkh_p, vkh_p, (ALF_TKE/ALF_EPS)) call ae1t0_l(mxyzp, hkh_p, hkh_p, (ALF_TKE/ALF_EPS)) endif enddo end subroutine diffuse ! ****************************************************************** subroutine strain(m1,m2,m3,ia,iz,ja,jz,ia_1,ja_1,iz1,jz1 & ,jd,up,vp,wp,vt3da,vt3db,vt3dc,vt3dd,vt3de & ,vt3df,vt3dg,vt3dh,vt3di,vt3dn,scr2,idiffk) implicit none integer, intent(in) :: m1 & , m2 & , m3 & , ia & , iz & , ja & , jz & , ia_1 & , ja_1 & , iz1 & , jz1 & , jd & , idiffk real, dimension(m1,m2,m3), intent(in) :: up,vp,wp real, dimension(m1,m2,m3), intent(inout):: vt3da, & vt3db, & vt3dc, & vt3dd, & vt3de, & vt3df, & vt3dg, & vt3dh, & vt3di, & vt3dn, & scr2 !local variables: integer :: i,j,k call grad(m1, m2, m3, ia , iz1, ja , jz , up, vt3da, 'XDIR', 'UPNT') call grad(m1, m2, m3, ia_1, iz , ja_1, jz , vp, vt3db, 'XDIR', 'VPNT') call grad(m1, m2, m3, ia_1, iz , ja , jz , wp, vt3df, 'XDIR', 'WPNT') call grad(m1, m2, m3, ia_1, iz , ja_1, jz , up, vt3dn, 'YDIR', 'UPNT') call grad(m1, m2, m3, ia , iz , ja , jz1, vp, vt3dc, 'YDIR', 'VPNT') call grad(m1, m2, m3, ia , iz , ja_1, jz , wp, vt3dg, 'YDIR', 'WPNT') call grad(m1, m2, m3, ia_1, iz , ja , jz , up, vt3dd, 'ZDIR', 'UPNT') call grad(m1, m2, m3, ia , iz , ja_1, jz , vp, vt3de, 'ZDIR', 'VPNT') if(idiffk.ge.3 .and. idiffk /= 7)then call grad(m1,m2,m3,ia,iz,ja,jz,wp,scr2,'ZDIR','WPNT') endif if (idiffk .le. 2 .or. idiffk == 7 .or. idiffk == 8) then do j = ja,jz do i = ia,iz do k = 2,m1-1 vt3dh(k,i,j) =2. * (vt3da(k,i,j) * vt3da(k,i,j) & + vt3dc(k,i,j) * vt3dc(k,i,j)) & + .0625 * (vt3db(k,i,j) + vt3db(k,i-1,j) & + vt3db(k,i,j-jd) + vt3db(k,i-1,j-jd) & + vt3dn(k,i,j) + vt3dn(k,i-1,j) & + vt3dn(k,i,j-jd) + vt3dn(k,i-1,j-jd)) ** 2 vt3di(k,i,j) = .0625 * ((vt3dd(k,i,j) + vt3dd(k-1,i,j) & + vt3dd(k,i-1,j) + vt3dd(k-1,i-1,j)) ** 2 & + (vt3de(k,i,j) + vt3de(k-1,i,j) & + vt3de(k,i,j-jd) + vt3de(k-1,i,j-jd)) ** 2) enddo enddo enddo else do j = ja,jz do i = ia,iz do k = 2,m1-1 vt3da(k,i,j) = 2. * vt3da(k,i,j) vt3dc(k,i,j) = 2. * vt3dc(k,i,j) scr2(k,i,j) = 2. * scr2(k,i,j) vt3db(k,i,j) = vt3db(k,i,j) + vt3dn(k,i,j) vt3dn(k,i,j) = vt3db(k,i,j) vt3dd(k,i,j) = vt3dd(k,i,j) + vt3df(k,i,j) vt3de(k,i,j) = vt3de(k,i,j) + vt3dg(k,i,j) vt3di(k,i,j) = 0.333333 & * (vt3da(k,i,j) + vt3dc(k,i,j) + scr2(k,i,j)) enddo enddo do k = 2,m1-1 vt3da(k,iz1,j) = 2. * vt3da(k,iz1,j) vt3db(k,ia_1,j) = vt3db(k,ia_1,j) + vt3dn(k,ia_1,j) vt3dn(k,ia_1,j) = vt3db(k,ia_1,j) vt3dd(k,ia_1,j) = vt3dd(k,ia_1,j) + vt3df(k,ia_1,j) enddo enddo do i = ia_1,iz do k = 2,m1-1 vt3dc(k,i,jz1) = 2. * vt3dc(k,i,jz1) vt3db(k,i,ja_1) = vt3db(k,i,ja_1) + vt3dn(k,i,ja_1) vt3dn(k,i,ja_1) = vt3db(k,i,ja_1) vt3de(k,i,ja_1) = vt3de(k,i,ja_1) + vt3dg(k,i,ja_1) enddo enddo do j = ja,jz do i = ia,iz do k = 2,m1-1 vt3dh(k,i,j) = .5 * ( & (vt3da(k,i,j) - vt3di(k,i,j)) ** 2 & + (vt3dc(k,i,j) - vt3di(k,i,j)) ** 2 & + ( scr2(k,i,j) - vt3di(k,i,j)) ** 2) & + .0625 * ((vt3db(k,i,j) + vt3db(k,i-1,j) & + vt3db(k,i,j-jd) + vt3db(k,i-1,j-jd)) ** 2 & + (vt3dd(k,i,j) + vt3dd(k,i-1,j) & + vt3dd(k-1,i,j) + vt3dd(k-1,i-1,j)) ** 2 & + (vt3de(k,i,j) + vt3de(k-1,i,j) & + vt3de(k,i,j-jd) + vt3de(k-1,i,j-jd)) ** 2) vt3di(k,i,j) = vt3dh(k,i,j) enddo enddo enddo endif return end subroutine strain ! ***************************************************************** subroutine bruvais(m1, m2, m3, ia, iz, ja, jz, theta, rtp, rv, rcp, & pp, pi0, en2, rtgt, lpw_R) use mem_scratch, only: & vctr11, & !INTENT(INOUT) vctr12, & !INTENT(INOUT) vctr32, & !INTENT(INOUT) vctr1, & !INTENT(INOUT) vctr2, & !INTENT(INOUT) vctr3, & !INTENT(INOUT) vctr4 !INTENT(INOUT) use micphys, only: & level, & !INTENT(in) ipris, & !INTENT(in) isnow, & !INTENT(in) igraup, & !INTENT(in) iaggr, & !INTENT(in) ihail !INTENT(in) use mem_grid, only: & zt, & !INTENT(in) nzp, & !INTENT(in) nz, & !INTENT(in) nzpmax !INTENT(in) use rconstants, only: & alvl, & !INTENT(in) rgas, & !INTENT(in) ep, & !INTENT(in) cp, & !INTENT(in) alvi, & !INTENT(in) p00, & !INTENT(in) cpor, & !INTENT(in) g !INTENT(in) implicit none ! Arguments: integer, INTENT(IN) :: m1, m2, m3, ia, iz, ja, jz real, INTENT(IN) :: theta(m1,m2,m3), pp(m1,m2,m3), pi0(m1,m2,m3), & rtp(m1,m2,m3), rcp(m1,m2,m3), rv(m1,m2,m3) real, INTENT(INOUT) :: en2(m1,m2,m3) real, INTENT(IN) :: rtgt(m2,m3) real, INTENT(IN) :: lpw_R(m2,m3) !local variables integer :: lpw(m2,m3) integer :: i, j, k, iweten, iwdiffk, ki, k2, k1 real :: c1, c2, c3, ci1, ci2, ci3, rvlsi, rvii real :: pi(nzpmax), temp(nzpmax), prt(nzpmax), rvls(nzpmax), rc(nzpmax) ! **(JP)** fatora expressoes logicas para fora dos lacos logical :: log1, log2, log3, log4 !lfr: Solving a problem with integer inside vtables lpw=int(lpw_R) ! calculate brunt-vaisalla frequency squared (en2) iweten = 1 iwdiffk = 0 c1 = alvl/rgas c2 = ep*alvl**2/(cp*rgas) c3 = alvl/cp ci1 = alvi/rgas ci2 = ep*alvi**2/(cp*rgas) ci3 = alvi/cp ! **(JP)** fatora expressoes logicas para fora dos lacos log1 = level>=1 log2 = level>=2 .and. iweten==1 log3 = (ipris>=1 .or. isnow>=1 .or. igraup>=1 .or. iaggr>=1 .or. ihail>=1) & .and. level==3 log4 = level==3 ! **(JP)** fim modificacao ! calculate potential temperature profile do j=ja,jz do i=ia,iz k2 = lpw(i,j) k1 = k2 - 1 do k=k1,m1 vctr11(k) = theta(k,i,j) vctr12(k) = theta(k,i,j) vctr32(k) = 0. enddo ! **(JP)** fatora expressoes logicas para fora dos lacos if (log1) then ! **(JP)** fim modificacao do k=k1,m1 vctr12(k) = vctr11(k)*(1. + 0.61*rv(k,i,j)) vctr32(k) = (rtp(k,i,j) - rv(k,i,j)) enddo endif ! check for saturation if level is 2 or greater. ! **(JP)** fatora expressoes logicas para fora dos lacos if (log2) then ! **(JP)** fim modificacao do k=k1,m1 pi(k) = (pp(k,i,j) + pi0(k,i,j))/cp temp(k) = theta(k,i,j)*pi(k) prt(k) = p00*pi(k)**cpor enddo call mrsl(m1,prt(:), temp(:), rvls(:)) do k=k2-1,m1 vctr2(k) = c1 vctr3(k) = c2 vctr4(k) = c3 enddo ki = m1 + 1 ! if any ice phase microphysics are activated .... ! **(JP)** fatora expressoes logicas para fora dos lacos if (log3) then ! **(JP)** fim modificacao ! find level of -20 c. assume ice saturation above this ! level. do k=k1,m1 if (temp(k)<=253.16) then ki = k go to 10 endif enddo ki = m1 + 1 10 continue call mrsi(m1-ki+1, prt(ki), temp(ki), rvls(ki)) !-srf 19/03/2005 !bug: Subscript out of range for array prt ! subscript=0, lower bound=1, upper bound=132, dimension=1 ! tempc < 253 sobre parte da Antartica. ! ki = max(ki,k2) !solucao 1 ! call mrsi(1,prt(ki-1),temp(ki-1),rvlsi) if (ki>1) call mrsi(1, prt(ki-1), temp(ki-1), rvlsi) ! solucao 2 !srf do k=ki,m1 vctr2(k) = ci1 vctr3(k) = ci2 vctr4(k) = ci3 enddo endif ! **(JP)** fatora expressoes logicas para fora dos lacos if (log4) then ! **(JP)** fim modificacao do k=k1,m1 rc(k) = rcp(k,i,j) enddo else do k=k1,m1 rc(k) = max(rv(k,i,j)/rvls(k) - 0.999, 0.) enddo endif endif do k=k2,m1-1 vctr1(k) = g/((zt(k+1) - zt(k-1))*rtgt(i,j)) enddo ! **(JP)** fatora expressoes logicas para fora dos lacos if (log2) then ! **(JP)** fim modificacao do k=k2,m1-1 if (rc(k)>0.) then rvii = rvls(k-1) if (k==ki) rvii = rvlsi en2(k,i,j) = vctr1(k)* & ((1. + vctr2(k)*rvls(k)/temp(k))/ & (1. + vctr3(k)*rvls(k)/temp(k)**2)* & ((vctr11(k+1) - vctr11(k-1))/vctr11(k) + & vctr4(k)/temp(k)*(rvls(k+1) - rvii)) - & (rtp(k+1,i,j) - rtp(k-1,i,j))) else en2(k,i,j) = vctr1(k)*((vctr12(k+1)-vctr12(k-1))/vctr12(k) - & (vctr32(k+1) - vctr32(k-1))) endif enddo else do k=k2,m1-1 en2(k,i,j) = vctr1(k)*((vctr12(k+1)-vctr12(k-1))/vctr12(k) - & (vctr32(k+1) - vctr32(k-1))) enddo endif ! **(JP)** remove para fora do laco, permitindo vetorizacao !en2(k1,i,j) = en2(k2,i,j) !en2(nzp,i,j)=en2(nz,i,j) ! **(JP)** fim da modificacao enddo enddo ! **(JP)** removido de dentro do laco do j=ja,jz do i=ia,iz en2(lpw(i,j)-1,i,j) = en2(lpw(i,j),i,j) enddo enddo do j=ja,jz do i=ia,iz en2(nzp,i,j) = en2(nz,i,j) enddo enddo ! **(JP)** fim da modificacao end subroutine bruvais ! ***************************************************************** subroutine mxdefm(m1, m2, m3, ia, iz, ja, jz, ibcon, jd, & vt3dh, vt3di, vt3dj, vt3dk, scr1, scr2, dn0, rtgt, dxt, dyt, lpw_R, mynum) ! +-------------------------------------------------------------+ ! \ this routine calculates the mixing coefficients with a \ ! \ smagorinsky-type deformational based k with an optional \ ! \ unstable brunt-vaisala enhancement and an optional \ ! \ richardson number modification. \ ! +-------------------------------------------------------------+ use mem_scratch, only: & vctr1, & !INTENT(INOUT) vctr2 !INTENT(INOUT) use mem_grid, only: & ngrid, & !INTENT(in) zm, & !INTENT(in) zt, & !INTENT(in) akminvar !INTENT(in) ! For new turb scheme use mem_turb, only: & csx, & !INTENT(in) idiffk, & !INTENT(in) rmin, & !INTENT(out) rmax, & !INTENT(out) zkhkm, & !INTENT(in) csz, & !INTENT(in) akmin !INTENT(in) use rconstants, only: & vonk implicit none ! Arguments: integer, intent(IN) :: m1, m2, m3, ia, iz, ja, jz integer, intent(IN) :: ibcon, jd, mynum !- EHE -> nao sao usadas!!! real, INTENT(INOUT) :: vt3dh(m1,m2,m3), vt3dk(m1,m2,m3), & scr1(m1,m2,m3), scr2(m1,m2,m3) real, INTENT(IN) :: dn0(m1,m2,m3), vt3di(m1,m2,m3), vt3dj(m1,m2,m3) real, INTENT(IN) :: rtgt(m2,m3), dxt(m2,m3) real, INTENT(IN) :: dyt(m2,m3) !- EHE -> nao e' usada!!! real, INTENT(IN) :: lpw_R(m2,m3) ! Local variables: integer :: lpw(m2,m3) integer :: i, j, k, irich, ienfl real :: csx2, sq300, enfl, rchmax, c1, c2, c3, c4, akm, ambda, vkz2 lpw=int(lpw_R) irich = 1 ienfl = 1 csx2 = csx(ngrid)*csx(ngrid) sq300 = 90000. if (idiffk(ngrid)==2 .or. idiffk(ngrid)==3) then rmin = -100. rmax = 1./zkhkm(ngrid) do j=ja,jz do i=ia,iz do k=lpw(i,j),m1-1 vt3dk(k,i,j) = max(min(vt3dj(k,i,j)/& max(vt3di(k,i,j), 1.e-15), rmax), rmin) enddo enddo enddo enfl = float(ienfl) rchmax = 1.0 + 9.0*float(irich) !--(DMK-CCATT)--------------------------------------------------------- !mudanca rams60 nova versao do k = 2, m1 !--(DMK-original)------------------------------------------------------ ! do k = lpw(i,j),m1 !--(DMK-CCATT)--------------------------------------------------------- vctr1(k) = csz(ngrid)*(zm(k) - zm(k-1)) vctr2(k) = vctr1(k)*vctr1(k) enddo endif if (idiffk(ngrid)==1 .or. idiffk(ngrid)==7.or. idiffk(ngrid)==8) then do j=ja,jz do i=ia,iz c2 = 1.0/(dxt(i,j)*dxt(i,j)) c3 = csx2*c2 akm = abs(akmin(ngrid))*0.075*c2**(0.666667) if (AKMIN(ngrid)<0.) then !ml/srf - for new turb scheme !---- !srf-define diferentes AKMINs para melhorar estabilidade !srf-sobre os Andes !!akm=extra2d(5,ngrid)%d2(i,j)*akm akm = akminvar(ngrid)%akmin2d(i,j)*akm endif do k=lpw(i,j),m1-1 scr2(k,i,j) = dn0(k,i,j)*max(akm, c3*sqrt(vt3dh(k,i,j))) enddo enddo enddo elseif (idiffk(ngrid)==2) then do j=ja,jz do i=ia,iz c1 = rtgt(i,j)*rtgt(i,j) c2 = 1.0/(dxt(i,j)*dxt(i,j)) c3 = csx2*c2 akm = abs(akmin(ngrid))*0.075*c2**(0.666667) c4 = vonk*vonk*c1 do k=lpw(i,j),m1-1 ! old csz*dz len scr1(k,i,j) = dn0(k,i,j)*c1*vctr2(k) ! asymptotic vertical scale length from bjorn with modifications: ! c3 is (csx*dx)^2, c1*vctr2(k) is (csz*dz)^2, sq300 is the square ! of 300 meters (used as a limit for horizontal grid spacing influence ! on vertical scale length), ambda is (asymptotic_vertical_length_scale)^2, ! and vkz2 is (vonk*height_above_surface)^2. ambda = max(c1*vctr2(k), min(sq300,c3)) vkz2 = c4*zt(k)*zt(k) scr1(k,i,j) = dn0(k,i,j)*vkz2/(vkz2/ambda + 1)*& (sqrt(vt3di(k,i,j)) + & enfl*sqrt(max(0., -vt3dj(k,i,j))))* & min(rchmax, sqrt(max(0.,(1.-zkhkm(ngrid)*vt3dk(k,i,j))))) scr2(k,i,j) = dn0(k,i,j)*max(akm, c3*sqrt(vt3dh(k,i,j))) vt3dh(k,i,j) = scr1(k,i,j)*zkhkm(ngrid) enddo enddo enddo !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! call friclyr(nzp,nxp,nyp,a(iscr1),a(iustarl),a(itstarl) ! + ,a(iustarw),a(itstarw),a(ipctlnd),a(itheta),a(irtgt)) !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc elseif (idiffk(ngrid)==3) then do j=ja,jz do i=ia,iz c1 = rtgt(i,j)*rtgt(i,j) do k=lpw(i,j),m1-1 scr1(k,i,j) = dn0(k,i,j)*c1*vctr2(k)*(sqrt(vt3dh(k,i,j)) + & enfl*sqrt(max(0.,-vt3dj(k,i,j))))* & min(rchmax, sqrt(max(0.,(1.-zkhkm(ngrid)*vt3dk(k,i,j))))) scr2(k,i,j) = scr1(k,i,j) vt3dh(k,i,j) = scr1(k,i,j)*zkhkm(ngrid) enddo enddo enddo !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! call friclyr(nzp,nxp,nyp,a(iscr1),a(iustarl),a(itstarl) ! + ,a(iustarw),a(itstarw),a(ipctlnd),a(itheta),a(irtgt)) !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc endif end subroutine mxdefm ! ***************************************************************** subroutine klbnd(m1,m2,m3,ibcon,jd,akay,dn0,lpw_R) implicit none integer, INTENT(IN) :: m1 & , m2 & , m3 & , ibcon & , jd real, dimension(m1,m2,m3), INTENT(INOUT) :: akay real, dimension(m1,m2,m3), INTENT(IN) :: dn0 real, dimension(m2,m3), INTENT(IN) :: lpw_R !local variables: integer ::i,j,k,k2 integer,dimension(m2,m3) :: lpw ! boundary conditions on a mixing coefficient lpw=int(lpw_R) do j = 1,m3 do i = 1,m2 k2=lpw(i,j) do k=1,lpw(i,j)-1 akay(k,i,j) = akay(k2,i,j) * dn0(k,i,j) / dn0(k2,i,j) enddo akay(m1,i,j) = akay(m1-1,i,j) * dn0(m1,i,j) & / dn0(m1-1,i,j) enddo enddo if (iand(ibcon,1) .ne. 0) then do j = 1,m3 do k = 1,m1 akay(k,1,j) = akay(k,2,j) enddo enddo endif if (iand(ibcon,2) .ne. 0) then do j = 1,m3 do k = 1,m1 akay(k,m2,j) = akay(k,m2-1,j) enddo enddo endif if (jd .eq. 1) then if (iand(ibcon,4) .ne. 0) then do i = 1,m2 do k = 1,m1 akay(k,i,1) = akay(k,i,2) enddo enddo endif if (iand(ibcon,8) .ne. 0) then do i = 1,m2 do k = 1,m1 akay(k,i,m3) = akay(k,i,m3-1) enddo enddo endif endif return end subroutine klbnd ! ***************************************************************** subroutine mxdefm_tracer(m1, m2, m3, ia, iz, ja, jz, ibcon, jd, & vt3dh, khtr, dn0, dxt, dyt, lpw_R, mynum) !ALF ! +-------------------------------------------------------------+ ! \ this routine calculates the mixing coefficients with a \ ! \ smagorinsky-type deformational based k \ ! +-------------------------------------------------------------+ ! khtr = coef. dif. horizontal para tracers ! frtr = fator de reducao do Akmin dos campos meteorologicos ! use mem_grid, only: & ngrid, & !INTENT(IN) !--(DMK-CCATT-INI)------------------------------------------------------- ngrids !--(DMK-CCATT-FIM)------------------------------------------------------- use mem_turb, only: & csx, & !INTENT(IN) akmin !INTENT(IN) !--(DMK-CCATT-INI)------------------------------------------------------- use monotonic_adv, ONLY: advmnt !INTENT(IN) !--(DMK-CCATT-FIM)------------------------------------------------------- implicit none ! Arguments: integer, intent(in) :: m1, m2, m3, ia, iz, ja, jz, ibcon, jd, mynum real, intent(in) :: lpw_R(m2,m3) real, intent(in) :: dxt(m2,m3), dyt(m2,m3) !dyt nao eh usada real, intent(in) :: vt3dh(m1,m2,m3), dn0(m1,m2,m3) real, intent(inout) :: khtr(m1,m2,m3) ! local variables: integer :: lpw(m2,m3) integer :: i, j, k real :: csx2, c2, c3, akm!, frtr !--(DMK-CCATT-INI)------------------------------------------------------- real :: frtr(ngrids) lpw=int(lpw_R) frtr = 0.1 if(advmnt > 0) frtr = 0.01 !--(DMK-CCATT-OLD)------------------------------------------------------- ! frtr = 0.05 !--(DMK-CCATT-FIM)------------------------------------------------------- csx2 = csx(ngrid)*csx(ngrid) do j=min(1,ja),max(jz,m3) do i=min(1,ia),max(iz,m2) c2 = 1.0/(dxt(i,j)*dxt(i,j)) c3 = csx2*c2 !--(DMK-CCATT-INI)------------------------------------------------------- akm = frtr(ngrid) * abs(akmin(ngrid)) * 0.075 * c2 ** (0.666667) !--(DMK-CCATT-OLD)------------------------------------------------------- ! akm = frtr*abs(akmin(ngrid))*0.075*c2**(0.666667) !--(DMK-CCATT-FIM)------------------------------------------------------- do k=lpw(i,j),m1-1 khtr(k,i,j) = dn0(k,i,j)*max(akm, c3*sqrt(vt3dh(k,i,j))) enddo enddo enddo end subroutine mxdefm_tracer ! ***************************************************************** !srf - OLAM subroutine bruvais_OLAM(m1, m2, m3, ia, iz, ja, jz, theta, rtp, rv, rcp, & pp, pi0, en2, rtgt, lpw_R) use mem_scratch, only: & vctr11, & !INTENT(INOUT) vctr12, & !INTENT(INOUT) vctr32, & !INTENT(INOUT) vctr1, & !INTENT(INOUT) vctr2, & !INTENT(INOUT) vctr3, & !INTENT(INOUT) vctr4 !INTENT(INOUT) use micphys, only: & level, & !INTENT(in) ipris, & !INTENT(in) isnow, & !INTENT(in) igraup, & !INTENT(in) iaggr, & !INTENT(in) ihail !INTENT(in) use mem_grid, only: & zt, & !INTENT(in) nzp, & !INTENT(in) nz, & !INTENT(in) nzpmax !INTENT(in) use rconstants, only: & alvl, & !INTENT(in) rgas, & !INTENT(in) ep, & !INTENT(in) cp, & !INTENT(in) alvi, & !INTENT(in) p00, & !INTENT(in) cpor, & !INTENT(in) g !INTENT(in) implicit none ! Arguments: integer, INTENT(IN) :: m1, m2, m3, ia, iz, ja, jz real, INTENT(IN) :: theta(m1,m2,m3), pp(m1,m2,m3), pi0(m1,m2,m3), & rtp(m1,m2,m3), rcp(m1,m2,m3), rv(m1,m2,m3) real, INTENT(INOUT) :: en2(m1,m2,m3) real, INTENT(IN) :: rtgt(m2,m3) real, INTENT(IN) :: lpw_R(m2,m3) !local variables integer :: i, j, k, iweten, iwdiffk, ki, k2, k1 real :: c1, c2, c3, ci1, ci2, ci3, rvlsi, rvii real :: pi(nzpmax), temp(nzpmax), prt(nzpmax), rvls(nzpmax), rc(nzpmax) ! **(JP)** fatora expressoes logicas para fora dos lacos logical :: log1, log2, log3, log4 integer :: lpw(m2,m3) lpw=int(lpw_R) !------------------------------------------------- !srf- 29/12/2008: adaptado da versao OLAM3.0 !------------------------------------------------- ! calculate brunt-vaisalla frequency squared (en2) !iweten = 1 ! !iwdiffk = 0 !c1 = alvl / rgas !c2 = ep * alvl ** 2 / (cp * rgas) !c3 = alvl / cp !ci1 = alvi / rgas !ci2 = ep * alvi ** 2 / (cp * rgas) !ci3 = alvi / cp ! **(JP)** fatora expressoes logicas para fora dos lacos log1 = level>=1 ! log2 = level .ge. 2 .and. iweten .eq. 1 ! log3 = (ipris .ge. 1 .or. isnow .ge. 1 .or. & ! igraup .ge. 1 .or. iaggr .ge. 1 .or. & ! ihail .ge. 1) .and. level .eq. 3 ! log4 = level .eq. 3 ! **(JP)** fim modificacao ! calculate potential temperature profile do j=ja,jz do i=ia,iz k2 = lpw(i,j) k1 = k2 - 1 do k=k1,m1 vctr11(k) = theta(k,i,j) vctr12(k) = theta(k,i,j) !vctr32(k) = 0. enddo ! **(JP)** fatora expressoes logicas para fora dos lacos if (log1) then ! **(JP)** fim modificacao do k=k1,m1 vctr12(k) = vctr11(k)*(1. + 0.61*rv(k,i,j)) !vctr32(k) = (rtp(k,i,j) - rv(k,i,j)) enddo endif ! check for saturation if level is 2 or greater. ! **(JP)** fatora expressoes logicas para fora dos lacos !if (log2) then ! **(JP)** fim modificacao ! do k = k1,m1 ! pi(k) = (pp(k,i,j) + pi0(k,i,j)) / cp ! temp(k) = theta(k,i,j) * pi(k) ! prt(k) = p00 * pi(k) ** cpor ! enddo ! call mrsl(m1,prt(:),temp(:),rvls(:)) ! do k = k2-1,m1 ! vctr2(k) = c1 ! vctr3(k) = c2 ! vctr4(k) = c3 ! enddo ! ki = m1 + 1 ! ! ! if any ice phase microphysics are activated .... ! ! **(JP)** fatora expressoes logicas para fora dos lacos ! if (log3) then ! ! **(JP)** fim modificacao ! ! find level of -20 c. assume ice saturation above this ! ! level. ! ! do k = k1,m1 ! if (temp(k) .le. 253.16) then ! ki = k ! go to 10 ! endif ! enddo ! ki = m1 + 1 !10 continue ! call mrsi(m1-ki+1,prt(ki),temp(ki),rvls(ki)) !-srf 19/03/2005 !bug: Subscript out of range for array prt ! subscript=0, lower bound=1, upper bound=132, dimension=1 ! tempc < 253 sobre parte da Antartica. ! ki = max(ki,k2) !solucao 1 ! call mrsi(1,prt(ki-1),temp(ki-1),rvlsi) ! if(ki > 1) call mrsi(1,prt(ki-1),temp(ki-1),rvlsi) ! solucao 2 !srf ! ! ! do k = ki,m1 ! vctr2(k) = ci1 ! vctr3(k) = ci2 ! vctr4(k) = ci3 ! enddo ! endif ! ! ! **(JP)** fatora expressoes logicas para fora dos lacos ! if (log4) then ! ! **(JP)** fim modificacao ! do k = k1,m1 ! rc(k) = rcp(k,i,j) ! enddo ! else ! do k = k1,m1 ! rc(k) = max(rv(k,i,j) / rvls(k) - .999,0.) ! enddo ! endif ! !endif do k = k2,m1-1 vctr1(k) = g/((zt(k+1) - zt(k-1))*rtgt(i,j)) enddo ! **(JP)** fatora expressoes logicas para fora dos lacos !if (log2) then ! **(JP)** fim modificacao ! do k = k2,m1-1 ! if (rc(k) .gt. 0.) then ! rvii = rvls(k-1) ! if (k .eq. ki) rvii = rvlsi ! en2(k,i,j) = vctr1(k) * ( & ! (1. + vctr2(k) * rvls(k) / temp(k)) & ! / (1. + vctr3(k) * rvls(k) / temp(k) ** 2) & ! * ((vctr11(k+1) - vctr11(k-1)) / vctr11(k) & ! + vctr4(k) / temp(k) * (rvls(k+1) - rvii)) & ! - (rtp(k+1,i,j) - rtp(k-1,i,j))) ! else ! en2(k,i,j) = vctr1(k)*((vctr12(k+1)-vctr12(k-1)) & ! / vctr12(k) - (vctr32(k+1) - vctr32(k-1))) ! endif ! ! enddo !else do k=k2,m1-1 en2(k,i,j) = vctr1(k)*((vctr12(k+1)-vctr12(k-1))/vctr12(k)) ! / vctr12(k) - (vctr32(k+1) - vctr32(k-1))) enddo !endif ! **(JP)** remove para fora do laco, permitindo vetorizacao !en2(k1,i,j) = en2(k2,i,j) !en2(nzp,i,j)=en2(nz,i,j) ! **(JP)** fim da modificacao enddo enddo ! **(JP)** removido de dentro do laco do j=ja,jz do i=ia,iz en2(lpw(i,j)-1,i,j) = en2(lpw(i,j),i,j) enddo enddo do j=ja,jz do i=ia,iz en2(nzp,i,j) = en2(nz,i,j) enddo enddo ! **(JP)** fim da modificacao end subroutine bruvais_OLAM ! ***************************************************************** !----------------------------------------------------------------------- ! Begin subroutine HCV_driver (IDIFF==7) !----------------------------------------------------------------------- subroutine HCV_driver(m1,m2,m3,ia,iz,ja,jz,mynum & ,up & ,vp & ,theta & ,rv & ,dn0 & ,sflux_t & ,sflux_r & ,sflux_u & ,sflux_v & ,rtgt & ,lpw_R & ,scr1 & ,vt3dh & ,vt3di & ,vt3dj & ,vt3dk ) use mem_grid, only: zt, dzt, dzm, ngrid !INTENT(IN) use mem_turb, only: zkhkm use rconstants, only: & tkmin, & !INTENT(IN) g, & !INTENT(IN) vonk !INTENT(IN) use mem_scratch, only: vctr30 !INTENT(OUT) implicit none real, intent(out) , dimension(m1,m2,m3) :: scr1,vt3dh,vt3di,vt3dj,vt3dk real, intent(in) , dimension(m2,m3) :: rtgt integer, intent(in) :: m1,m2,m3,ia,iz,ja,jz,mynum real, intent(in) , dimension(m1,m2,m3) :: theta, & rv, & up, & vp, & dn0 real, intent(in), dimension(m2,m3) :: lpw_R real, intent(in), dimension(m2,m3) :: sflux_t, & sflux_r, & sflux_u, & sflux_v integer :: i,j,k,k2,kzi,lpw integer,dimension(m2,m3) :: kzi_2d real :: sbf, zl, wstar, ustar, h do j=ja,jz do i=ia,iz lpw=int(lpw_R(i,j)) k2=lpw ! --------------- ! Fill column vector with virtual potential temperature ! --------------- ! do k=k2,m1-1 do k=1,m1-1 vctr30(k) = theta(k,i,j) * (1. + .61 * rv(k,i,j)) enddo ! ! --------------- ! Compute velocity scale for neutra/stable layer [ustar] ! --------------- ustar = max(0.1,sqrt(sqrt(sflux_u(i,j)**2 + sflux_v(i,j)**2))) ! --------------- ! Compute surface buoyancy flux [sbf] and Monin Obukhov height [zl] ! --------------- sbf = (g/vctr30(k2)) & * (sflux_t(i,j) * (1. + .61 * rv(k2,i,j)) & + sflux_r(i,j) * .61 * theta(k2,i,j)) if(sbf .ne. 0.) then zl = (-ustar** 3)/(vonk*sbf) else zl = 500. endif ! ! --------------- ! Compute boundary layer depth [h] ! Compute boundary layer depth [h] for CLE (Zilitinkevich's, 1972) if (zl <= 500.0 .and. zl > 0) then h = 2.4E3 * (ustar**(3.0/2.0)) else !Compute boundary layer depth [h] for CLC - using gradient Richardon number [Ri] call get_richgrad(h,k2,m1,zt,rtgt(i,j),vt3dj(1,i,j),vt3di(1,i,j)& ,vt3dk(1,i,j),zkhkm(ngrid),i,j) endif ! --------------- ! Compute convective velocity scale [wstar] / boundary layer depth [vt2da] ! --------------- wstar = (max(0., sbf * h)) ** (1./3.) ! --------------- ! Find the highest model level zt(k) that is at or below boundary layer ! height h. ! --------------- do k=k2,m1-1 kzi = k if (zt(k+1)*rtgt(i,j) .gt. h) exit enddo kzi_2d(i,j)=kzi ! --------------- ! get KZZ ! --------------- vt3dh(:,i,j)=1.e-3 scr1 (:,i,j)=1.e-3 do k=k2,kzi !****************************************************************************** ! ! Call subroutine kcgcv2 for to calculate KZZ - [scr1(k,i,j] ! !****************************************************************************** call kcgcv2(scr1(k,i,j), ustar, wstar, zl, zt(k)*rtgt(i,j),h, i,j) ! Kmom = scr1(k,i,j) scr1(k,i,j) = max(1.e-3,scr1(k,i,j) * dn0(k,i,j)) ! For Ksclr = Kmom * 3.0 vt3dh(k,i,j)=scr1(k,i,j)!*3.0 - done outside of this subroutine enddo enddo enddo scr1 (1,:,:)=scr1 (2,:,:) vt3dh(1,:,:)=vt3dh(2,:,:) !write(mynum,*),"kzz=",maxval(scr1),minval(scr1),maxval(kzi_2d),maxloc(kzi_2d)& !;call flush(mynum) !i=7 !j=15 !if(mynum==3) then ! do k=1,kzi_2d(i,j)+1 ! write(mynum,*)"k kzz=",k,scr1(k,i,j) ! call flush(3) ! enddo !endif return end subroutine HCV_driver !****************************************************************** ! ! Begin subroutine get_richgrad for calculate boundary layer depth [h] for CLC ! !****************************************************************************** subroutine get_richgrad(h,k2,m1,zt,rtgt,vt3dj,vt3di,vt3dk,zkhkm,i,j) implicit none integer, intent (in) :: k2,m1,i,j real, intent(out) :: h real, intent(in) :: zkhkm,rtgt real, dimension(m1), intent(in) :: zt,vt3dj,vt3di real, dimension(m1), intent(inout) :: vt3dk real :: rmin,rmax,agrich,rchmax integer k,irich irich = 1 rmin = -100. rmax = 1. / zkhkm do k = k2,m1-1 vt3dk(k) = max(min(vt3dj(k) & / max(vt3di(k),1.e-15),rmax),rmin) enddo rchmax = 1.0 + 9.0 * float(irich) h=zt(k2)*rtgt do k = k2,m1-1 agrich=min(rchmax,sqrt(max(0.,(1.-zkhkm*vt3dk(k))))) if(agrich .le. 1.e-5) then h=max(rtgt*zt(k2), rtgt*(zt(k)+zt(k-1))*.5) exit endif enddo return end subroutine get_richgrad ! ***************************************************************** ! Begin subroutine kcgcv2 ! The goal is to provide exchange turbulent coefficients (KZZ) for the B-RAMS system. ! ***************************************************************** SUBROUTINE kcgcv2 (KZZ, ustar, wstar, zl, z, h, i, j) !################################################################ !# # !# PURPOUSE: This routine was designed to run in RMAS, see # !# reference 4. # !# # !# # !# Programer: Haroldo Fraga de Campos Velho # !# # !# Permanet Address: # !# # !# LAC-INPE # !# P.O. Box 515 # !# 12.201-970 - Sao Jose dos Campos (SP) # !# Brasil # !# # !# Fax: +55 012 345-6375 # !# Phone: +55 012 345-6356 # !# E-mail: haroldo@lac.inpe.br # !# Home-page: http://www.lac.inpe.br/~haroldo/ # !# # !# # !# Date: Abril 15, 2000 # !# Last alteration: June 19, 2006 # !# # !# # !# REFERENCES: # !# # !# 1. G.A. Degrazia, H.F. Campos Velho, J.C. Carvalho (1997): # !# "Nonlocal Exchange Coefficients for the Convective # !# Boundary Layer Derived from Spectral Properties", # !# Beitrage zur Physik der Atmosphare, Vol. 70, No. 1, # !# pp. 57-64. # !# # !# 2. G.A. Degrazia, O.L.L. Moraes (1992): "A Model for Eddy # !# Diffusivity in a Stable Boundary Layer", Boundary- # !# Layer Meteorology, Vol. 58, pp. 205-124. # !# # !# 3. G.A. Degrazia, D. Anfossi, J.C. Carvalho, T. Tirabassi, # !# H.F. Campos Velho (2000): "Turbulence Parameterization # !# for PBL Dispersion Models in All Stability Conditions", # !# Atmospheric Environment (to appear) # !# # !# 4. R.A. Pielke, W.R. Cotton, R.L. Walko, C.J. Tremback, W.A. # !# Lyons, L.D. Grasso, M.E. Nicholls, M.D. Moran, D.A. # !# Wesley, T.J. Lee, J. Copeland (1992): "A Comprehensive # !# Meteorological Modeling System - RAMS", Meteorologycal # !# and Atmospheric Physics, Vol. 49, pp. 69-91. # !# # !################################################################ !# # !# SUBROUTINES, FUNCTIONS REQUIRED: none # !# # !# INPUTS PARAMETERS: # !# Z: vertical coordinate # !# H: boundary layer height # !# ZL: Obukhov length # !# USTAR: velocity scale for neutral/stable layer # !# WSTAR: velocity scale for convective layer # !# # !# OUTPUT PARAMETERS: # !# KZZ: Vertical eddy diffusivity # !# # !# # !# LOCAL VARIABLES: # !# LAMB: local Monin-Obukhov length, see reference 2. # !# Q: stability function, Eq. 2.17 in reference 1. # !# PSI: Dissipation function, Eq. 2.19 in reference 1. # !# FC: Coriolis parameter. # !# CORR: Correction factor. # !# ZLM: Average Monin-Obukhov length. # !# KN: Vertical eddy diffusivity from wind shear. # !# # !################################################################ IMPLICIT NONE REAL, INTENT(OUT) :: KZZ !scr1(k,i,j) REAL, INTENT(IN) :: zl, wstar, h, z, ustar INTEGER :: i, j REAL LAMB, Q, PSI, FC, CORR, KN, ZIL, ZH REAL AUX1, AUX2, AUXL REAL ALP1, ALP2 REAL RNU, DEN ! Nondimensional vertical coordinate ZH = min(1.,Z/H) ! Coriolis parameter (see ref. 3) FC = 1.0E-04 ! Contribution from wind shear to turbulence (see Eq(28) in ref. 3) AUX1 = 0.4*(1-ZH)**0.85 * (USTAR*Z) AUX2 = (1.0 + 15.0*(FC*Z/USTAR))**(4.0/3.0) KN = AUX1/AUX2 IF(ABS(ZL).GT.500) THEN KZZ = KN RETURN ENDIF IF(ZL.GT.0.0) THEN ! 1 - Stable parameterization (see ref. 2) ! 1.1 - Set experimental parameters IF(ZL.le.500.0) THEN ALP1 = 1.5 ALP2 = 1.0 ENDIF ! 1.2 - Calculate the local Monin-Obukhov lenght (see Eq(3) in ref. 2) AUXL = (3.0*ALP1)/2.0 - ALP2 LAMB = ABS(ZL)*(1.0 - ZH)**AUXL ! 1.3 - Stable eddy diffusivity (see Eq(29) in ref. 3) RNU = 0.4 * (1.-ZH)**(3.0/4.0) DEN = 1. + 3.7* (Z/LAMB) KZZ = RNU/DEN * (USTAR*Z) RETURN ELSE ! 2 - Convective parameterization ! The stability parameter ZIL is the ratio between the ! convective boundary layer (CBL) height and Monin-Obukhov ! length. For a full convective CBL: 50 < ZIL < 360. ! See comment in page 60 in reference 1. ! Correction factor (see Eq. (14) in ref. 3) CORR = SQRT(0.01*(H /(-ZL))) !See Eq. 14 4 ! Here will be used the same value in KAMM model. ZIL = 136.0 ! 2.1 - Stability function: Eq. (2.17) in reference 1. Q = 1.0 - EXP(-4*ZH) - 0.0003*EXP(8*ZH) ! 2.2 - Eddy diffusivity: Eq. (2.18) in reference 1. AUX2 = 4.0/3.0 ! 2.3 - Eddy diffusivity: Eq. (26) in reference 3. KZZ = 0.16* (WSTAR*H) * CORR * Q**AUX2 ! Shear effect contribution (see Eq. (15), ref. 3) KZZ = KZZ + KN ENDIF END SUBROUTINE kcgcv2 !-----------------------------------------------------------