MODULE dynldf_bilapg !!====================================================================== !! *** MODULE dynldf_bilapg *** !! Ocean dynamics: lateral viscosity trend !!====================================================================== !! History : OPA ! 1997-07 (G. Madec) Original code !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module !! 2.0 ! 2004-08 (C. Talandier) New trends organization !!---------------------------------------------------------------------- #if defined key_ldfslp || defined key_esopa !!---------------------------------------------------------------------- !! 'key_ldfslp' Rotation of mixing tensor !!---------------------------------------------------------------------- !! dyn_ldf_bilapg : update the momentum trend with the horizontal part !! of the horizontal s-coord. bilaplacian diffusion !! ldfguv : !!---------------------------------------------------------------------- USE oce ! ocean dynamics and tracers USE dom_oce ! ocean space and time domain USE ldfdyn_oce ! ocean dynamics lateral physics USE zdf_oce ! ocean vertical physics USE trdmod ! ocean dynamics trends USE trdmod_oce ! ocean variables trends USE ldfslp ! iso-neutral slopes available USE in_out_manager ! I/O manager USE lib_mpp ! MPP library USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE prtctl ! Print control IMPLICIT NONE PRIVATE PUBLIC dyn_ldf_bilapg ! called by step.F90 !FTRANS zfuw zfvw zdiu zdiv :I :z REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zfvw , zdiu, zdiv ! 2D workspace (ldfguv) !FTRANS zdju zdj1u zdjv zdj1v :I :z REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zdju, zdj1u, zdjv, zdj1v ! 2D workspace (ldfguv) !! * Control permutation of array indices # include "oce_ftrans.h90" # include "dom_oce_ftrans.h90" # include "ldfdyn_oce_ftrans.h90" # include "zdf_oce_ftrans.h90" # include "ldfslp_ftrans.h90" !! * Substitutions # include "domzgr_substitute.h90" # include "ldfdyn_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OPA 3.3 , NEMO Consortium (2010) !! $Id$ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS INTEGER FUNCTION dyn_ldf_bilapg_alloc() !!---------------------------------------------------------------------- !! *** ROUTINE dyn_ldf_bilapg_alloc *** !!---------------------------------------------------------------------- ALLOCATE( zfuw(jpi,jpk) , zfvw (jpi,jpk) , zdiu(jpi,jpk) , zdiv (jpi,jpk) , & & zdju(jpi,jpk) , zdj1u(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_bilapg_alloc ) ! IF( dyn_ldf_bilapg_alloc /= 0 ) CALL ctl_warn('dyn_ldf_bilapg_alloc: failed to allocate arrays') END FUNCTION dyn_ldf_bilapg_alloc SUBROUTINE dyn_ldf_bilapg( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE dyn_ldf_bilapg *** !! !! ** Purpose : Compute the before trend of the horizontal momentum !! diffusion and add it to the general trend of momentum equation. !! !! ** Method : The lateral momentum diffusive trends is provided by a !! a 4th order operator rotated along geopotential surfaces. It is !! computed using before fields (forward in time) and geopotential !! slopes computed in routine inildf. !! -1- compute the geopotential harmonic operator applied to !! (ub,vb) and multiply it by the eddy diffusivity coefficient !! (done by a call to ldfgpu and ldfgpv routines) The result is in !! (zwk1,zwk2) arrays. Applied the domain lateral boundary conditions !! by call to lbc_lnk. !! -2- applied to (zwk1,zwk2) the geopotential harmonic operator !! by a second call to ldfgpu and ldfgpv routines respectively. The !! result is in (zwk3,zwk4) arrays. !! -3- Add this trend to the general trend (ta,sa): !! (ua,va) = (ua,va) + (zwk3,zwk4) !! 'key_trddyn' defined: the trend is saved for diagnostics. !! !! ** Action : - Update (ua,va) arrays with the before geopotential !! biharmonic mixing trend. !! - save the trend in (zwk3,zwk4) ('key_trddyn') !!---------------------------------------------------------------------- USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released USE wrk_nemo, ONLY: zwk1 => wrk_3d_3 , zwk2 => wrk_3d_4 ! 3D workspace USE oce , ONLY: zwk3 => ta , zwk4 => sa ! ta, sa used as 3D workspace !! DCSE_NEMO: need additional directives for renamed module variables !FTRANS zwk1 :I :I :z !FTRANS zwk2 :I :I :z !FTRANS zwk3 :I :I :z !FTRANS zwk4 :I :I :z ! INTEGER, INTENT( in ) :: kt ! ocean time-step index ! INTEGER :: ji, jj, jk ! dummy loop indices !!---------------------------------------------------------------------- IF( wrk_in_use(3, 3,4) ) THEN CALL ctl_stop('dyn_ldf_bilapg: requested workspace arrays unavailable') ; RETURN ENDIF IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn_ldf_bilapg : horizontal biharmonic operator in s-coordinate' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~' zwk1(:,:,:) = 0.e0_wp ; zwk3(:,:,:) = 0.e0_wp zwk2(:,:,:) = 0.e0_wp ; zwk4(:,:,:) = 0.e0_wp ! ! allocate dyn_ldf_bilapg arrays IF( dyn_ldf_bilapg_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_ldf_bilapg: failed to allocate arrays') ENDIF ! Laplacian of (ub,vb) multiplied by ahm ! -------------------------------------- CALL ldfguv( ub, vb, zwk1, zwk2, 1 ) ! rotated harmonic operator applied to (ub,vb) ! ! and multiply by ahmu, ahmv (output in (zwk1,zwk2) ) CALL lbc_lnk( zwk1, 'U', -1. ) ; CALL lbc_lnk( zwk2, 'V', -1. ) ! Lateral boundary conditions ! Bilaplacian of (ub,vb) ! ---------------------- CALL ldfguv( zwk1, zwk2, zwk3, zwk4, 2 ) ! rotated harmonic operator applied to (zwk1,zwk2) ! ! (output in (zwk3,zwk4) ) ! Update the momentum trends ! -------------------------- #if defined key_z_first DO jj = 2, jpjm1 ! add the diffusive trend to the general momentum trends DO ji = 2, jpim1 DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1 #else DO jj = 2, jpjm1 ! add the diffusive trend to the general momentum trends DO jk = 1, jpkm1 DO ji = 2, jpim1 #endif ua(ji,jj,jk) = ua(ji,jj,jk) + zwk3(ji,jj,jk) va(ji,jj,jk) = va(ji,jj,jk) + zwk4(ji,jj,jk) END DO END DO END DO ! IF( wrk_not_released(3, 3,4) ) CALL ctl_stop('dyn_ldf_bilapg: failed to release workspace arrays') ! END SUBROUTINE dyn_ldf_bilapg SUBROUTINE ldfguv( pu, pv, plu, plv, kahm ) !!---------------------------------------------------------------------- !! *** ROUTINE ldfguv *** !! !! ** Purpose : Apply a geopotential harmonic operator to (pu,pv) !! (defined at u- and v-points) and multiply it by the eddy !! viscosity coefficient (if kahm=1). !! !! ** Method : The harmonic operator rotated along geopotential !! surfaces is applied to (pu,pv) using the slopes of geopotential !! surfaces computed in inildf routine. The result is provided in !! (plu,plv) arrays. It is computed in 2 stepv: !! !! First step: horizontal part of the operator. It is computed on !! ========== pu as follows (idem on pv) !! horizontal fluxes : !! zftu = e2u*e3u/e1u di[ pu ] - e2u*uslp dk[ mi(mk(pu)) ] !! zftv = e1v*e3v/e2v dj[ pu ] - e1v*vslp dk[ mj(mk(pu)) ] !! take the horizontal divergence of the fluxes (no divided by !! the volume element : !! plu = di-1[ zftu ] + dj-1[ zftv ] !! !! Second step: vertical part of the operator. It is computed on !! =========== pu as follows (idem on pv) !! vertical fluxes : !! zftw = e1t*e2t/e3w * (wslpi^2+wslpj^2) dk-1[ pu ] !! - e2t * wslpi di[ mi(mk(pu)) ] !! - e1t * wslpj dj[ mj(mk(pu)) ] !! take the vertical divergence of the fluxes add it to the hori- !! zontal component, divide the result by the volume element and !! if kahm=1, multiply by the eddy diffusivity coefficient: !! plu = aht / (e1t*e2t*e3t) { plu + dk[ zftw ] } !! else: !! plu = 1 / (e1t*e2t*e3t) { plu + dk[ zftw ] } !! !! ** Action : !! plu, plv : partial harmonic operator applied to !! pu and pv (all the components except !! second order vertical derivative term) !! 'key_trddyn' defined: the trend is saved for diagnostics. !!---------------------------------------------------------------------- USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released USE wrk_nemo, ONLY: ziut => wrk_2d_1 , zjuf => wrk_2d_2 , zjvt => wrk_2d_3 USE wrk_nemo, ONLY: zivf => wrk_2d_4 #if ! defined key_z_first USE wrk_nemo, ONLY: zdku => wrk_2d_5 , zdk1u => wrk_2d_6 USE wrk_nemo, ONLY: zdkv => wrk_2d_7 , zdk1v => wrk_2d_8 #endif USE timing , ONLY: timing_start, timing_stop !! !FTRANS pu :I :I :z !FTRANS pv :I :I :z !FTRANS plu :I :I :z !FTRANS plv :I :I :z !! DCSE_NEMO: work around deficiency in ftrans ! REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pu , pv ! 1st call: before horizontal velocity ! ! 2nd call: ahm x these fields REAL(wp), INTENT(in ) :: pu(jpi,jpj,jpkorig) , pv(jpi,jpj,jpkorig) ! REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: plu, plv ! partial harmonic operator applied to ! ! pu and pv (all the components except ! ! second order vertical derivative term) REAL(wp), INTENT( out) :: plu(jpi,jpj,jpkorig), plv(jpi,jpj,jpkorig) ! partial harmonic operator applied to INTEGER , INTENT(in ) :: kahm ! =1 1st call ; =2 2nd call ! INTEGER :: ji, jj, jk ! dummy loop indices INTEGER :: jif, jjf ! dummy loop indices over full domain REAL(wp) :: zabe1 , zabe2 , zcof1 , zcof2 ! local scalar REAL(wp) :: zcoef0, zcoef3, zcoef4 ! - - REAL(wp) :: zbur, zbvr, zmkt, zmkf, zuav, zvav ! - - REAL(wp) :: zuwslpi, zuwslpj, zvwslpi, zvwslpj ! - - #if defined key_z_first ! Can use scalars instead of work arrays when built with z-first REAL(wp) :: zdku, zdkv, zdk1u, zdk1v #endif !!---------------------------------------------------------------------- CALL timing_start('ldfguv') IF( wrk_in_use(2, 1,2,3,4,5,6,7,8) ) THEN CALL ctl_stop('dyn:ldfguv: requested workspace arrays unavailable') ; RETURN END IF CALL timing_start('ldfguv_1st') #if defined key_z_first ! ! ********** ! ! First step ! ! ! ********** ! DO jj = 2, jpjm1 DO ji = 2, jpim1 ! Treat jk = 1 separately as is special case jk = 1 plu(ji,jj,jk) = & ! ------------- ziut (ji+1, jj) - tmask(ji+1,jj,jk) * & ( (e2t(ji+1,jj) * fse3t(ji+1,jj,jk) / e1t(ji+1,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) ) & ! + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & + (-e2t(ji+1,jj) / 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) ) ) * & ( ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + & ( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) + & ! +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + & ( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) & ) ) - & ! ------------- ziut (ji,jj ) + tmask(ji,jj,jk) * & ( (e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) ) & ! + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & + (-e2t(ji,jj) / MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1._wp ) & * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) ) * & ( ( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji,jj,jk+1) + & ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) + & ! +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) ( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji,jj,jk+1) + & ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) & ) ) + & ! ------------- zjuf (ji ,jj) - fmask(ji,jj,jk) * & ( (e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) ) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) ) & ! + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) & + (-e1f(ji,jj) /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) ) ) * & ( & ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) & + ( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj ,jk+1) & ! +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) + ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) & + ( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj,jk+1) & ) ) - & ! ------------- zjuf (ji,jj-1) fmask(ji,jj-1,jk) * & ( (e1f(ji,jj-1) * fse3f(ji,jj-1,jk) / e2f(ji,jj-1) ) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) ) & ! + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) & + (-e1f(ji,jj-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) ) ) * & ( & ( pu(ji,jj,jk) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) & + ( pu(ji,jj-1 ,jk) - pu(ji,jj-1 ,jk+1) ) * umask(ji,jj-1 ,jk+1) & ! +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) + ( pu(ji,jj,jk) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) & + ( pu(ji,jj-1 ,jk) - pu(ji,jj-1 ,jk+1) ) * umask(ji,jj-1,jk+1) & ) ) plv(ji,jj,jk) = & ! ------------- zivf (ji,jj ) - fmask(ji,jj,jk) * & ( (e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) ) & ! + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj) & + ((-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. )) & * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * ( & ( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) + & ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + & ! +zdk1u(ji,jj) + zdku (ji+1,jj) ) ) ( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) + & ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) & ) ) - & ! ------------- zivf (ji-1,jj) + fmask(ji-1,jj,jk) * & ( (e2f(ji-1,jj) * fse3f(ji-1,jj,jk) / e1f(ji-1,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) ) & ! + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj) & + ((-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. )) & * 0.5 * ( uslp(ji-1,jj+1,jk) + uslp(ji-1,jj,jk) )) * ( & ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) + & ( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) + & ! +zdk1u(ji,jj) + zdku (ji+1,jj) ) ) ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) + & ( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) & ) ) + & ! ------------- zjvt (ji,jj+1) - tmask(ji,jj+1,jk) * & ( (e1t(ji,jj+1) * fse3t(ji,jj+1,jk) / e2t(ji,jj+1)) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) ) & ! + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj) & + ((-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. ) ) & * 0.5 * ( vslp(ji,jj,jk) + vslp(ji,jj+1,jk) )) * ( & ( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj ,jk+1) + & ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) + & ! +zdk1u(ji,jj-1) + zdku (ji,jj) ) ) ( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj ,jk+1) + & ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) & ) ) - & ! ------------- zjvt (ji,jj ) tmask(ji,jj,jk) * & ( (e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) ) & ! + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj) & + ((-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. ) ) & * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * ( & ( pu(ji,jj-1,jk) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + & ( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj,jk+1) + & ! +zdk1u(ji,jj-1) + zdku (ji,jj) ) ) ( pu(ji,jj-1,jk) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + & ( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj,jk+1) & ) ) DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1 ! plu(ji,jj,jk) = ziut (ji+1,jj) - & ! ziut (ji,jj ) + & ! zjuf (ji ,jj) - & ! zjuf (ji,jj-1) plu(ji,jj,jk) = & ! ------------- ziut (ji+1, jj ) - tmask(ji+1,jj,jk) * & ( (e2t(ji+1,jj) * fse3t(ji+1,jj,jk) / e1t(ji+1,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) ) & ! + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & + (-e2t(ji+1,jj) / 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) ) ) * & ( ( pu(ji+1,jj,jk-1) - pu(ji+1,jj,jk ) ) * umask(ji+1,jj,jk) + & ( pu(ji ,jj,jk ) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) + & ! +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) ( pu(ji+1,jj,jk ) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + & ( pu(ji ,jj,jk-1) - pu(ji ,jj,jk ) ) * umask(ji ,jj,jk) & ) ) - & ! ------------- ziut (ji , jj ) + tmask(ji,jj,jk) * & ( (e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) ) & ! + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & + (-e2t(ji,jj) / MAX(umask(ji-1,jj,jk)+umask(ji,jj,jk+1) & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk), 1._wp ) & * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) ) * & ( ( pu(ji ,jj,jk-1) - pu(ji ,jj,jk ) ) * umask(ji,jj,jk) + & ( pu(ji-1,jj,jk ) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) + & ! +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) ( pu(ji ,jj,jk ) - pu(ji ,jj,jk+1) ) * umask(ji,jj,jk+1) + & ( pu(ji-1,jj,jk-1) - pu(ji-1,jj,jk ) ) * umask(ji-1,jj,jk) & ) ) + & ! ------------- zjuf (ji , jj ) - fmask(ji,jj,jk) * & ( (e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) ) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) ) & ! + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) & + (-e1f(ji,jj) /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) ) ) * & ( & (pu(ji,jj+1,jk-1) - pu(ji,jj+1,jk ) ) * umask(ji,jj+1,jk) + & (pu(ji,jj ,jk ) - pu(ji,jj ,jk+1) ) * umask(ji,jj,jk+1) + & ! +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) (pu(ji,jj+1,jk ) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) + & (pu(ji,jj ,jk-1) - pu(ji,jj ,jk ) ) * umask(ji,jj,jk) & ) ) - & ! ------------- zjuf (ji , jj-1) fmask(ji,jj-1,jk) * & ( (e1f(ji,jj-1) * fse3f(ji,jj-1,jk) / e2f(ji,jj-1) ) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) ) & ! + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) & + (-e1f(ji,jj-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) ) ) * & ( & (pu(ji,jj,jk-1) - pu(ji,jj,jk ) ) * umask(ji,jj,jk) + & (pu(ji,jj-1,jk ) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + & ! +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) (pu(ji,jj,jk ) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) + & (pu(ji,jj-1,jk-1) - pu(ji,jj-1 ,jk ) ) * umask(ji,jj-1,jk) & ) ) ! plv(ji,jj,jk) = zivf (ji,jj ) - & ! zivf (ji-1,jj) + & ! zjvt (ji,jj+1) - & ! zjvt (ji,jj ) plv(ji,jj,jk) = & ! ------------- zivf (ji,jj ) - fmask(ji,jj,jk) * & ( (e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) ) & ! + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj) & + ((-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. )) & * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * ( & ( pu(ji ,jj,jk-1) - pu(ji ,jj,jk ) ) * umask(ji ,jj,jk ) + & ( pu(ji+1,jj,jk ) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + & ! +zdk1u(ji,jj) + zdku (ji+1,jj) ) ) ( pu(ji ,jj,jk ) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) + & ( pu(ji+1,jj,jk-1) - pu(ji+1,jj,jk ) ) * umask(ji+1,jj,jk ) & ) ) - & ! ------------- zivf (ji-1,jj) + fmask(ji-1,jj,jk) * & ( (e2f(ji-1,jj) * fse3f(ji-1,jj,jk) / e1f(ji-1,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) ) & ! + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj) & + ((-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. )) & * 0.5 * ( uslp(ji-1,jj+1,jk) + uslp(ji-1,jj,jk) )) * ( & ( pu(ji-1 ,jj,jk-1) - pu(ji-1 ,jj,jk ) ) * umask(ji-1 ,jj,jk ) + & ( pu(ji,jj,jk ) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) + & ! +zdk1u(ji,jj) + zdku (ji+1,jj) ) ) ( pu(ji-1 ,jj,jk ) - pu(ji-1 ,jj,jk+1) ) * umask(ji-1 ,jj,jk+1) + & ( pu(ji,jj,jk-1) - pu(ji,jj,jk ) ) * umask(ji,jj,jk ) & ) ) + & ! ------------- zjvt (ji,jj+1) - tmask(ji,jj+1,jk) * & ( (e1t(ji,jj+1) * fse3t(ji,jj+1,jk) / e2t(ji,jj+1)) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) ) & ! + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj) & + ((-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. ) ) & * 0.5 * ( vslp(ji,jj,jk) + vslp(ji,jj+1,jk) )) * & ( & ( pu(ji,jj,jk-1) - pu(ji,jj,jk ) ) * umask(ji,jj,jk) + & ( pu(ji,jj+1 ,jk ) - pu(ji,jj+1 ,jk+1) ) * umask(ji,jj+1,jk+1) + & ! +zdk1u(ji,jj-1) + zdku (ji,jj) ) ) ( pu(ji,jj,jk ) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) + & ( pu(ji,jj+1 ,jk-1) - pu(ji,jj+1 ,jk ) ) * umask(ji,jj+1,jk) & ) ) - & ! ------------- zjvt (ji,jj ) tmask(ji,jj,jk) * & ( (e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) ) & ! + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj) & + ((-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. ) ) & * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * & ( & ( pu(ji,jj-1,jk-1) - pu(ji,jj-1,jk ) ) * umask(ji,jj-1,jk) + & ( pu(ji,jj ,jk ) - pu(ji,jj ,jk+1) ) * umask(ji,jj,jk+1) + & ! +zdk1u(ji,jj-1) + zdku (ji,jj) ) ) ( pu(ji,jj-1,jk ) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + & ( pu(ji,jj ,jk-1) - pu(ji,jj ,jk ) ) * umask(ji,jj,jk) ) ) END DO END DO ! END DO #else ! ! ********** ! ! =============== DO jk = 1, jpkm1 ! First step ! ! Horizontal slab ! ! ********** ! ! =============== ! I.1 Vertical gradient of pu and pv at level jk and jk+1 ! ------------------------------------------------------- ! surface boundary condition: zdku(jk=1)=zdku(jk=2) ! zdkv(jk=1)=zdkv(jk=2) zdk1u(:,:) = ( pu(:,:,jk) - pu(:,:,jk+1) ) * umask(:,:,jk+1) zdk1v(:,:) = ( pv(:,:,jk) - pv(:,:,jk+1) ) * vmask(:,:,jk+1) IF( jk == 1 ) THEN zdku(:,:) = zdk1u(:,:) zdkv(:,:) = zdk1v(:,:) ELSE zdku(:,:) = ( pu(:,:,jk-1) - pu(:,:,jk) ) * umask(:,:,jk) zdkv(:,:) = ( pv(:,:,jk-1) - pv(:,:,jk) ) * vmask(:,:,jk) ENDIF ! -----f----- ! I.2 Horizontal fluxes on U | ! ------------------------=== t u t ! | ! i-flux at t-point -----f----- DO jj = 1, jpjm1 DO ji = 2, jpi zabe1 = e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj) zmkt = 1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. ) zcof1 = -e2t(ji,jj) * zmkt & * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) ziut(ji,jj) = tmask(ji,jj,jk) * & ( zabe1 * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) ) & + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) END DO END DO ! j-flux at f-point DO jj = 1, jpjm1 DO ji = 1, jpim1 zabe2 = e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) zmkf = 1./MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) & + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1. ) zcof2 = -e1f(ji,jj) * zmkf & * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) zjuf(ji,jj) = fmask(ji,jj,jk) * & ( zabe2 * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) ) & + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) & +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) END DO END DO ! | t | ! I.3 Horizontal fluxes on V | | ! ------------------------=== f---v---f ! | | ! i-flux at f-point | t | DO jj = 1, jpjm1 DO ji = 1, jpim1 zabe1 = e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj) zmkf = 1./MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1. ) zcof1 = -e2f(ji,jj) * zmkf & * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) zivf(ji,jj) = fmask(ji,jj,jk) * & ( zabe1 * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) ) & + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj) & +zdk1u(ji,jj) + zdku (ji+1,jj) ) ) END DO END DO ! j-flux at t-point DO jj = 2, jpj DO ji = 1, jpim1 zabe2 = e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj) zmkt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. ) zcof2 = -e1t(ji,jj) * zmkt & * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) zjvt(ji,jj) = tmask(ji,jj,jk) * & ( zabe2 * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) ) & + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj) & +zdk1u(ji,jj-1) + zdku (ji,jj) ) ) END DO END DO ! I.4 Second derivative (divergence) (not divided by the volume) ! --------------------- DO jj = 2, jpjm1 DO ji = 2, jpim1 plu(ji,jj,jk) = ziut (ji+1,jj) - ziut (ji,jj ) & + zjuf (ji ,jj) - zjuf (ji,jj-1) plv(ji,jj,jk) = zivf (ji,jj ) - zivf (ji-1,jj) & + zjvt (ji,jj+1) - zjvt (ji,jj ) END DO END DO ! ! =============== END DO ! End of slab ! ! =============== #endif CALL timing_stop('ldfguv_1st','section') !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, CALL timing_start('ldfguv_2nd') ! ! ************ ! ! =============== DO jj = 2, jpjm1 ! Second step ! ! Horizontal slab ! ! ************ ! ! =============== ! II.1 horizontal (pu,pv) gradients ! --------------------------------- #if defined key_z_first DO ji = 2, jpi DO jk = 1, mbkmax(ji,jj) ! jpk ! i-gradient of u at jj zdiu (ji,jk) = tmask(ji,jj ,jk) * ( pu(ji,jj ,jk) - pu(ji-1,jj ,jk) ) ! j-gradient of u and v at jj zdju (ji,jk) = fmask(ji,jj ,jk) * ( pu(ji,jj+1,jk) - pu(ji ,jj ,jk) ) zdjv (ji,jk) = tmask(ji,jj ,jk) * ( pv(ji,jj ,jk) - pv(ji ,jj-1,jk) ) ! j-gradient of u and v at jj+1 zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( pu(ji,jj ,jk) - pu(ji ,jj-1,jk) ) zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pv(ji,jj+1,jk) - pv(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) * ( pv(ji+1,jj,jk) - pv(ji ,jj ,jk) ) END DO END DO #else DO jk = 1, jpk DO ji = 2, jpi ! i-gradient of u at jj zdiu (ji,jk) = tmask(ji,jj ,jk) * ( pu(ji,jj ,jk) - pu(ji-1,jj ,jk) ) ! j-gradient of u and v at jj zdju (ji,jk) = fmask(ji,jj ,jk) * ( pu(ji,jj+1,jk) - pu(ji ,jj ,jk) ) zdjv (ji,jk) = tmask(ji,jj ,jk) * ( pv(ji,jj ,jk) - pv(ji ,jj-1,jk) ) ! j-gradient of u and v at jj+1 zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( pu(ji,jj ,jk) - pu(ji ,jj-1,jk) ) zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pv(ji,jj+1,jk) - pv(ji ,jj ,jk) ) END DO END DO DO jk = 1, jpk DO ji = 1, jpim1 ! i-gradient of v at jj zdiv (ji,jk) = fmask(ji,jj ,jk) * ( pv(ji+1,jj,jk) - pv(ji ,jj ,jk) ) END DO END DO #endif ! II.2 Vertical fluxes ! -------------------- ! Surface and bottom vertical fluxes set to zero #if defined key_z_first DO ji=1, jpi zfuw(ji, 1 ) = 0.e0 zfvw(ji, 1 ) = 0.e0 zfuw(ji,mbkmax(ji,jj):jpk) = 0.e0 zfvw(ji,mbkmax(ji,jj):jpk) = 0.e0 END DO #else zfuw(:, 1 ) = 0.e0 zfvw(:, 1 ) = 0.e0 zfuw(:,jpk) = 0.e0 zfvw(:,jpk) = 0.e0 #endif ! interior (2=