MODULE dynldf_iso !!====================================================================== !! *** MODULE dynldf_iso *** !! Ocean dynamics: lateral viscosity trend !!====================================================================== !! History : OPA ! 97-07 (G. Madec) Original code !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module !! - ! 2004-08 (C. Talandier) New trends organization !! 2.0 ! 2005-11 (G. Madec) s-coordinate: horizontal diffusion !!---------------------------------------------------------------------- #if defined key_ldfslp || defined key_esopa !!---------------------------------------------------------------------- !! 'key_ldfslp' slopes of the direction of mixing !!---------------------------------------------------------------------- !! dyn_ldf_iso : update the momentum trend with the horizontal part !! of the lateral diffusion using isopycnal or horizon- !! tal s-coordinate laplacian operator. !!---------------------------------------------------------------------- USE oce ! ocean dynamics and tracers USE dom_oce ! ocean space and time domain USE ldfdyn_oce ! ocean dynamics lateral physics USE ldftra_oce ! ocean tracer lateral physics USE zdf_oce ! ocean vertical physics USE trdmod ! ocean dynamics trends USE trdmod_oce ! ocean variables trends USE ldfslp ! iso-neutral slopes USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE in_out_manager ! I/O manager USE lib_mpp ! MPP library USE prtctl ! Print control IMPLICIT NONE PRIVATE PUBLIC dyn_ldf_iso ! called by step.F90 PUBLIC dyn_ldf_iso_alloc ! called by nemogcm.F90 !FTRANS zdiu zdju zdiv zdjv zdj1u zdj1v zfuw zfvw :I :z REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u ! 2D workspace (dyn_ldf_iso) REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfvw, zdiv, zdjv, zdj1v ! - - !! * Control permutation of array indices # include "oce_ftrans.h90" # include "dom_oce_ftrans.h90" # include "ldfdyn_oce_ftrans.h90" # include "ldftra_oce_ftrans.h90" # include "ldfslp_ftrans.h90" # include "zdf_oce_ftrans.h90" !! * Substitutions # include "domzgr_substitute.h90" # include "ldfdyn_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OPA 3.3 , NEMO Consortium (2011) !! $Id$ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS INTEGER FUNCTION dyn_ldf_iso_alloc() !!---------------------------------------------------------------------- !! *** ROUTINE dyn_ldf_iso_alloc *** !!---------------------------------------------------------------------- ALLOCATE( zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) , & & zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc ) ! IF( dyn_ldf_iso_alloc /= 0 ) CALL ctl_warn('dyn_ldf_iso_alloc: array allocate failed.') END FUNCTION dyn_ldf_iso_alloc SUBROUTINE dyn_ldf_iso( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE dyn_ldf_iso *** !! !! ** Purpose : Compute the before trend of the rotated laplacian !! operator of lateral momentum diffusion except the diagonal !! vertical term that will be computed in dynzdf module. Add it !! to the general trend of momentum equation. !! !! ** Method : !! The momentum lateral diffusive trend is provided by a 2nd !! order operator rotated along neutral or geopotential surfaces !! (in s-coordinates). !! It is computed using before fields (forward in time) and isopyc- !! nal or geopotential slopes computed in routine ldfslp. !! Here, u and v components are considered as 2 independent scalar !! fields. Therefore, the property of splitting divergent and rota- !! tional part of the flow of the standard, z-coordinate laplacian !! momentum diffusion is lost. !! horizontal fluxes associated with the rotated lateral mixing: !! u-component: !! ziut = ( ahmt + ahmb0 ) e2t * e3t / e1t di[ ub ] !! - ahmt e2t * mi-1(uslp) dk[ mi(mk(ub)) ] !! zjuf = ( ahmf + ahmb0 ) e1f * e3f / e2f dj[ ub ] !! - ahmf e1f * mi(vslp) dk[ mj(mk(ub)) ] !! v-component: !! zivf = ( ahmf + ahmb0 ) e2t * e3t / e1t di[ vb ] !! - ahmf e2t * mj(uslp) dk[ mi(mk(vb)) ] !! zjvt = ( ahmt + ahmb0 ) e1f * e3f / e2f dj[ ub ] !! - ahmt e1f * mj-1(vslp) dk[ mj(mk(vb)) ] !! take the horizontal divergence of the fluxes: !! diffu = 1/(e1u*e2u*e3u) { di [ ziut ] + dj-1[ zjuf ] } !! diffv = 1/(e1v*e2v*e3v) { di-1[ zivf ] + dj [ zjvt ] } !! Add this trend to the general trend (ua,va): !! ua = ua + diffu !! CAUTION: here the isopycnal part is with a coeff. of aht. This !! should be modified for applications others than orca_r2 (!!bug) !! !! ** Action : !! Update (ua,va) arrays with the before geopotential biharmonic !! mixing trend. !! Update (avmu,avmv) to accompt for the diagonal vertical component !! of the rotated operator in dynzdf module !!---------------------------------------------------------------------- USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released #if ! defined key_z_first ! Then these workspace arrays can be left as 2D USE wrk_nemo, ONLY: zjvt => wrk_2d_3 ! 2D workspace USE wrk_nemo, ONLY: zivf => wrk_2d_4 ! 2D workspace USE wrk_nemo, ONLY: ziut => wrk_2d_1 , zjuf => wrk_2d_2 USE wrk_nemo, ONLY: zdku => wrk_2d_5 , zdkv => wrk_2d_6 USE wrk_nemo, ONLY: zdk1u => wrk_2d_7 , zdk1v => wrk_2d_8 #endif USE timing, ONLY: timing_start, timing_stop ! INTEGER, INTENT( in ) :: kt ! ocean time-step index ! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: zabe1, zabe2, zcof1, zcof2 ! local scalars REAL(wp) :: zmskt, zmskf, zbu, zbv, zuah, zvah ! - - REAL(wp) :: zcoef0, zcoef3, zcoef4, zmkt, zmkf ! - - REAL(wp) :: zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj ! - - REAL(wp) :: zcof, zrecip #if defined key_z_first REAL(wp) :: zdku, zdk1u, zdki1u, zdk1i1u, zdkim1u, zdk1im1u REAL(wp) :: zdkj1u, zdk1j1u, zdkjm1u, zdk1jm1u REAL(wp) :: zdkv, zdk1v, zdki1v, zdk1i1v, zdkim1v, zdk1im1v REAL(wp) :: zdkj1v, zdk1j1v, zdkjm1v, zdk1jm1v REAL(wp) :: aziut, aziuti1, azjuf, azjufjm1 REAL(wp) :: azivf, azivfim1, azjvtj1, azjvt #endif !!---------------------------------------------------------------------- ! ziut zjvt zjuf zivf :I :I :z !!$#if defined key_z_first !!$ REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ziut, zjuf, zivf, zjvt !!$ ALLOCATE( ziut(jpi,jpj,jpk), zjuf(jpi,jpj,jpk), & !!$ zivf(jpi,jpj,jpk), zjvt(jpi,jpj,jpk) ) !!$#endif CALL timing_start('dyn_ldf_iso') IF( wrk_in_use(2, 1,2,3,4,5,6,7,8) ) THEN CALL ctl_stop('dyn_ldf_iso: requested workspace arrays unavailable') ; RETURN END IF IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or ' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate horizontal diffusive operator' ! ! allocate dyn_ldf_bilap arrays IF( dyn_ldf_iso_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays') ENDIF ! s-coordinate: Iso-level diffusion on momentum but not on tracer IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN ! #if defined key_z_first DO jj = 2, jpjm1 ! set the slopes of iso-level DO ji = fs_2, fs_jpim1 DO jk = 1, mbkmax(ji,jj) ! jpk #else DO jk = 1, jpk ! set the slopes of iso-level DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. #endif uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk) vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk) wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5 wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5 END DO END DO END DO ! Lateral boundary conditions on the slopes CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) !!bug IF( kt == nit000 ) then IF(lwp) WRITE(numout,*) ' max slop: u', SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)), & & ' wi', sqrt(MAXVAL(wslpi)) , ' wj', sqrt(MAXVAL(wslpj)) endif !!end ENDIF !CALL timing_start('dyn_ldf_iso_hslab') #if defined key_z_first ! Vertical u- and v-shears at level jk and jk+1 ! --------------------------------------------- ! surface boundary condition: zdku(jk=1)=zdku(jk=2) ! zdkv(jk=1)=zdkv(jk=2) !!$ DO jj = 1, jpj, 1 !!$ DO ji = 1, jpi, 1 !!$ !!$ ! jk=1 special case !!$ !zdk1u(ji,jj,1) = ( ub(ji,jj,1) -ub(ji,jj,2) ) * umask(ji,jj,2) !!$ zdk1v(ji,jj,1) = ( vb(ji,jj,1) -vb(ji,jj,2) ) * vmask(ji,jj,2) !!$ !zdku(ji,jj,1) = zdk1u(ji,jj,1) !!$ zdkv(ji,jj,1) = zdk1v(ji,jj,1) !!$ !!$ DO jk = 2, jpkm1, 1 !!$ !zdk1u(ji,jj,jk) = ( ub(ji,jj,jk) -ub(ji,jj,jk+1) ) * umask(ji,jj,jk+1) !!$ zdk1v(ji,jj,jk) = ( vb(ji,jj,jk) -vb(ji,jj,jk+1) ) * vmask(ji,jj,jk+1) !!$ !!$ !zdku(ji,jj,jk) = ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) * umask(ji,jj,jk) !!$ zdkv(ji,jj,jk) = ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) * vmask(ji,jj,jk) !!$ END DO !jk !!$ END DO !!$ END DO !!$ ! -----f----- !!$ ! Horizontal fluxes on U | !!$ ! --------------------=== t u t !!$ ! | !!$ ! i-flux at t-point -----f----- !!$ !!$ IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u) !!$ DO jj = 2, jpjm1 !!$ DO ji = 2, jpi !!$ DO jk = 1, jpkm1, 1 !!$ zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji,jj) !!$ !!$ zmskt = 1._wp/MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & !!$ & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1._wp ) !!$ !!$ zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5_wp * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) !!$ !!$ ziut(ji,jj,jk) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) & !!$ & + zcof1 * ( zdku + zdk1im1u & !!$ & +zdk1u + zdkim1u ) ) * tmask(ji,jj,jk) !!$ END DO !!$ END DO !!$ END DO !!$ ELSE ! other coordinate system (zco or sco) : e3t !!$ DO jj = 2, jpjm1 !!$ DO ji = 2, jpi !!$ DO jk = 1, jpkm1, 1 !!$ zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj) !!$ !!$ zmskt = 1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & !!$ & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. ) !!$ !!$ zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) !!$ !!$ ziut(ji,jj,jk) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) & !!$ & + zcof1 * ( zdku(ji,jj,jk) + zdk1u(ji-1,jj,jk) & !!$ & +zdk1u(ji,jj,jk) + zdku (ji-1,jj,jk) ) ) * tmask(ji,jj,jk) !!$ !!$ END DO !!$ END DO !!$ END DO !!$ ENDIF ! j-flux at f-point ! BLOCKABLE(ji,jj,jk) ! BLOCKING SIZE (0) !!$ DO jj = 1, jpjm1, 1 !!$! BLOCKING SIZE (0) !!$ DO ji = 1, jpim1, 1 !!$! BLOCKING SIZE (4) !!$ DO jk = 1, jpkm1, 1 !!$ zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) !!$ !!$ zmskf = 1./MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) & !!$ & + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1. ) !!$ !!$ zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) !!$ !!$ zjuf(ji,jj,jk) = ( zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) & !!$ & + zcof2 * ( zdku (ji,jj+1,jk) + zdk1u(ji,jj,jk) & !!$ & +zdk1u(ji,jj+1,jk) + zdku (ji,jj,jk) ) ) * fmask(ji,jj,jk) !!$ END DO !!$ END DO !!$ END DO ! | t | ! Horizontal fluxes on V | | ! --------------------=== f---v---f ! | | ! | t | !!$ IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u) !!$ DO jj = 2, jpj !!$ DO ji = 1, jpim1 !!$ DO jk = 1, jpkm1, 1 !!$ !!$ ! i-flux at f-point | t | !!$ zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj) !!$ !!$ zmskf = 1._wp/MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & !!$ + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1._wp ) !!$ !!$ zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5_wp * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) !!$ !!$ zivf(ji,jj,jk) = ( zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) & !!$ + zcof1 * ( zdkv(ji,jj,jk) + zdk1v(ji+1,jj,jk) & !!$ +zdk1v(ji,jj,jk) + zdkv(ji+1,jj,jk) ) ) * fmask(ji,jj,jk) !!$ !!$ !!$ ! j-flux at t-point !!$ zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj) !!$ !!$ zmskt = 1._wp/MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & !!$ + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1._wp ) !!$ !!$ zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5_wp * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) !!$ !!$ zjvt(ji,jj,jk) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) & !!$ + zcof2 * ( zdkv(ji,jj-1,jk) + zdk1v(ji,jj,jk) & !!$ +zdk1v(ji,jj-1,jk) + zdkv (ji,jj,jk) ) ) * tmask(ji,jj,jk) !!$ END DO !!$ END DO !!$ END DO !!$ ELSE ! other coordinate system (zco or sco) : e3t !!$ DO jj = 2, jpj !!$ DO ji = 1, jpim1 !!$ DO jk = 1, jpkm1 !!$ !!$ ! i-flux at f-point !!$ zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj) !!$ !!$ zmskf = 1._wp/MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & !!$ + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1._wp ) !!$ !!$ zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5_wp * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) !!$ !!$ zivf(ji,jj,jk) = ( zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) & !!$ + zcof1 * ( zdkv(ji,jj,jk) + zdk1v(ji+1,jj,jk) & !!$ +zdk1v(ji,jj,jk) + zdkv(ji+1,jj,jk) ) ) * fmask(ji,jj,jk) !!$ ! j-flux at t-point !!$ zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj) !!$ !!$ zmskt = 1._wp/MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & !!$ + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1._wp ) !!$ !!$ zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5_wp * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) !!$ !!$ zjvt(ji,jj,jk) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) & !!$ + zcof2 * ( zdkv (ji,jj-1,jk) + zdk1v(ji,jj,jk) & !!$ +zdk1v(ji,jj-1,jk) + zdkv (ji,jj,jk) ) ) * tmask(ji,jj,jk) !!$ END DO !!$ END DO !!$ END DO !!$ !!$ ENDIF IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u) DO jj = 2, jpjm1 DO ji = 2, jpim1 !DIR$ SHORTLOOP ! SAFE_ADDRESS allows the compiler to generate code that will vectorise despite the ! If condition on jk>1. This does mean that there will be out-of-bounds *reads* ! when jk is 1 but that doesn't matter. !DIR$ SAFE_ADDRESS DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1, 1 ! Vertical u- and v-shears at level jk and jk+1 ! --------------------------------------------- ! surface boundary condition: zdku(jk=1)=zdku(jk=2) ! zdkv(jk=1)=zdkv(jk=2) zdk1u = ( ub(ji ,jj ,jk) -ub(ji ,jj ,jk+1) ) * umask(ji ,jj,jk+1) zdk1i1u = ( ub(ji+1,jj ,jk) -ub(ji+1,jj ,jk+1) ) * umask(ji+1,jj,jk+1) zdk1j1u = ( ub(ji ,jj+1,jk) -ub(ji ,jj+1,jk+1) ) * umask(ji ,jj+1,jk+1) zdk1im1u = ( ub(ji-1,jj ,jk) -ub(ji-1,jj ,jk+1) ) * umask(ji-1,jj,jk+1) zdk1jm1u = ( ub(ji ,jj-1,jk) -ub(ji ,jj-1,jk+1) ) * umask(ji ,jj-1,jk+1) IF(jk > 1)THEN zdku = ( ub(ji ,jj ,jk-1) - ub(ji ,jj ,jk) ) * umask(ji,jj,jk) zdki1u = ( ub(ji+1,jj ,jk-1) - ub(ji+1,jj ,jk) ) * umask(ji+1,jj,jk) zdkim1u= ( ub(ji-1,jj ,jk-1) - ub(ji-1,jj ,jk) ) * umask(ji-1,jj,jk) zdkj1u = ( ub(ji ,jj+1,jk-1) - ub(ji ,jj+1,jk) ) * umask(ji,jj+1,jk) zdkjm1u= ( ub(ji ,jj-1,jk-1) - ub(ji ,jj-1,jk) ) * umask(ji,jj-1,jk) ELSE zdku = zdk1u zdki1u = zdk1i1u zdkim1u= zdk1im1u zdkj1u = zdk1j1u zdkjm1u= zdk1jm1u END IF zdku = zdku + zdk1u zdki1u = zdki1u + zdk1i1u zdkim1u = zdkim1u + zdk1im1u zdkj1u = zdkj1u +zdk1j1u zdkjm1u = zdk1jm1u + zdkjm1u ! volume elements zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ! horizontal component of isopycnal momentum diffusive trends ! -----f----- ! Horizontal fluxes on U | ! --------------------=== t u t ! | ! i-flux at t-point -----f----- ! z-coordinate - partial steps : min(e3u) aziuti1 = ( ((fsahmt(ji+1,jj,jk)+ahmb0) * e2t(ji+1,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji+1,jj)) * & ( ub(ji+1,jj,jk) - ub(ji,jj,jk) ) + (- aht0 * e2t(ji+1,jj) * & (1._wp/MAX( umask(ji,jj,jk )+umask(ji+1,jj,jk+1) & + umask(ji,jj,jk+1)+umask(ji+1,jj,jk ), 1._wp )) * 0.5 * & ( uslp(ji,jj,jk) + uslp(ji+1,jj,jk) )) * & ( zdku + zdki1u ) ) * tmask(ji+1,jj,jk) aziut = ( ((fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * & ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) + (- aht0 * e2t(ji,jj) * & (1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. )) * 0.5 * & ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )) * & ( zdku + zdkim1u ) ) * tmask(ji,jj,jk) ! j-flux at f-point - identical to non-ln_zps azjuf = ( (( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)) * & ( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) + (- aht0 * e1f(ji,jj) * & (1./MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) & + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1. )) * 0.5 * & ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )) * & ( zdku + zdkj1u ) ) * fmask(ji,jj,jk) azjufjm1 = ( (( fsahmf(ji,jj-1,jk) + ahmb0 ) * e1f(ji,jj-1) * fse3f(ji,jj-1,jk) / e2f(ji,jj-1)) * & ( ub(ji,jj,jk) - ub(ji,jj-1,jk) ) + (- aht0 * e1f(ji,jj-1) * & (1./MAX( umask(ji,jj,jk )+umask(ji,jj-1,jk+1) & + umask(ji,jj,jk+1)+umask(ji,jj-1,jk ), 1. )) * 0.5 * & ( vslp(ji+1,jj-1,jk) + vslp(ji,jj-1,jk) )) * & ( zdku + zdkjm1u ) ) * fmask(ji,jj-1,jk) ! Second derivative (divergence) and add to the general trend ! ----------------------------------------------------------- zuah =( & ! ziut (ji+1,jj,jk) - & aziuti1 - & ! ziut (ji ,jj,jk) + & aziut + & ! zjuf (ji ,jj,jk) - & azjuf - & ! zjuf (ji,jj-1,jk) & azjufjm1 & ) / zbu ! add the trends to the general trends ua (ji,jj,jk) = ua (ji,jj,jk) + zuah END DO !DIR$ SHORTLOOP ! SAFE_ADDRESS allows the compiler to generate code that will vectorise despite the ! If condition on jk>1. This does mean that there will be out-of-bounds *reads* ! when jk is 1 but that doesn't matter. !DIR$ SAFE_ADDRESS DO jk = 1, mbkmax(ji,jj)-1, 1 ! jpkm1, 1 zdk1v = ( vb(ji ,jj ,jk) -vb(ji ,jj ,jk+1) ) * vmask(ji ,jj,jk+1) zdk1i1v = ( vb(ji+1,jj ,jk) -vb(ji+1,jj ,jk+1) ) * vmask(ji+1,jj,jk+1) zdk1im1v = ( vb(ji-1,jj ,jk) -vb(ji-1,jj ,jk+1) ) * vmask(ji-1,jj,jk+1) zdk1j1v = ( vb(ji ,jj+1,jk) -vb(ji ,jj+1,jk+1) ) * vmask(ji ,jj+1,jk+1) zdk1jm1v = ( vb(ji ,jj-1,jk) -vb(ji ,jj-1,jk+1) ) * vmask(ji ,jj-1,jk+1) IF(jk > 1)THEN zdkv = ( vb(ji ,jj ,jk-1) - vb(ji ,jj ,jk) ) * vmask(ji,jj,jk) zdki1v = ( vb(ji+1,jj ,jk-1) - vb(ji+1,jj ,jk) ) * vmask(ji+1,jj,jk) zdkim1v= ( vb(ji-1,jj ,jk-1) - vb(ji-1,jj ,jk) ) * vmask(ji-1,jj,jk) zdkj1v = ( vb(ji ,jj+1,jk-1) - vb(ji ,jj+1,jk) ) * vmask(ji,jj+1,jk) zdkjm1v= ( vb(ji ,jj-1,jk-1) - vb(ji ,jj-1,jk) ) * vmask(ji,jj-1,jk) ELSE zdkv = zdk1v zdki1v = zdk1i1v zdkim1v= zdk1im1v zdkj1v = zdk1j1v zdkjm1v= zdk1jm1v END IF zdkv = zdkv + zdk1v zdki1v = zdk1i1v + zdki1v zdkim1v = zdkim1v +zdk1im1v zdkj1v = zdk1j1v + zdkj1v zdkjm1v = zdkjm1v +zdk1jm1v ! | t | ! Horizontal fluxes on V | | ! --------------------=== f---v---f ! | | ! | t | ! i-flux at f-point - identical to non-ln_zps case azivf = ( (( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / & e1f(ji,jj)) * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) & + (- aht0 * e2f(ji,jj) / & MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1._wp ) * & 0.5_wp * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * & ( zdkv + zdki1v ) ) * fmask(ji,jj,jk) azivfim1 = ( (( fsahmf(ji-1,jj,jk) + ahmb0 ) * e2f(ji-1,jj) * & fse3f(ji-1,jj,jk) / e1f(ji-1,jj)) * ( vb(ji,jj,jk) - vb(ji-1,jj,jk) ) & + (- aht0 * e2f(ji-1,jj) / & MAX( vmask(ji,jj,jk )+vmask(ji-1,jj,jk+1) & + vmask(ji,jj,jk+1)+vmask(ji-1,jj,jk ), 1._wp ) * 0.5_wp * & ( uslp(ji-1,jj+1,jk) + uslp(ji-1,jj,jk) )) * & ( zdkim1v + zdkv ) ) * fmask(ji-1,jj,jk) ! j-flux at t-point - min(e3u) instead of e3t azjvtj1 = ( & ((fsahmt(ji,jj+1,jk)+ahmb0) * e1t(ji,jj+1) * MIN( fse3v(ji,jj+1,jk), fse3v(ji,jj,jk) ) / & e2t(ji,jj+1)) * ( vb(ji,jj+1,jk) - vb(ji,jj,jk) ) & + (- aht0 * e1t(ji,jj+1) / & MAX( vmask(ji,jj,jk )+vmask(ji,jj+1,jk+1) & + vmask(ji,jj,jk+1)+vmask(ji,jj+1,jk ), 1._wp ) * & 0.5_wp * ( vslp(ji,jj,jk) + vslp(ji,jj+1,jk) )) * & ( zdkv + zdkj1v ) & ) * tmask(ji,jj+1,jk) azjvt = ( & ((fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / & e2t(ji,jj)) * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) & + (- aht0 * e1t(ji,jj) /MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1._wp ) * 0.5_wp * & ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * & ( zdkjm1v + zdkv ) ) * tmask(ji,jj,jk) ! volume elements zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ! Second derivative (divergence) and add to the general trend ! ----------------------------------------------------------- zvah =( & ! zivf (ji ,jj ,jk) - & azivf - & ! zivf (ji-1,jj ,jk) + & azivfim1 + & ! zjvt (ji ,jj+1,jk) - & azjvtj1 - & ! zjvt (ji ,jj ,jk) & azjvt & ) / zbv ! add the trend to the general trend va (ji,jj,jk) = va (ji,jj,jk) + zvah END DO END DO END DO ELSE ! other coordinate system (zco or sco) : e3t DO jj = 2, jpjm1 DO ji = 2, jpim1 !DIR$ SHORTLOOP ! SAFE_ADDRESS allows the compiler to generate code that will vectorise despite the ! If condition on jk>1. This does mean that there will be out-of-bounds *reads* ! when jk is 1 but that doesn't matter. !DIR$ SAFE_ADDRESS DO jk = 1, mbkmax(ji,jj)-1, 1 ! jpkm1, 1 ! Vertical u- and v-shears at level jk and jk+1 ! --------------------------------------------- ! surface boundary condition: zdku(jk=1)=zdku(jk=2) ! zdkv(jk=1)=zdkv(jk=2) zdk1u = ( ub(ji ,jj ,jk) -ub(ji ,jj ,jk+1) ) * umask(ji ,jj,jk+1) zdk1i1u = ( ub(ji+1,jj ,jk) -ub(ji+1,jj ,jk+1) ) * umask(ji+1,jj,jk+1) zdk1j1u = ( ub(ji ,jj+1,jk) -ub(ji ,jj+1,jk+1) ) * umask(ji ,jj+1,jk+1) zdk1im1u = ( ub(ji-1,jj ,jk) -ub(ji-1,jj ,jk+1) ) * umask(ji-1,jj,jk+1) zdk1jm1u = ( ub(ji ,jj-1,jk) -ub(ji ,jj-1,jk+1) ) * umask(ji ,jj-1,jk+1) IF(jk > 1)THEN zdku = ( ub(ji ,jj ,jk-1) - ub(ji ,jj ,jk) ) * umask(ji,jj,jk) zdki1u = ( ub(ji+1,jj ,jk-1) - ub(ji+1,jj ,jk) ) * umask(ji+1,jj,jk) zdkim1u= ( ub(ji-1,jj ,jk-1) - ub(ji-1,jj ,jk) ) * umask(ji-1,jj,jk) zdkj1u = ( ub(ji ,jj+1,jk-1) - ub(ji ,jj+1,jk) ) * umask(ji,jj+1,jk) zdkjm1u= ( ub(ji ,jj-1,jk-1) - ub(ji ,jj-1,jk) ) * umask(ji,jj-1,jk) ELSE zdku = zdk1u zdki1u = zdk1i1u zdkim1u= zdk1im1u zdkj1u = zdk1j1u zdkjm1u= zdk1jm1u END IF zdku = zdku + zdk1u zdki1u = zdki1u + zdk1i1u zdkim1u = zdkim1u + zdk1im1u zdkj1u = zdkj1u +zdk1j1u zdkjm1u = zdk1jm1u + zdkjm1u ! volume element zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ! horizontal component of isopycnal momentum diffusive trends ! -----f----- ! Horizontal fluxes on U | ! --------------------=== t u t ! | ! i-flux at t-point -----f----- ! other coordinate system (zco or sco) : e3t aziuti1 = ( ((fsahmt(ji+1,jj,jk)+ahmb0) * e2t(ji+1,jj) * fse3t(ji+1,jj,jk) / e1t(ji+1,jj)) * & ( ub(ji+1,jj,jk) - ub(ji,jj,jk) ) + (- aht0 * e2t(ji+1,jj) * & (1./MAX( umask(ji,jj,jk )+umask(ji+1,jj,jk+1) & + umask(ji,jj,jk+1)+umask(ji+1,jj,jk ), 1. )) * 0.5 * & ( uslp(ji,jj,jk) + uslp(ji+1,jj,jk) )) * & ( zdku + zdki1u ) ) * tmask(ji+1,jj,jk) aziut = ( ((fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * & ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) + (- aht0 * e2t(ji,jj) * & (1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. )) * 0.5 * & ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )) * & ( zdku + zdkim1u ) ) * tmask(ji,jj,jk) azjuf = ( (( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)) * & ( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) + (- aht0 * e1f(ji,jj) * & (1./MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) & + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1. )) * 0.5 * & ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )) * & ( zdku + zdkj1u ) ) * fmask(ji,jj,jk) azjufjm1 = ( (( fsahmf(ji,jj-1,jk) + ahmb0 ) * e1f(ji,jj-1) * fse3f(ji,jj-1,jk) / e2f(ji,jj-1)) * & ( ub(ji,jj,jk) - ub(ji,jj-1,jk) ) + (- aht0 * e1f(ji,jj-1) * & (1./MAX( umask(ji,jj,jk )+umask(ji,jj-1,jk+1) & + umask(ji,jj,jk+1)+umask(ji,jj-1,jk ), 1. )) * 0.5 * & ( vslp(ji+1,jj-1,jk) + vslp(ji,jj-1,jk) )) * & ( zdku + zdkjm1u ) ) * fmask(ji,jj-1,jk) ! Second derivative (divergence) and add to the general trend ! ----------------------------------------------------------- zuah =( & ! ziut (ji+1,jj,jk) - & aziuti1 - & ! ziut (ji ,jj,jk) + & aziut + & ! zjuf (ji ,jj,jk) - & azjuf - & ! zjuf (ji,jj-1,jk) & azjufjm1 & ) / zbu ! add the trend to the general trend ua (ji,jj,jk) = ua (ji,jj,jk) + zuah END DO !DIR$ SHORTLOOP ! SAFE_ADDRESS allows the compiler to generate code that will vectorise despite the ! If condition on jk>1. This does mean that there will be out-of-bounds *reads* ! when jk is 1 but that doesn't matter. !DIR$ SAFE_ADDRESS DO jk = 1, mbkmax(ji,jj)-1, 1 ! jpkm1, 1 zdk1v = ( vb(ji ,jj ,jk) -vb(ji ,jj ,jk+1) ) * vmask(ji ,jj,jk+1) zdk1i1v = ( vb(ji+1,jj ,jk) -vb(ji+1,jj ,jk+1) ) * vmask(ji+1,jj,jk+1) zdk1im1v = ( vb(ji-1,jj ,jk) -vb(ji-1,jj ,jk+1) ) * vmask(ji-1,jj,jk+1) zdk1j1v = ( vb(ji ,jj+1,jk) -vb(ji ,jj+1,jk+1) ) * vmask(ji ,jj+1,jk+1) zdk1jm1v = ( vb(ji ,jj-1,jk) -vb(ji ,jj-1,jk+1) ) * vmask(ji ,jj-1,jk+1) IF(jk > 1)THEN zdkv = ( vb(ji ,jj ,jk-1) - vb(ji ,jj ,jk) ) * vmask(ji,jj,jk) zdki1v = ( vb(ji+1,jj ,jk-1) - vb(ji+1,jj ,jk) ) * vmask(ji+1,jj,jk) zdkim1v= ( vb(ji-1,jj ,jk-1) - vb(ji-1,jj ,jk) ) * vmask(ji-1,jj,jk) zdkj1v = ( vb(ji ,jj+1,jk-1) - vb(ji ,jj+1,jk) ) * vmask(ji,jj+1,jk) zdkjm1v= ( vb(ji ,jj-1,jk-1) - vb(ji ,jj-1,jk) ) * vmask(ji,jj-1,jk) ELSE zdkv = zdk1v zdki1v = zdk1i1v zdkim1v= zdk1im1v zdkj1v = zdk1j1v zdkjm1v= zdk1jm1v END IF zdkv = zdkv + zdk1v zdki1v = zdk1i1v + zdki1v zdkim1v = zdkim1v +zdk1im1v zdkj1v = zdk1j1v + zdkj1v zdkjm1v = zdkjm1v +zdk1jm1v azivf = ( (( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / & e1f(ji,jj)) * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) & + (- aht0 * e2f(ji,jj) / & MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1._wp ) * & 0.5_wp * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * & ( zdkv + zdki1v ) ) * fmask(ji,jj,jk) azivfim1 = ( (( fsahmf(ji-1,jj,jk) + ahmb0 ) * e2f(ji-1,jj) * & fse3f(ji-1,jj,jk) / e1f(ji-1,jj)) * ( vb(ji,jj,jk) - vb(ji-1,jj,jk) ) & + (- aht0 * e2f(ji-1,jj) / & MAX( vmask(ji,jj,jk )+vmask(ji-1,jj,jk+1) & + vmask(ji,jj,jk+1)+vmask(ji-1,jj,jk ), 1._wp ) * 0.5_wp * & ( uslp(ji-1,jj+1,jk) + uslp(ji-1,jj,jk) )) * & ( zdkim1v + zdkv ) ) * fmask(ji-1,jj,jk) azjvtj1 = ( & ((fsahmt(ji,jj+1,jk)+ahmb0) * e1t(ji,jj+1) * fse3t(ji,jj+1,jk) / & e2t(ji,jj+1)) * ( vb(ji,jj+1,jk) - vb(ji,jj,jk) ) & + (- aht0 * e1t(ji,jj+1) / & MAX( vmask(ji,jj,jk )+vmask(ji,jj+1,jk+1) & + vmask(ji,jj,jk+1)+vmask(ji,jj+1,jk ), 1._wp ) * & 0.5_wp * ( vslp(ji,jj,jk) + vslp(ji,jj+1,jk) )) * & ( zdkv + zdkj1v ) & ) * tmask(ji,jj+1,jk) azjvt = ( ((fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)) * & ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) & + (- aht0 * e1t(ji,jj) /MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1._wp ) * 0.5_wp * & ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * & ( zdkjm1v + zdkv ) ) * tmask(ji,jj,jk) ! volume elements zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ! Second derivative (divergence) and add to the general trend ! ----------------------------------------------------------- zvah =( & ! zivf (ji ,jj ,jk) - & azivf - & ! zivf (ji-1,jj ,jk) + & azivfim1 + & ! zjvt (ji ,jj+1,jk) - & azjvtj1 - & ! zjvt (ji ,jj ,jk) & azjvt & ) / zbv ! add the trend to the general trend va (ji,jj,jk) = va (ji,jj,jk) + zvah END DO END DO END DO END IF ! ln_zps #else ! ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== ! Vertical u- and v-shears at level jk and jk+1 ! --------------------------------------------- ! surface boundary condition: zdku(jk=1)=zdku(jk=2) ! zdkv(jk=1)=zdkv(jk=2) zdk1u(:,:) = ( ub(:,:,jk) -ub(:,:,jk+1) ) * umask(:,:,jk+1) zdk1v(:,:) = ( vb(:,:,jk) -vb(:,:,jk+1) ) * vmask(:,:,jk+1) IF( jk == 1 ) THEN zdku(:,:) = zdk1u(:,:) zdkv(:,:) = zdk1v(:,:) ELSE zdku(:,:) = ( ub(:,:,jk-1) - ub(:,:,jk) ) * umask(:,:,jk) zdkv(:,:) = ( vb(:,:,jk-1) - vb(:,:,jk) ) * vmask(:,:,jk) ENDIF ! -----f----- ! Horizontal fluxes on U | ! --------------------=== t u t ! | ! i-flux at t-point -----f----- IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u) DO jj = 2, jpjm1 DO ji = fs_2, jpi ! vector opt. zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji,jj) zmskt = 1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. ) zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) ziut(ji,jj) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) & & + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & & +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk) END DO END DO ELSE ! other coordinate system (zco or sco) : e3t DO jj = 2, jpjm1 DO ji = fs_2, jpi ! vector opt. zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj) zmskt = 1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. ) zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) ziut(ji,jj) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) & & + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & & +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk) END DO END DO ENDIF ! j-flux at f-point DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) zmskf = 1./MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) & & + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1. ) zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) zjuf(ji,jj) = ( zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) & & + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) & & +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) * fmask(ji,jj,jk) END DO END DO ! | t | ! Horizontal fluxes on V | | ! --------------------=== f---v---f ! | | ! i-flux at f-point | t | DO jj = 2, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj) zmskf = 1./MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & & + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1. ) zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) zivf(ji,jj) = ( zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) & & + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj) & & +zdk1v(ji,jj) + zdkv (ji+1,jj) ) ) * fmask(ji,jj,jk) END DO END DO ! j-flux at t-point IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u) DO jj = 2, jpj DO ji = 1, fs_jpim1 ! vector opt. zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj) zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & & + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. ) zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) zjvt(ji,jj) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) & & + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) & & +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk) END DO END DO ELSE ! other coordinate system (zco or sco) : e3t DO jj = 2, jpj DO ji = 1, fs_jpim1 ! vector opt. zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj) zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & & + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. ) zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) zjvt(ji,jj) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) & & + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) & & +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk) END DO END DO ENDIF ! Second derivative (divergence) and add to the general trend ! ----------------------------------------------------------- DO jj = 2, jpjm1 DO ji = 2, jpim1 !! Question vectop possible??? !!bug ! volume elements zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ! horizontal component of isopycnal momentum diffusive trends zuah =( ziut (ji+1,jj) - ziut (ji,jj ) + & & zjuf (ji ,jj) - zjuf (ji,jj-1) ) / zbu zvah =( zivf (ji,jj ) - zivf (ji-1,jj) + & & zjvt (ji,jj+1) - zjvt (ji,jj ) ) / zbv ! add the trends to the general trends ua (ji,jj,jk) = ua (ji,jj,jk) + zuah va (ji,jj,jk) = va (ji,jj,jk) + zvah END DO END DO ! ! =============== END DO ! End of slab ! ! =============== #endif !CALL timing_stop('dyn_ldf_iso_hslab','section') ! print sum trends (used for debugging) IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' ldfh - Ua: ', mask1=umask, & & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) !CALL timing_start('dyn_ldf_iso_vslab') #if defined key_z_first ! ! =============== DO jj = 2, jpjm1 ! Vertical slab ! ! =============== ! I. vertical trends associated with the lateral mixing ! ===================================================== ! (excluding the vertical flux proportional to dk[t] ! I.1 horizontal momentum gradient ! -------------------------------- DO ji = 2, jpi DO jk = 1, mbkmax(ji,jj) ! jpk ! i-gradient of u at jj zdiu (ji,jk) = tmask(ji,jj ,jk) * ( ub(ji,jj ,jk) - ub(ji-1,jj ,jk) ) ! j-gradient of u and v at jj zdju (ji,jk) = fmask(ji,jj ,jk) * ( ub(ji,jj+1,jk) - ub(ji ,jj ,jk) ) zdjv (ji,jk) = tmask(ji,jj ,jk) * ( vb(ji,jj ,jk) - vb(ji ,jj-1,jk) ) ! j-gradient of u and v at jj+1 zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj ,jk) - ub(ji ,jj-1,jk) ) zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji ,jj ,jk) ) END DO END DO DO ji = 1, jpim1 DO jk = 1, mbkmax(ji,jj) ! jpk ! i-gradient of v at jj zdiv (ji,jk) = fmask(ji,jj ,jk) * ( vb(ji+1,jj,jk) - vb(ji ,jj ,jk) ) END DO END DO ! I.2 Vertical fluxes ! ------------------- ! Surface and bottom vertical fluxes set to zero !!$ DO ji = 1, jpi !!$ zfuw(ji, 1 ) = 0.e0 !!$ zfvw(ji, 1 ) = 0.e0 !!$ zfuw(ji,jpk) = 0.e0 !!$ zfvw(ji,jpk) = 0.e0 !!$ END DO ! interior (2=