! $Id: step3d_t.F 1576 2014-07-03 16:15:16Z gcambon $ ! !====================================================================== ! CROCO is a branch of ROMS developped at IRD and INRIA, in France ! The two other branches from UCLA (Shchepetkin et al) ! and Rutgers University (Arango et al) are under MIT/X style license. ! CROCO specific routines (nesting) are under CeCILL-C license. ! ! CROCO website : http://www.croco-ocean.org !====================================================================== ! #include "cppdefs.h" #ifdef SOLVE3D subroutine step3d_t (tile) ! implicit none integer tile # ifdef TRACERS integer trd, omp_get_thread_num # include "param.h" # include "private_scratch.h" # include "ocean3d.h" # include "compute_tile_bounds.h" trd=omp_get_thread_num() call step3d_t_tile (Istr,Iend,Jstr,Jend, & A2d(1,1,trd), A2d(1,2,trd), A2d(1,3,trd), & A2d(1,4,trd), A2d(1,5,trd), A2d(1,6,trd), & A2d(1,7,trd), A2d(1,8,trd), A2d(1,9,trd), & A3d(1,1,trd)) # endif /* TRACERS */ return end # ifdef TRACERS subroutine step3d_t_tile (Istr,Iend,Jstr,Jend, & FX,FE, WORK, FC,CF,BC,DC,EC,GC, swdk) # ifdef SUBSTANCE USE comsubstance, ONLY : nv_grav, nvp, ws_part, nv_sand # ifdef MUSTANG USE comMUSTANG, ONLY : flx_s2w_CROCO, flx_w2s_CROCO USE comMUSTANG, ONLY : flx_w2s_sum_CROCO USE comMUSTANG, ONLY : corflux, corfluy USE comsubstance, ONLY : irkm_var_assoc, nvpc, nvp # ifdef key_sand2D USE comsubstance, ONLY : l_subs2D # endif # endif # endif implicit none # include "param.h" # include "grid.h" # include "ocean3d.h" # include "mixing.h" # include "climat.h" # include "scalars.h" # include "sources.h" # include "forces.h" # ifdef DIAGNOSTICS_TS # include "diagnostics.h" # endif # ifdef DIAGNOSTICS_PV # include "diags_pv.h" # endif integer Istr,Iend,Jstr,Jend, itrc, i,j,k, indx, kmld & ,imin,imax,jmin,jmax,iAkt,nadv # if defined PSOURCE || defined PSOURCE_MASS & ,is # endif # ifdef PSOURCE_MASS & ,iis,jjs # endif # ifdef SUBSTANCE integer isubstep,nbsubstep real :: cflm,cfl_loc,dz_cflm,ws_cflm # endif # ifdef AGRIF # include "zoom.h" # endif real FX(PRIVATE_2D_SCRATCH_ARRAY), & FE(PRIVATE_2D_SCRATCH_ARRAY), cff, & WORK(PRIVATE_2D_SCRATCH_ARRAY), epsil, & FC(PRIVATE_1D_SCRATCH_ARRAY,0:N), & CF(PRIVATE_1D_SCRATCH_ARRAY,0:N), & BC(PRIVATE_1D_SCRATCH_ARRAY,0:N), & DC(PRIVATE_1D_SCRATCH_ARRAY,0:N), & EC(PRIVATE_1D_SCRATCH_ARRAY,0:N), & GC(PRIVATE_1D_SCRATCH_ARRAY,0:N), & swdk(PRIVATE_2D_SCRATCH_ARRAY,0:N) real cff1,cff2,gama,dRz,hbltmp,sig,dXmax,dEmax, & dpth,smax,amax,amaxx # ifdef TS_VADV_FCT real tmax,tmin,Ppos,Qpos,Pneg,Qneg,lmdmax # endif parameter (epsil=1.E-16) ! # ifdef MASKING # define SWITCH * # else # define SWITCH ! # endif ! # if defined TS_HADV_UP5 || defined TS_HADV_C6 || \ defined TS_HADV_WENO5 || defined BIO_HADV_WENO5 || \ defined TS_VADV_WENO5 ! !-------------------------------------------------------------------- ! Definition of flux operators: 1st, 2nd, 3rd, 4th, 5th or 6th order, ! used in UP5 and C6 advection schemes (and order degradation near ! land masks). cdiff is part of laplacian diffusion in flux1 (used ! near mask): ! 0 --> flux1=flux2 (second order C2 advection scheme) ! 1 --> flux1 gives 1st order monotonic UP1 advection scheme !-------------------------------------------------------------------- ! REAL :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2 REAL :: ua, vel, cdiff, cdif REAL :: flux1, flux2, flux3, flux4, flux5, flux6 REAL :: flx2, flx3, flx4, flx5 REAL :: mask0, mask1, mask2, mask3 flux2(q_im1, q_i, ua, cdiff) = 0.5*( q_i + q_im1 ) flux1(q_im1, q_i, ua, cdiff) = flux2(q_im1, q_i, ua, cdiff) - & 0.5*cdiff*sign(1.,ua)*(q_i-q_im1) flux4(q_im2, q_im1, q_i, q_ip1, ua) = & ( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0 flux3(q_im2, q_im1, q_i, q_ip1, ua) = & flux4(q_im2, q_im1, q_i, q_ip1, ua) + & sign(1.,ua)*((q_ip1 - & q_im2)-3.*(q_i-q_im1))/12.0 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & ( 37.*(q_i+q_im1) - 8.*(q_ip1+q_im2) & +(q_ip2+q_im3) )/60.0 flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) & -sign(1.,ua)*( & (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0 # endif # if defined TS_HADV_WENO5 || defined BIO_HADV_WENO5 || defined TS_VADV_WENO5 REAL :: flux3_weno, flux5_weno # endif ! # include "compute_auxiliary_bounds.h" nadv = 3 !<-- do not remove: used in compute_tracer_fluxes.h ! !====================================================================== ! ! Compute horizontal advection ! !====================================================================== ! do itrc=1,NT do k=1,N # ifdef BIO_HADV_WENO5 # if defined TEMPERATURE && defined SALINITY if (itrc .le. 2) then # elif defined TEMPERATURE || defined SALINITY if (itrc .eq. 1) then # endif # endif ! !---------------------------------------------------------- ! Compute fluxes FX and FE !---------------------------------------------------------- ! # undef PREDICTOR # include "compute_horiz_tracer_fluxes.h" ! !---------------------------------------------------------- ! WENO5 advection for passive tracers (biology ...) !---------------------------------------------------------- ! # ifdef BIO_HADV_WENO5 else !<-- if (itrc .gt. 2) # define FLUX5 flux5_weno # define FLUX3 flux3_weno # define FLUX2 flux1 # define UP5_MASKING ! cdif=1. # include "t3dadv_order5.h" ! # undef FLUX5 # undef FLUX3 # undef FLUX2 # undef UP5_MASKING endif # endif ! # if defined MUSTANG && defined MUSTANG_CORFLUX if (k.eq.1 .and. & (itrc.ge.itsubs1+nv_grav .and. & itrc.le.itsubs1+nv_grav+nv_sand-1) ! isand1, isand2 & .or. & (itrc.ge.itsubs1+nvpc .and. & itrc.le.itsubs1+nvp-1) ! nvpc+1, nvp & ) then do j = Jstr, Jend do i = Istr, Iend+1 FX(i, j) = FX(i, j) * & corflux(itrc-itsubs1+1, i, j) enddo enddo do j = Jstr, Jend+1 do i = Istr, Iend FE(i, j) = FE(i, j) * & corfluy(itrc-itsubs1+1, i, j) enddo enddo endif # endif ! !---------------------------------------------------------- ! Apply point sources for river runoff simulations !---------------------------------------------------------- ! # ifdef PSOURCE do is=1,Nsrc # ifdef MPI i=Isrc_mpi(is,mynode) j=Jsrc_mpi(is,mynode) # else i=Isrc(is) j=Jsrc(is) # endif if (Istr.le.i .and. i.le.Iend+1 & .and. Jstr.le.j .and. j.le.Jend+1) then if (Dsrc(is).eq.0) then if (Lsrc(is,itrc)) then FX(i,j)=Huon(i,j,k)*Tsrc(is,k,itrc) # if defined MUSTANG && defined key_sand2D ! To do: sum FX for k=1 # endif # ifdef MASKING else if (rmask(i,j).eq.0 .and. rmask(i-1,j).eq.1) then FX(i,j)=Huon(i,j,k)*t(i-1,j,k,3,itrc) elseif (rmask(i,j).eq.1. .and. rmask(i-1,j).eq.0) then FX(i,j)=Huon(i,j,k)*t(i ,j,k,3,itrc) endif # endif endif elseif (Dsrc(is).eq.1) then if (Lsrc(is,itrc)) then FE(i,j)=Hvom(i,j,k)*Tsrc(is,k,itrc) # ifdef MASKING else if (rmask(i,j).eq.0 .and. rmask(i,j-1).eq.1) then FE(i,j)=Hvom(i,j,k)*t(i,j-1,k,3,itrc) elseif (rmask(i,j).eq.1 .and. rmask(i,j-1).eq.0) then FE(i,j)=Hvom(i,j,k)*t(i,j ,k,3,itrc) endif # endif endif endif endif enddo # endif /* PSOURCE */ ! !---------------------------------------------------------- ! Finalize horizontal advection: compute flux divergences !---------------------------------------------------------- ! do j=Jstr,Jend do i=Istr,Iend # if !(defined MUSTANG && defined key_sand2D) t(i,j,k,nnew,itrc)=Hz_bak(i,j,k)*t(i,j,k,nstp,itrc) & -dt*pm(i,j)*pn(i,j)*( FX(i+1,j)-FX(i,j) & +FE(i,j+1)-FE(i,j)) # else if (l_subs2D(itrc-itsubs1+1)) then if (k.eq.1) then t(i,j,k,nnew,itrc)= & ( Hz_bak(i,j,k)*t(i,j,k,nstp,itrc) & -dt*pm(i,j)*pn(i,j)*(FX(i+1,j)-FX(i,j) & +FE(i,j+1)-FE(i,j)) & +dt*flx_s2w_CROCO(i,j,itrc) ) & /(Hz(i,j,k)+dt*flx_w2s_CROCO(i,j,itrc)) flx_w2s_sum_CROCO(i,j,itrc)=flx_w2s_CROCO(i,j,itrc) & *dt*t(i,j,1,nnew,itrc) endif else t(i,j,k,nnew,itrc)= & Hz_bak(i,j,k)*t(i,j,k,nstp,itrc) & -dt*pm(i,j)*pn(i,j)*( FX(i+1,j)-FX(i,j) & +FE(i,j+1)-FE(i,j)) endif # endif enddo enddo !--> discard FX,FE # if defined AGRIF && defined AGRIF_CONSERV_TRA if (Agrif_Root()) then MYFX(IstrR:IendR,JstrR:JendR,k,itrc)= & dt*FX(IstrR:IendR,JstrR:JendR) MYFY(IstrR:IendR,JstrR:JendR,k,itrc)= & dt*FE(IstrR:IendR,JstrR:JendR) else MYFX(IstrR:IendR,JstrR:JendR,k,itrc)= & MYFX(IstrR:IendR,JstrR:JendR,k,itrc) & +dt*FX(IstrR:IendR,JstrR:JendR) MYFY(IstrR:IendR,JstrR:JendR,k,itrc)= & MYFY(IstrR:IendR,JstrR:JendR,k,itrc) & +dt*FE(IstrR:IendR,JstrR:JendR) endif # endif ! !---------------------------------------------------------- ! Store diagnostic transport term !---------------------------------------------------------- ! # ifdef DIAGNOSTICS_TS do j=Jstr,Jend do i=Istr,Iend TXadv(i,j,k,itrc)=-(FX(i+1,j)-FX(i,j)) SWITCH rmask(i,j) TYadv(i,j,k,itrc)=-(FE(i,j+1)-FE(i,j)) SWITCH rmask(i,j) enddo enddo ! # if defined TS_HADV_WENO5 || defined TS_HADV_UP5 || defined TS_HADV_UP3 ! ! Compute centered scheme fluxes FX and FE to separate ! advection and diffusion parts of the transport term ! # undef PREDICTOR # if defined TS_HADV_WENO5 # define TS_HADV_C6 # undef TS_HADV_WENO5 # include "compute_horiz_tracer_fluxes.h" # undef TS_HADV_C6 # define TS_HADV_WENO5 # elif defined TS_HADV_UP5 # define TS_HADV_C6 # undef TS_HADV_UP5 # include "compute_horiz_tracer_fluxes.h" # undef TS_HADV_C6 # define TS_HADV_UP5 # elif defined TS_HADV_UP3 # define TS_HADV_C4 # undef TS_HADV_UP3 # include "compute_horiz_tracer_fluxes.h" # undef TS_HADV_C4 # define TS_HADV_UP3 # endif do j=Jstr,Jend do i=Istr,Iend THmix(i,j,k,itrc)=( TXadv(i,j,k,itrc)+TYadv(i,j,k,itrc) & +(FX(i+1,j)-FX(i,j)) & +(FE(i,j+1)-FE(i,j)) ) SWITCH rmask(i,j) TXadv(i,j,k,itrc)=-(FX(i+1,j)-FX(i,j)) SWITCH rmask(i,j) TYadv(i,j,k,itrc)=-(FE(i,j+1)-FE(i,j)) SWITCH rmask(i,j) enddo enddo # else do j=Jstr,Jend do i=Istr,Iend THmix(i,j,k,itrc)=0. enddo enddo # endif # endif /* DIAGNOSTICS_TS */ enddo !<-- k enddo !<-- itrc ! !---------------------------------------------------------- ! Compute fraction of the solar shortwave flux "swdk" ! penetrating to grid level depth (at vertical w-points). ! (swdk is used later in this routine) !---------------------------------------------------------- ! # if defined LMD_SKPP || defined LMD_BKPP || defined GLS_MIXING # define wrk1 FX # define wrk2 FE do k=1,N-1 do j=Jstr,Jend do i=Istr,Iend wrk1(i,j)=z_w(i,j,k)-z_w(i,j,N) enddo enddo call lmd_swfrac_tile (Istr,Iend,Jstr,Jend,1.,wrk1,wrk2) do j=Jstr,Jend do i=Istr,Iend swdk(i,j,k)=wrk2(i,j) enddo enddo enddo # undef wrk1 # undef wrk2 # endif /* LMD_SKPP || LMD_BKPP || GLS_MIXING */ ! !---------------------------------------------------------- ! Add tracer divergence due to cell-centered point sources !---------------------------------------------------------- ! #ifdef PSOURCE_MASS do itrc=1,NT do is=1,Nsrc # ifdef MPI iis=Isrc_mpi(is,mynode) jjs=Jsrc_mpi(is,mynode) # else iis=Isrc(is) jjs=Jsrc(is) # endif if (Istr.le.iis .and. iis.le.Iend+1 & .and. Jstr.le.jjs .and. jjs.le.Jend+1) then cff1=dt*pm(iis,jjs)*pn(iis,jjs) do k=1,N if (Lsrc(is,itrc)) then cff2=Tsrc(is,k,itrc) else cff2=t(iis,jjs,k,3,itrc) endif t(iis,jjs,k,nnew,itrc)=t(iis,jjs,k,nnew,itrc)+ & cff1*Qsrc(is,k)*cff2 enddo endif enddo ! source enddo ! itrc # endif /* PSOURCE_MASS */ ! !====================================================================== ! ! Compute vertical advection ! !====================================================================== ! do j=Jstr,Jend do itrc=1,NT ! !---------------------------------------------------------- ! Compute vertical fluxes FC !---------------------------------------------------------- ! # undef PREDICTOR # include "compute_vert_tracer_fluxes.h" # ifdef TS_HADV_TEST FC(Istr:Iend,0:N) = 0.0 # endif !---------------------------------------------------------- ! Apply flux correction for monotonicity !---------------------------------------------------------- ! # ifdef TS_VADV_FCT # include "FCT.h" # endif ! !---------------------------------------------------------- ! Apply vertical advective flux divergence !---------------------------------------------------------- ! # if defined MUSTANG && defined key_sand2D if (.not.l_subs2D(itrc-itsubs1+1)) then # endif do k=1,N do i=Istr,Iend t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc) & -CF(i,0)*(FC(i,k)-FC(i,k-1)) enddo enddo !--> discard FC # ifdef DIAGNOSTICS_TS do k=1,N do i=Istr,Iend TVadv(i,j,k,itrc) = -(FC(i,k)-FC(i,k-1)) & SWITCH rmask(i,j) enddo enddo # endif /* DIAGNOSTICS_TS */ # if defined MUSTANG && defined key_sand2D endif # endif ! ! !====================================================================== ! ! Compute surface and Bottom forcing ! !====================================================================== ! !---------------------------------------------------------- ! Add top and bottom fluxes !---------------------------------------------------------- ! do i=Istr,Iend FC(i,N)=dt*stflx(i,j,itrc) FC(i,0)=-dt*btflx(i,j,itrc) enddo ! !---------------------------------------------------------- ! Add solar radiation flux in temperature equation ! Also compute the nonlocal transport flux for unstable ! (convective) forcing conditions into matrix DC when using ! the Large et al. 1994 KPP scheme. !---------------------------------------------------------- ! # ifdef TEMPERATURE if (itrc.eq.itemp) then do k=1,N-1 do i=Istr,Iend FC(i,k)=0. # if defined LMD_SKPP || defined LMD_BKPP \ || defined GLS_MIXING & +dt*srflx(i,j)*swdk(i,j,k) # ifdef LMD_NONLOCAL & -dt*ghats(i,j,k)*(stflx(i,j,itemp)-srflx(i,j)) # endif # endif enddo enddo endif # endif /* TEMPERATURE */ # ifdef SALINITY if (itrc.eq.isalt) then do k=1,N-1 do i=Istr,Iend FC(i,k)=0. # if defined LMD_SKPP || defined LMD_BKPP # ifdef LMD_NONLOCAL & -dt*ghats(i,j,k)*stflx(i,j,isalt) # endif # endif enddo enddo endif # endif /* SALINITY */ ! # if defined SALINITY || defined TEMPERATURE # if defined SALINITY && defined TEMPERATURE if (itrc.eq.itemp .or. itrc.eq.isalt) then # elif defined SALINITY if (itrc.eq.isalt) then # elif defined TEMPERATURE if (itrc.eq.itemp) then # endif do k=1,N do i=Istr,Iend t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)+FC(i,k ) & -FC(i,k-1) # if defined DIAGNOSTICS_TS || defined DIAGNOSTICS_PV TForc(i,j,k,itrc)=(FC(i,k)-FC(i,k-1)) & /(dt*pm(i,j)*pn(i,j)) & SWITCH rmask(i,j) # endif /* DIAGNOSTICS_TS */ enddo enddo # if defined DIAGNOSTICS_TS || defined DIAGNOSTICS_PV else do k=1,N do i=Istr,Iend TForc(i,j,k,itrc)=0. enddo enddo # endif /* DIAGNOSTICS_TS */ endif # endif /* defined SALINITY || TEMPERATURE */ ! ! !====================================================================== ! ! Compute vertical mixing (implicit step) ! ! If lateral rotated diffusion is chosen, this step is performed ! after lateral diffusion to include implicit correction to lateral ! diffusion ! ! If MUSTANG is used, mixing is done here for all passive tracers ! but sediment tracers include a settling process with substeps ! !====================================================================== ! # ifdef VADV_ADAPT_IMP do i=Istr,Iend DC(i,0)=dt*pn(i,j)*pm(i,j) enddo # endif # ifdef TS_MIX_IMP # if defined TEMPERATURE && defined SALINITY if (itrc .le. 2) then ! active tracers # elif defined TEMPERATURE || defined SALINITY if (itrc .le. 1) then ! active tracers # endif do k=1,N do i=istr,iend t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc) / Hz(i,j,k) enddo enddo else ! passive tracers ! (sed. and others) # ifdef MUSTANG ! if (itrc.ge.itsubs1+nv_grav .and. & itrc.le.itsubs1+nvp) then ! sediment tracers # ifdef key_sand2D if (.not.l_subs2D(itrc-itsubs1+1)) then # endif # undef TS_MIX_IMP # define CD EC # define NBSUBSTEP GC # include "t3dmix_tridiagonal_settling.h" # undef NBSUBSTEP # undef CD # define TS_MIX_IMP # ifdef key_sand2D endif # endif else ! other tracers # undef TS_MIX_IMP # include "t3dmix_tridiagonal.h" # define TS_MIX_IMP endif ! # else /* MUSTANG */ ! ! passive tracers # ifdef BIO_HADV_WENO5 # undef TS_MIX_IMP # include "t3dmix_tridiagonal.h" # define TS_MIX_IMP # else do k=1,N do i=istr,iend t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc) / Hz(i,j,k) enddo enddo # endif /* BIO_HADV_WENO5 */ # endif /* MUSTANG */ ! endif ! itrc <= 2 ! # else /* TS_MIX_IMP */ # ifndef MUSTANG ! ! all tracers # include "t3dmix_tridiagonal.h" ! # else if (itrc.ge.itsubs1+nv_grav .and. & itrc.le.itsubs1+nvp) then ! sediment tracers # ifdef key_sand2D if (.not.l_subs2D(itrc-itsubs1+1)) then # endif # define CD EC # define NBSUBSTEP GC # include "t3dmix_tridiagonal_settling.h" # undef NBSUBSTEP # undef CD # ifdef key_sand2D endif # endif else ! all other tracers # include "t3dmix_tridiagonal.h" endif # endif /* MUSTANG */ # endif /* TS_MIX_IMP */ ! ! Test of constant tracers ! # ifdef CONST_TRACERS do k=1,N do i=Istr,Iend t(i,j,k,nnew,itrc)=t(i,j,k,nstp,itrc) enddo enddo # if defined DIAGNOSTICS_TS || defined DIAGNOSTICS_PV do k=1,N do i=Istr,Iend TVmix(i,j,k,itrc)=0.0 # ifdef MASKING & * rmask(i,j) # endif enddo enddo # endif /* DIAGNOSTICS_TS */ # endif /* CONST_TRACERS */ enddo ! <-- itrc ! !---------------------------------------------------------- ! Compute vertical gradient in vertical T-diffusion ! coefficient for floats random walk. !---------------------------------------------------------- ! # if defined FLOATS && defined RANDOM_WALK do k=1,N do i=Istr,Iend dAktdz(i,j,k)=(Akt(i,j,k,1)-Akt(i,j,k-1,1))/Hz(i,j,k) enddo enddo # endif enddo ! <-- j ! !====================================================================== ! Set lateral boundary conditions; nudge toward tracer climatology; ! apply land-sea mask and exchange periodic boundary conditions. !====================================================================== ! do itrc=1,NT call t3dbc_tile (Istr,Iend,Jstr,Jend, nnew,itrc, WORK) enddo ! !====================================================================== ! Compute biological fluxes !====================================================================== ! # ifdef BIOLOGY call biology_tile (Istr,Iend,Jstr,Jend) # endif ! !====================================================================== ! Compute inverse vertical Rho gradients for use in t3dmix !====================================================================== ! # if !defined TS_MIX_S && \ (defined TS_DIF4 || defined TS_DIF2 || defined SPONGE_DIF2) # ifdef TS_MIX_GEO ! ! --- Geopotential diffusion: idRz contains idZ values --- ! do k=1,N-1 do j=jstr,jend do i=istr,iend idRz(i,j,k) = 1./(z_r(i,j,k+1)-z_r(i,j,k)) enddo enddo enddo do j=jstr,jend do i=istr,iend idRz(i,j,N) = 0. idRz(i,j,0) = 0. enddo enddo # elif defined TS_MIX_ISO ! ! --- Isopycnal diffusion: idRz is inverse vertical rho gradient --- ! do k=1,N-1 do j=jstr,jend do i=istr,iend # ifdef SPLIT_EOS dpth=z_w(i,j,N)-0.5*(z_r(i,j,k+1)+z_r(i,j,k)) dRz =rho1(i,j,k+1)-rho1(i,j,k) & +(qp1(i,j,k+1)- qp1(i,j,k)) & *dpth*(1.- qp2*dpth) # else dRz = rho1(i,j,k+1)-rho1(i,j,k) # endif cff = 1./min(dRz, -1.E-14) ! minimum stratification imposed cff1 = 1./(z_r(i,j,k+1)-z_r(i,j,k)) # ifdef LMD_SKPP2005 hbltmp=hbls(i,j,3-nstp) # elif defined LMD_SKPP1994 || defined GLS_MIXING hbltmp=hbl(i,j) # else hbltmp=10. # endif gama = 1. ! mixing becomes sig = (z_w(i,j,N)-z_w(i,j,k))/max(hbltmp,10.) ! isosigma within if (sig .lt. 1.) gama = sig*sig*(3.-2.*sig) ! PBL depth (>10m) dXmax = max(abs(dRdx(i,j,k )),abs(dRdx(i+1,j,k )), & abs(dRdx(i,j,k+1)),abs(dRdx(i+1,j,k+1)),1E-14) dEmax = max(abs(dRde(i,j,k )),abs(dRde(i,j+1,k )), & abs(dRde(i,j,k+1)),abs(dRde(i,j+1,k+1)),1E-14) smax=min(Rslope_max,Gslope_max*min(pm(i,j),pn(i,j))/cff1) ! idRz(i,j,k) = max(cff, -smax*gama*cff1 / dXmax, & -smax*gama*cff1 / dEmax ) enddo enddo enddo do j=jstr,jend do i=istr,iend idRz(i,j,N) = 0. idRz(i,j,0) = 0. enddo enddo # endif /* TS_MIX_ISO */ ! # if defined EW_PERIODIC || defined NS_PERIODIC || defined MPI call exchange_w3d_tile ( istr,iend,jstr,jend, & idRz(START_2D_ARRAY,0) ) # endif ! # endif /* TS_DIF2 || TS_DIF4 || SPONGE_DIF2 */ ! ! do itrc=1,NT !<-- itrc ! !====================================================================== ! Compute Nudging terms !====================================================================== ! # if defined TNUDGING && defined TCLIMATOLOGY # ifdef ZONAL_NUDGING if (iic.eq.ntstart .or. mod(iic,10).eq.0) then call zonavg_3d(istr,iend,jstr,jend, & t(START_2D_ARRAY,1,nnew,itrc),tzon(:,:,itrc)) endif if (iic.eq.ntstart) then call zonavg_3d(istr,iend,jstr,jend, & tclm(START_2D_ARRAY,1,itrc),tclmzon(:,:,itrc)) endif # endif /* ZONAL_NUDGING */ do k=1,N do j=Jstr,Jend do i=Istr,Iend t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc) & +dt*Tnudgcof(i,j,k,itrc)* # ifdef ZONAL_NUDGING & (tclmzon(j,k,itrc)-tzon(j,k,itrc)) # else & (tclm(i,j,k,itrc)-t(i,j,k,nnew,itrc)) # endif /* ZONAL_NUDGING */ # if defined DIAGNOSTICS_TS || defined DIAGNOSTICS_PV TForc(i,j,k,itrc)=(TForc(i,j,k,itrc) & +Tnudgcof(i,j,k,itrc) & *(tclm(i,j,k,itrc) # ifdef ZONAL_NUDGING & -tzon(j,k,itrc)) # else & -t(i,j,k,nnew,itrc)) # endif /* ZONAL_NUDGING */ & *(Hz(i,j,k)/(pm(i,j)*pn(i,j)))) & SWITCH rmask(i,j) # endif /* DIAGNOSTICS_TS */ enddo enddo enddo # endif /* TNUDGING && TCLIMATOLOGY */ ! !====================================================================== ! Compute the tendency term of tracer diagnostics ! Divide all diagnostic terms by the cell volume ! (Hz(i,j,k,itrc)/(pm(i,j).*pn(i,j)). There after the unit ! of diagnostic terms will be: (unit of tracers)* s-1. ! ! Note: the Horizontal mixing term is computed in t3dmix ! where Trate is updated accordingly !====================================================================== ! # if defined DIAGNOSTICS_TS || defined DIAGNOSTICS_PV do k=1,N do j=Jstr,Jend do i=Istr,Iend Trate(i,j,k,itrc)=(Hz(i,j,k)*t(i,j,k,nnew,itrc) & -Hz_bak(i,j,k)*t(i,j,k,nstp,itrc)) & /(dt*pm(i,j)*pn(i,j)) & SWITCH rmask(i,j) ! # ifdef DIAGNOSTICS_TS_ADV TXadv(i,j,k,itrc)=TXadv(i,j,k,itrc)+ & t(i,j,k,3,itrc)*(Huon(i+1,j,k)-Huon(i,j,k)) TYadv(i,j,k,itrc)=TYadv(i,j,k,itrc)+ & t(i,j,k,3,itrc)*(Hvom(i,j+1,k)-Hvom(i,j,k)) TVadv(i,j,k,itrc)=TVadv(i,j,k,itrc)- & t(i,j,k,3,itrc)*(Huon(i+1,j,k)-Huon(i,j,k)+ & Hvom(i,j+1,k)-Hvom(i,j,k)) # endif /* DIAGNOSTICS_TS_ADV */ ! cff=pm(i,j)*pn(i,j)/Hz(i,j,k) # ifdef DIAGNOSTICS_TS TXadv(i,j,k,itrc)=TXadv(i,j,k,itrc)*cff TYadv(i,j,k,itrc)=TYadv(i,j,k,itrc)*cff TVadv(i,j,k,itrc)=TVadv(i,j,k,itrc)*cff # endif TForc(i,j,k,itrc)=TForc(i,j,k,itrc)*cff Trate(i,j,k,itrc)=Trate(i,j,k,itrc)*cff TVmix(i,j,k,itrc)=TVmix(i,j,k,itrc)*cff THmix(i,j,k,itrc)=THmix(i,j,k,itrc)*cff # if defined DIAGNOSTICS_TSVAR && !defined SPONGE_DIF2 \ && !defined TS_DIF2 && !defined TS_DIF4 # include "finalize_diagnostics_tsadv.h" # endif enddo enddo enddo ! !---------------------------------------------------------- ! Compute tracer diagnostics averaged over the MLD !---------------------------------------------------------- ! # ifdef DIAGNOSTICS_TS_MLD # define T_mld_nnew FX # define T_mld_nstp FE do j=Jstr,Jend do i=Istr,Iend TXadv_mld(i,j,itrc)=0. TYadv_mld(i,j,itrc)=0. TVadv_mld(i,j,itrc)=0. TVmix_mld(i,j,itrc)=0. THmix_mld(i,j,itrc)=0. TForc_mld(i,j,itrc)=0. Trate_mld(i,j,itrc)=0. Tentr_mld(i,j,itrc)=0. T_mld_nnew(i,j)=0. T_mld_nstp(i,j)=0. enddo enddo do j=Jstr,Jend do i=Istr,Iend # ifdef LMD_SKPP kmld=kbl(i,j) # else kmld=N-5 # endif do k=N,kmld,-1 cff=Hz(i,j,k)/(z_w(i,j,N)-z_w(i,j,kmld-1)) TXadv_mld(i,j,itrc)=TXadv_mld(i,j,itrc)+ & TXadv(i,j,k,itrc)*cff TYadv_mld(i,j,itrc)=TYadv_mld(i,j,itrc)+ & TYadv(i,j,k,itrc)*cff TVadv_mld(i,j,itrc)=TVadv_mld(i,j,itrc)+ & TVadv(i,j,k,itrc)*cff TVmix_mld(i,j,itrc)=TVmix_mld(i,j,itrc)+ & TVmix(i,j,k,itrc)*cff TForc_mld(i,j,itrc)=TForc_mld(i,j,itrc)+ & TForc(i,j,k,itrc)*cff T_mld_nnew(i,j) =T_mld_nnew(i,j)+ & t(i,j,k,nnew,itrc)*cff enddo enddo enddo # if (!defined TS_DIF2 && !defined TS_DIF4 && !defined SPONGE_DIF2) \ || defined DIAGNOSTICS_DEBUG ! ! Compute entrainement/detrainement term. ! If diffusion terms are computed in t3dmix routines, entrainement ! terms are also computed in those routines ! do j=Jstr,Jend do i=Istr,Iend # if defined LMD_SKPP || defined GLS_MIXING if (kbl_nstp(i,j).eq.0) kbl_nstp(i,j)=kbl(i,j) kmld=kbl_nstp(i,j) # else kmld=N-5 # endif do k=N,kmld,-1 cff=Hz_bak(i,j,k)/(z_w(i,j,N)-z_w(i,j,kmld-1)) T_mld_nstp(i,j)=T_mld_nstp(i,j)+ & t(i,j,k,nstp,itrc)*cff enddo if (itrc .eq. NT) kbl_nstp(i,j)=kbl(i,j) enddo enddo do j=Jstr,Jend do i=Istr,Iend Trate_mld(i,j,itrc)=(T_mld_nnew(i,j)-T_mld_nstp(i,j))/dt Tentr_mld(i,j,itrc)=Trate_mld(i,j,itrc)- & TXadv_mld(i,j,itrc)- & TYadv_mld(i,j,itrc)- & TVadv_mld(i,j,itrc)- & TVmix_mld(i,j,itrc)- & TForc_mld(i,j,itrc) enddo enddo # endif /* TS_DIF2 && TS_DIF4 && SPONGE_DIF2 */ # undef T_mld_nnew # undef T_mld_nstp # endif /* DIAGNOSTICS_TS_MLD */ # ifdef DIAGNOSTICS_DEBUG if (istr.eq.1 .and. jstr.eq.1 # ifdef TEMPERATURE & .and. itrc.eq.itemp # endif /* TEMPERATURE */ & ) then i=5 j=5 # if defined DIAGNOSTICS_TS_MLD cff=Trate_mld(i,j,itrc)- & TXadv_mld(i,j,itrc)- & TYadv_mld(i,j,itrc)- & TVadv_mld(i,j,itrc)- & TVmix_mld(i,j,itrc)- & THmix_mld(i,j,itrc)- & Tentr_mld(i,j,itrc)- & TForc_mld(i,j,itrc) MPI_master_only write(stdout,'(A,2x,1pe10.3)') & ' STEP3D_T: T budget closure MLD :', cff ! # endif /* DIAGNOSTICS_TS_MLD */ cff=Trate(i,j,N-5,itrc)- & TXadv(i,j,N-5,itrc)- & TYadv(i,j,N-5,itrc)- & TVadv(i,j,N-5,itrc)- & TVmix(i,j,N-5,itrc)- & THmix(i,j,N-5,itrc)- & TForc(i,j,N-5,itrc) MPI_master_only write(stdout,'(A,2x,1pe10.3)') & ' STEP3D_T: T budget closure at k=N-5:', cff endif # endif /* DIAGNOSTICS_DEBUG */ # endif /* DIAGNOSTICS_TS */ ! ! !====================================================================== ! Set PHYSICAL lateral boundary conditions. !====================================================================== ! !---------------------------------------------------------- ! Apply land mask !---------------------------------------------------------- ! # ifdef MASKING # ifdef EW_PERIODIC # define I_RANGE Istr,Iend # else # define I_RANGE IstrR,IendR # endif # ifdef NS_PERIODIC # define J_RANGE Jstr,Jend # else # define J_RANGE JstrR,JendR # endif do k=1,N do j=J_RANGE do i=I_RANGE t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)*rmask(i,j) enddo enddo enddo # undef I_RANGE # undef J_RANGE # endif /* MASKING */ ! !---------------------------------------------------------- ! Exchange periodic boundaries and computational margins. !---------------------------------------------------------- ! # if defined EW_PERIODIC || defined NS_PERIODIC || defined MPI # ifdef THREE_GHOST_POINTS_TS call exchange_r3d_3pts_tile (Istr,Iend,Jstr,Jend, & t(START_2D_ARRAY,1,nnew,itrc)) # else call exchange_r3d_tile (Istr,Iend,Jstr,Jend, & t(START_2D_ARRAY,1,nnew,itrc)) # endif # endif enddo ! <-- itrc return end # endif /* TRACERS */ ! !====================================================================== ! # if defined BIO_HADV_WENO5 || \ defined TS_HADV_WENO5 || defined TS_VADV_WENO5 || \ defined UV_HADV_WENO5 || defined UV_VADV_WENO5 || \ defined W_HADV_WENO5 || defined W_VADV_WENO5 || \ defined BEDLOAD_WENO5 ! function flux5_weno(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) ! ! WENO5: Weighted Essentially Non-Oscillatory scheme ! with 5th-order accuracy ! ! This function computes tracer reconstruction at the grid cell's left ! edge (u-point i-1/2 or v-point j-1/2). WENO5 uses a convex combination ! of the polynomials reconstructed on the three ENO3 stencils in order ! to achieve higher accuracy on smooth profiles. Both left and right ! combinations of mirror symetric stencils around i-1/2 (j-1/2) points ! are computed. For stability, the upstream stencil combination, ! identified by the sign of ua (u or v), is selected. ! ! The scheme includes improvements from Borges et al., 2008 (WENO_Z), i.e., ! smoothness indicators of higher order (with new non-oscillatory weights) ! that provide a scheme with less dissipation, higher resolution and better ! monotonicity preservation than the classical WENO5. In the meantime, WENO_Z ! removes the need to tune parameter Eps for added dissipation (needed in the ! original scheme due to suboptimal performance on critical points). ! ! References: ! Guang-Shan JIANG and Chi-Wang SHU, 1996: Efficient Implementation of ! Weighted ENO Schemes. JOURNAL OF COMPUTATIONAL PHYSICS 126, 202–228. ! ! Rong WANG and Raymond J. SPITERI, 2007: Linear instability of the ! fifth-order WENO method, SIAM J. Numer. Anal., 45, 1871-1901. ! ! Borges R., M. Carmona, B. Costa, W.S. Don, 2008: An improved weighted ! essentially non-oscillatory scheme for hyperbolic conservation laws. ! Journal of Computational Physics 227 (2008) 3191–3211 ! ! Implementation: P. Marchesiello and J. Demange 2013 ! implicit none REAL :: flux5_weno REAL :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua REAL :: IS0, IS1, IS2 REAl :: d0, d1, d2 REAl :: a0, a1, a2 REAL :: w0, w1, w2 REAL :: p0, p1, p2 REAL :: Eps, cff1, cff2, T5 # define WENO_Z ! Epsilon value for computing weigths # ifdef WENO_Z Eps = 1.e-40 # else ! --> from 1.e-7 (smoother) to 1.e-5 (sharper) Eps = 1.e-6 # endif ! Weigths coefficients d0=1./10. d1=6./10. d2=3./10. ! other coefficients cff1=13./12. cff2=1./6. ! if (ua .ge. 0.) then ! Take upstream stencils for stability ! ! === Reconstruction at u_i-1/2 using LEFT stencils === ! ! Nonlinear Smoothness Indicators IS0 = cff1*(q_im3 - 2.*q_im2 + q_im1)**2 & + 0.25*(q_im3 - 4.*q_im2 + 3*q_im1)**2 IS1 = cff1*(q_im2 - 2.*q_im1 + q_i )**2 & + 0.25*(q_im2 - q_i )**2 IS2 = cff1*(q_im1 - 2.*q_i + q_ip1)**2 & + 0.25*(3.*q_im1 - 4.*q_i + q_ip1)**2 # ifdef WENO_Z ! Non-normalized stencil weights ! with parameter T5 of new smoothness indicator T5 = abs(IS2-IS0) a0 = d0*(1+T5/(Eps+IS0)) a1 = d1*(1+T5/(Eps+IS1)) a2 = d2*(1+T5/(Eps+IS2)) # else a0 = d0/(Eps+IS0)**2 a1 = d1/(Eps+IS1)**2 a2 = d2/(Eps+IS2)**2 # endif ! Normalized Weigths w0 = a0/(a0+a1+a2) w1 = a1/(a0+a1+a2) w2 = a2/(a0+a1+a2) ! Polynomials p0 = cff2*(2.*q_im3 - 7.*q_im2 + 11.*q_im1) p1 = cff2*( -q_im2 + 5.*q_im1 + 2.*q_i ) p2 = cff2*(2.*q_im1 + 5.*q_i - q_ip1) ! Combination flux5_weno = w0*p0 + w1*p1 + w2*p2 else ! ! === Reconstruction at u_i-1/2 using RIGHT stencils === ! ! Nonlinear Smoothness Indicators IS0 = cff1*(q_ip2 - 2.*q_ip1 + q_i )**2 & + 0.25*(q_ip2 - 4.*q_ip1 + 3*q_i )**2 IS1 = cff1*(q_ip1 - 2.*q_i + q_im1)**2 & + 0.25*(q_ip1 - q_im1)**2 IS2 = cff1*(q_i - 2.*q_im1 + q_im2)**2 & + 0.25*(3.*q_i - 4.*q_im1 + q_im2)**2 # ifdef WENO_Z ! Non-normalized stencil weights ! with parameter T5 of new smoothness indicator T5 = abs(IS2-IS0) a0 = d0*(1+T5/(Eps+IS0)) a1 = d1*(1+T5/(Eps+IS1)) a2 = d2*(1+T5/(Eps+IS2)) # else a0 = d0/(Eps+IS0)**2 a1 = d1/(Eps+IS1)**2 a2 = d2/(Eps+IS2)**2 # endif ! Normalized Weigths w0 = a0/(a0+a1+a2) w1 = a1/(a0+a1+a2) w2 = a2/(a0+a1+a2) ! Polynomials p0 = cff2*(2.*q_ip2 - 7.*q_ip1 + 11.*q_i ) p1 = cff2*( -q_ip1 + 5.*q_i + 2.*q_im1) p2 = cff2*(2.*q_i + 5.*q_im1 - q_im2) ! Combination flux5_weno = w0*p0 + w1*p1 + w2*p2 endif return end ! !====================================================================== ! function flux3_weno( q_im2, q_im1, q_i, q_ip1, ua) ! ! WENO3: Weighted Essentially Non-Oscillatory of 3rd-order Accuracy ! ! This function computes tracer reconstruction at the grid cell's left ! edge (u-point i-1/2 or v-point j-1/2). ! ! Implementation: P. Marchesiello and J. Demange 2013 ! implicit none REAL :: flux3_weno REAL :: q_im2, q_im1, q_i, q_ip1, ua REAL :: IS0, IS1 REAl :: a0, a1 REAL :: w0, w1 REAL :: p0, p1 REAL :: Eps, d0,d1, T3 ! ! Epsilon value for computing weigths # ifdef WENO_Z Eps = 1.e-40 # else ! --> from 1.e-7 (smoother) to 1.e-5 (sharper) Eps = 1.e-6 # endif d0=1./3. d1=2./3. ! if (ua .ge. 0.) then ! Take upstream stencils for stability IS0 = (q_im1-q_im2)**2 IS1 = (q_im1-q_i)**2 # ifdef WENO_Z T3 = abs(IS1-IS0) a0 = d0*(1+T3/(Eps+IS0)) a1 = d1*(1+T3/(Eps+IS1)) # else a0 = 1./(3.*(Eps+IS0)**2) a1 = 2./(3.*(Eps+IS1)**2) # endif w0 = a0/(a0+a1) w1 = a1/(a0+a1) p0 = 1./2.*(3.*q_im1-q_im2) p1 = 1./2.*(q_im1+q_i) flux3_weno = w0*p0 + w1*p1 else IS0 = (q_i-q_ip1)**2 IS1 = (q_im1-q_i)**2 # ifdef WENO_Z T3 = abs(IS1-IS0) a0 = d0*(1+T3/(Eps+IS0)) a1 = d1*(1+T3/(Eps+IS1)) # else a0 = 1./(3.*(Eps+IS0)**2) a1 = 2./(3.*(Eps+IS1)**2) # endif w0 = a0/(a0+a1) w1 = a1/(a0+a1) p0 = 1./2.*(3.*q_i-q_ip1) p1 = 1./2.*(q_im1+q_i) flux3_weno = w0*p0 + w1*p1 endif return end # endif /* WENO */ ! !====================================================================== ! #else subroutine step3d_t_empty return end #endif /* SOLVE3D */