- Timestamp:
- 2015-11-30T20:55:41+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r5944 r5956 11 11 !! 3.3 ! 2010-10 (G. Nurser, C. Harris, G. Madec) add Griffies operator 12 12 !! - ! 2010-11 (F. Dupond, G. Madec) bug correction in slopes just below the ML 13 !! 3.7 ! 2013-12 (F. Lemarie, G. Madec) add limiter on triad slopes 13 14 !!---------------------------------------------------------------------- 14 #if defined key_ldfslp || defined key_esopa 15 15 16 !!---------------------------------------------------------------------- 16 !! 'key_ldfslp' Rotation of lateral mixing tensor17 !!----------------------------------------------------------------------18 !! ldf_slp_grif : calculates the triads of isoneutral slopes (Griffies operator)19 17 !! ldf_slp : calculates the slopes of neutral surface (Madec operator) 18 !! ldf_slp_triad : calculates the triads of isoneutral slopes (Griffies operator) 20 19 !! ldf_slp_mxl : calculates the slopes at the base of the mixed layer (Madec operator) 21 20 !! ldf_slp_init : initialization of the slopes computation … … 23 22 USE oce ! ocean dynamics and tracers 24 23 USE dom_oce ! ocean space and time domain 25 USE ldftra_oce ! lateral diffusion: traceur 26 USE ldfdyn_oce ! lateral diffusion: dynamics 24 USE ldfdyn ! lateral diffusion: eddy viscosity coef. 27 25 USE phycst ! physical constants 28 26 USE zdfmxl ! mixed layer depth … … 30 28 ! 31 29 USE in_out_manager ! I/O manager 30 USE prtctl ! Print control 32 31 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 33 USE prtctl ! Print control 32 USE lib_mpp ! distribued memory computing library 33 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 34 34 USE wrk_nemo ! work arrays 35 35 USE timing ! Timing 36 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)37 36 38 37 IMPLICIT NONE 39 38 PRIVATE 40 39 41 PUBLIC ldf_slp ! routine called by step.F90 42 PUBLIC ldf_slp_grif ! routine called by step.F90 43 PUBLIC ldf_slp_init ! routine called by opa.F90 44 45 LOGICAL , PUBLIC, PARAMETER :: lk_ldfslp = .TRUE. !: slopes flag 46 ! !! Madec operator 47 ! Arrays allocated in ldf_slp_init() routine once we know whether we're using the Griffies or Madec operator 40 PUBLIC ldf_slp ! routine called by step.F90 41 PUBLIC ldf_slp_triad ! routine called by step.F90 42 PUBLIC ldf_slp_init ! routine called by nemogcm.F90 43 44 LOGICAL , PUBLIC :: l_ldfslp = .FALSE. !: slopes flag 45 46 LOGICAL , PUBLIC :: ln_traldf_iso = .TRUE. !: iso-neutral direction 47 LOGICAL , PUBLIC :: ln_traldf_triad = .FALSE. !: griffies triad scheme 48 49 LOGICAL , PUBLIC :: ln_triad_iso = .FALSE. !: pure horizontal mixing in ML 50 LOGICAL , PUBLIC :: ln_botmix_triad = .FALSE. !: mixing on bottom 51 REAL(wp), PUBLIC :: rn_sw_triad = 1._wp !: =1 switching triads ; =0 all four triads used 52 REAL(wp), PUBLIC :: rn_slpmax = 0.01_wp !: slope limit 53 54 LOGICAL , PUBLIC :: l_grad_zps = .FALSE. !: special treatment for Horz Tgradients w partial steps (triad operator) 55 56 ! !! Classic operator (Madec) 48 57 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: uslp, wslpi !: i_slope at U- and W-points 49 58 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: vslp, wslpj !: j-slope at V- and W-points 50 ! !! Griffies operator59 ! !! triad operator (Griffies) 51 60 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wslp2 !: wslp**2 from Griffies quarter cells 52 61 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) :: triadi_g, triadj_g !: skew flux slopes relative to geopotentials 53 62 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) :: triadi , triadj !: isoneutral slopes relative to model-coordinate 54 55 ! !! Madec operator 56 ! Arrays allocated in ldf_slp_init() routine once we know whether we're using the Griffies or Madec operator 63 ! !! both operators 64 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ah_wslp2 !: ah * slope^2 at w-point 65 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akz !: stabilizing vertical diffusivity 66 67 ! !! Madec operator 57 68 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: omlmask ! mask of the surface mixed layer at T-pt 58 69 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: uslpml, wslpiml ! i_slope at U- and W-points just below the mixed layer … … 63 74 !! * Substitutions 64 75 # include "domzgr_substitute.h90" 65 # include "ldftra_substitute.h90"66 # include "ldfeiv_substitute.h90"67 76 # include "vectopt_loop_substitute.h90" 68 77 !!---------------------------------------------------------------------- 69 !! NEMO/OPA 4.0 , NEMO Consortium (201 1)78 !! NEMO/OPA 4.0 , NEMO Consortium (2014) 70 79 !! $Id$ 71 80 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 105 114 INTEGER :: ii0, ii1 ! temporary integer 106 115 INTEGER :: ij0, ij1 ! temporary integer 107 REAL(wp) :: zeps, zm1_g, zm1_2g, z1_16, zcofw ! local scalars116 REAL(wp) :: zeps, zm1_g, zm1_2g, z1_16, zcofw, z1_slpmax ! local scalars 108 117 REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - - 109 118 REAL(wp) :: zcj, zfj, zav, zbv, zaj, zbj ! - - 110 119 REAL(wp) :: zck, zfk, zbw ! - - 111 REAL(wp) :: zdepv, zdepu ! - -120 REAL(wp), POINTER, DIMENSION(:,: ) :: zslpml_hmlpu, zslpml_hmlpv 112 121 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz, zww 113 122 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdzr 114 123 REAL(wp), POINTER, DIMENSION(:,:,:) :: zgru, zgrv 115 REAL(wp), POINTER, DIMENSION(:,: ) :: zhmlpu, zhmlpv116 124 !!---------------------------------------------------------------------- 117 125 ! … … 119 127 ! 120 128 CALL wrk_alloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv ) 121 CALL wrk_alloc( jpi,jpj, zhmlpu, zhmlpv ) 122 123 IF ( ln_traldf_iso .OR. ln_dynldf_iso ) THEN 124 125 zeps = 1.e-20_wp !== Local constant initialization ==! 126 z1_16 = 1.0_wp / 16._wp 127 zm1_g = -1.0_wp / grav 128 zm1_2g = -0.5_wp / grav 129 ! 130 zww(:,:,:) = 0._wp 131 zwz(:,:,:) = 0._wp 132 ! 133 DO jk = 1, jpk !== i- & j-gradient of density ==! 134 DO jj = 1, jpjm1 135 DO ji = 1, fs_jpim1 ! vector opt. 136 zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj ,jk) - prd(ji,jj,jk) ) 137 zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji ,jj+1,jk) - prd(ji,jj,jk) ) 138 END DO 139 END DO 140 END DO 141 IF( ln_zps ) THEN ! partial steps correction at the bottom ocean level 142 DO jj = 1, jpjm1 143 DO ji = 1, jpim1 144 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 145 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 146 END DO 147 END DO 148 ENDIF 149 IF( ln_zps .AND. ln_isfcav ) THEN ! partial steps correction at the bottom ocean level 150 DO jj = 1, jpjm1 151 DO ji = 1, jpim1 129 CALL wrk_alloc( jpi,jpj, zslpml_hmlpu, zslpml_hmlpv ) 130 131 zeps = 1.e-20_wp !== Local constant initialization ==! 132 z1_16 = 1.0_wp / 16._wp 133 zm1_g = -1.0_wp / grav 134 zm1_2g = -0.5_wp / grav 135 z1_slpmax = 1._wp / rn_slpmax 136 ! 137 zww(:,:,:) = 0._wp 138 zwz(:,:,:) = 0._wp 139 ! 140 DO jk = 1, jpk !== i- & j-gradient of density ==! 141 DO jj = 1, jpjm1 142 DO ji = 1, fs_jpim1 ! vector opt. 143 zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj ,jk) - prd(ji,jj,jk) ) 144 zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji ,jj+1,jk) - prd(ji,jj,jk) ) 145 END DO 146 END DO 147 END DO 148 IF( ln_zps ) THEN ! partial steps correction at the bottom ocean level 149 DO jj = 1, jpjm1 150 DO ji = 1, jpim1 151 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 152 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 153 END DO 154 END DO 155 ENDIF 156 IF( ln_zps .AND. ln_isfcav ) THEN ! partial steps correction at the bottom ocean level 157 DO jj = 1, jpjm1 158 DO ji = 1, jpim1 152 159 IF ( miku(ji,jj) > 1 ) zgru(ji,jj,miku(ji,jj)) = grui(ji,jj) 153 160 IF ( mikv(ji,jj) > 1 ) zgrv(ji,jj,mikv(ji,jj)) = grvi(ji,jj) 154 END DO 155 END DO 156 ENDIF 157 ! 158 !== Local vertical density gradient at T-point == ! (evaluated from N^2) 159 ! interior value 160 DO jk = 2, jpkm1 161 ! ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point 162 ! ! trick: tmask(ik ) = 0 => all pn2 = 0 => zdzr = 0 163 ! ! else tmask(ik+1) = 0 => pn2(ik+1) = 0 => zdzr divides by 1 164 ! ! tmask(ik+1) /= 0 => all pn2 /= 0 => zdzr divides by 2 165 ! ! NB: 1/(tmask+1) = (1-.5*tmask) substitute a / by a * ==> faster 166 zdzr(:,:,jk) = zm1_g * ( prd(:,:,jk) + 1._wp ) & 167 & * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask(:,:,jk+1) ) 168 END DO 169 ! surface initialisation 170 zdzr(:,:,1) = 0._wp 171 IF ( ln_isfcav ) THEN 172 ! if isf need to overwrite the interior value at at the first ocean point 173 DO jj = 1, jpjm1 174 DO ji = 1, jpim1 175 zdzr(ji,jj,1:mikt(ji,jj)) = 0._wp 176 END DO 177 END DO 178 END IF 179 ! 180 ! !== Slopes just below the mixed layer ==! 181 CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr ) ! output: uslpml, vslpml, wslpiml, wslpjml 182 183 184 ! I. slopes at u and v point | uslp = d/di( prd ) / d/dz( prd ) 185 ! =========================== | vslp = d/dj( prd ) / d/dz( prd ) 186 ! 187 IF ( ln_isfcav ) THEN 188 DO jj = 2, jpjm1 189 DO ji = fs_2, fs_jpim1 ! vector opt. 190 zhmlpu(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji+1,jj ), 5._wp) - 0.5 * ( risfdep(ji,jj) + risfdep(ji+1,jj ) ) 191 zhmlpv(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji ,jj+1), 5._wp) - 0.5 * ( risfdep(ji,jj) + risfdep(ji ,jj+1) ) 192 END DO 193 END DO 194 ELSE 195 DO jj = 2, jpjm1 196 DO ji = fs_2, fs_jpim1 ! vector opt. 197 zhmlpu(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji+1,jj ), 5._wp) 198 zhmlpv(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji ,jj+1), 5._wp) 199 END DO 200 END DO 201 END IF 202 ! PM + GM can be optimized by using : zuslp_hmlpu(ji,jj)= uslpml(ji,jj) / zhmlpu (ji,jj) 203 204 DO jk = 2, jpkm1 !* Slopes at u and v points 205 DO jj = 2, jpjm1 206 DO ji = fs_2, fs_jpim1 ! vector opt. 207 ! ! horizontal and vertical density gradient at u- and v-points 208 zau = zgru(ji,jj,jk) / e1u(ji,jj) 209 zav = zgrv(ji,jj,jk) / e2v(ji,jj) 210 zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj ,jk) ) 211 zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji ,jj+1,jk) ) 212 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 213 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 214 zbu = MIN( zbu, -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,jk)* ABS( zau ) ) 215 zbv = MIN( zbv, -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav ) ) 216 ! ! uslp and vslp output in zwz and zww, resp. 217 zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj ,jk) ) 218 zfj = MAX( omlmask(ji,jj,jk), omlmask(ji ,jj+1,jk) ) 219 ! thickness of water column between surface and level k at u/v point 220 zdepu = 0.5_wp * ( ( fsdept (ji,jj,jk) + fsdept (ji+1,jj ,jk) ) & 221 - ( risfdep(ji,jj) + risfdep(ji+1,jj) ) - fse3u(ji,jj,miku(ji,jj)) ) 222 zdepv = 0.5_wp * ( ( fsdept (ji,jj,jk) + fsdept (ji,jj+1,jk) ) & 223 - ( risfdep(ji,jj) + risfdep(ji,jj+1) ) - fse3v(ji,jj,mikv(ji,jj)) ) 224 ! 225 zwz(ji,jj,jk) = ( 1. - zfi) * zau / ( zbu - zeps ) & 226 & + zfi * uslpml(ji,jj) * zdepu / zhmlpu(ji,jj) 227 zwz(ji,jj,jk) = zwz(ji,jj,jk) * wumask(ji,jj,jk) 228 zww(ji,jj,jk) = ( 1. - zfj) * zav / ( zbv - zeps ) & 229 & + zfj * vslpml(ji,jj) * zdepv / zhmlpv(ji,jj) 230 zww(ji,jj,jk) = zww(ji,jj,jk) * wvmask(ji,jj,jk) 231 232 161 END DO 162 END DO 163 ENDIF 164 ! 165 zdzr(:,:,1) = 0._wp !== Local vertical density gradient at T-point == ! (evaluated from N^2) 166 DO jk = 2, jpkm1 167 ! ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point 168 ! ! trick: tmask(ik ) = 0 => all pn2 = 0 => zdzr = 0 169 ! ! else tmask(ik+1) = 0 => pn2(ik+1) = 0 => zdzr divides by 1 170 ! ! umask(ik+1) /= 0 => all pn2 /= 0 => zdzr divides by 2 171 ! ! NB: 1/(tmask+1) = (1-.5*tmask) substitute a / by a * ==> faster 172 zdzr(:,:,jk) = zm1_g * ( prd(:,:,jk) + 1._wp ) & 173 & * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask(:,:,jk+1) ) 174 END DO 175 ! 176 ! !== Slopes just below the mixed layer ==! 177 CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr ) ! output: uslpml, vslpml, wslpiml, wslpjml 178 179 ! I. slopes at u and v point | uslp = d/di( prd ) / d/dz( prd ) 180 ! =========================== | vslp = d/dj( prd ) / d/dz( prd ) 181 ! 182 IF ( ln_isfcav ) THEN 183 DO jj = 2, jpjm1 184 DO ji = fs_2, fs_jpim1 ! vector opt. 185 zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / ( MAX(hmlpt(ji,jj), hmlpt(ji+1,jj ), 5._wp) & 186 & - 0.5 * ( risfdep(ji,jj) + risfdep(ji+1,jj ) ) 187 zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / ( MAX(hmlpt(ji,jj), hmlpt(ji ,jj+1), 5._wp) & 188 & - 0.5 * ( risfdep(ji,jj) + risfdep(ji ,jj+1) ) 189 END DO 190 END DO 191 ELSE 192 DO jj = 2, jpjm1 193 DO ji = fs_2, fs_jpim1 ! vector opt. 194 zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji+1,jj ), 5._wp) 195 zslpml_hmlpu(ji,jj) = vslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji ,jj+1), 5._wp) 196 END DO 197 END DO 198 END IF 199 200 DO jk = 2, jpkm1 !* Slopes at u and v points 201 DO jj = 2, jpjm1 202 DO ji = fs_2, fs_jpim1 ! vector opt. 203 ! ! horizontal and vertical density gradient at u- and v-points 204 zau = zgru(ji,jj,jk) * r1_e1u(ji,jj) 205 zav = zgrv(ji,jj,jk) * r1_e2v(ji,jj) 206 zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj ,jk) ) 207 zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji ,jj+1,jk) ) 208 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 209 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 210 zbu = MIN( zbu, - z1_slpmax * ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,jk)* ABS( zau ) ) 211 zbv = MIN( zbv, - z1_slpmax * ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav ) ) 212 ! ! uslp and vslp output in zwz and zww, resp. 213 zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 214 zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 215 ! thickness of water column between surface and level k at u/v point 216 zdepu = 0.5_wp * ( ( fsdept (ji,jj,jk) + fsdept (ji+1,jj ,jk) ) & 217 - ( risfdep(ji,jj) + risfdep(ji+1,jj) ) - fse3u(ji,jj,miku(ji,jj)) ) 218 zdepv = 0.5_wp * ( ( fsdept (ji,jj,jk) + fsdept (ji,jj+1,jk) ) & 219 - ( risfdep(ji,jj) + risfdep(ji,jj+1) ) - fse3v(ji,jj,mikv(ji,jj)) ) 220 ! 221 zwz(ji,jj,jk) = ( ( 1._wp - zfi) * zau / ( zbu - zeps ) & 222 & + zfi * zdepu * zslpml_hmlpu(ji,jj) ) * umask(ji,jj,jk) 223 zww(ji,jj,jk) = ( ( 1._wp - zfj) * zav / ( zbv - zeps ) & 224 & + zfj * zdepv * zslpml_hmlpv(ji,jj) ) * vmask(ji,jj,jk) 233 225 !!gm modif to suppress omlmask.... (as in Griffies case) 234 ! 235 ! 236 ! 237 ! 238 ! 239 ! 240 ! 226 ! ! ! jk must be >= ML level for zf=1. otherwise zf=0. 227 ! zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp ) 228 ! zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp ) 229 ! zci = 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) ) 230 ! zcj = 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) ) 231 ! zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk) 232 ! zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk) 241 233 !!gm end modif 242 END DO 243 END DO 244 END DO 245 CALL lbc_lnk( zwz, 'U', -1. ) ; CALL lbc_lnk( zww, 'V', -1. ) ! lateral boundary conditions 246 ! 247 ! !* horizontal Shapiro filter 248 DO jk = 2, jpkm1 249 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 250 DO ji = 2, jpim1 251 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 252 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 253 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 254 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 255 & + 4.* zwz(ji ,jj ,jk) ) 256 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 257 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 258 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 259 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 260 & + 4.* zww(ji,jj ,jk) ) 261 END DO 262 END DO 263 DO jj = 3, jpj-2 ! other rows 264 DO ji = fs_2, fs_jpim1 ! vector opt. 265 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 266 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 267 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 268 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 269 & + 4.* zwz(ji ,jj ,jk) ) 270 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 271 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 272 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 273 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 274 & + 4.* zww(ji,jj ,jk) ) 275 END DO 276 END DO 277 ! !* decrease along coastal boundaries 278 DO jj = 2, jpjm1 279 DO ji = fs_2, fs_jpim1 ! vector opt. 280 uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk ) ) * 0.5_wp & 281 & * ( umask(ji,jj ,jk) + umask(ji,jj ,jk+1) ) * 0.5_wp & 282 & * umask(ji,jj,jk-1) 283 vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk ) ) * 0.5_wp & 284 & * ( vmask(ji ,jj,jk) + vmask(ji ,jj,jk+1) ) * 0.5_wp & 285 & * vmask(ji,jj,jk-1) 286 END DO 287 END DO 288 END DO 289 290 291 ! II. slopes at w point | wslpi = mij( d/di( prd ) / d/dz( prd ) 292 ! =========================== | wslpj = mij( d/dj( prd ) / d/dz( prd ) 293 ! 294 DO jk = 2, jpkm1 295 DO jj = 2, jpjm1 296 DO ji = fs_2, fs_jpim1 ! vector opt. 297 ! !* Local vertical density gradient evaluated from N^2 298 zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) * wmask(ji,jj,jk) 299 ! !* Slopes at w point 300 ! ! i- & j-gradient of density at w-points 301 zci = MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk ) & 302 & + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps ) * e1t(ji,jj) 303 zcj = MAX( vmask(ji,jj-1,jk ) + vmask(ji,jj,jk-1) & 304 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk ) , zeps ) * e2t(ji,jj) 305 zai = ( zgru (ji-1,jj,jk ) + zgru (ji,jj,jk-1) & 306 & + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk ) ) / zci 307 zaj = ( zgrv (ji,jj-1,jk ) + zgrv (ji,jj,jk-1) & 308 & + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk ) ) / zcj 309 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 310 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 311 zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zai ) ) 312 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zaj ) ) 313 ! ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 314 zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) ) ! zfk=1 in the ML otherwise zfk=0 315 zck = ( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj) - fsdepw(ji,jj,mikt(ji,jj)), 10._wp ) 316 zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) & 317 & + zck * wslpiml(ji,jj) * zfk ) * wmask(ji,jj,jk) 318 zww(ji,jj,jk) = ( zaj / ( zbj - zeps ) * ( 1._wp - zfk ) & 319 & + zck * wslpjml(ji,jj) * zfk ) * wmask(ji,jj,jk) 234 END DO 235 END DO 236 END DO 237 CALL lbc_lnk( zwz, 'U', -1. ) ; CALL lbc_lnk( zww, 'V', -1. ) ! lateral boundary conditions 238 ! 239 ! !* horizontal Shapiro filter 240 DO jk = 2, jpkm1 241 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 242 DO ji = 2, jpim1 243 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 244 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 245 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 246 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 247 & + 4.* zwz(ji ,jj ,jk) ) 248 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 249 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 250 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 251 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 252 & + 4.* zww(ji,jj ,jk) ) 253 END DO 254 END DO 255 DO jj = 3, jpj-2 ! other rows 256 DO ji = fs_2, fs_jpim1 ! vector opt. 257 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 258 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 259 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 260 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 261 & + 4.* zwz(ji ,jj ,jk) ) 262 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 263 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 264 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 265 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 266 & + 4.* zww(ji,jj ,jk) ) 267 END DO 268 END DO 269 ! !* decrease along coastal boundaries 270 DO jj = 2, jpjm1 271 DO ji = fs_2, fs_jpim1 ! vector opt. 272 uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk ) ) * 0.5_wp & 273 & * ( umask(ji,jj ,jk) + umask(ji,jj ,jk+1) ) * 0.5_wp 274 vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk ) ) * 0.5_wp & 275 & * ( vmask(ji ,jj,jk) + vmask(ji ,jj,jk+1) ) * 0.5_wp 276 END DO 277 END DO 278 END DO 279 280 281 ! II. slopes at w point | wslpi = mij( d/di( prd ) / d/dz( prd ) 282 ! =========================== | wslpj = mij( d/dj( prd ) / d/dz( prd ) 283 ! 284 DO jk = 2, jpkm1 285 DO jj = 2, jpjm1 286 DO ji = fs_2, fs_jpim1 ! vector opt. 287 ! !* Local vertical density gradient evaluated from N^2 288 zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) 289 ! !* Slopes at w point 290 ! ! i- & j-gradient of density at w-points 291 zci = MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk ) & 292 & + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps ) * e1t(ji,jj) 293 zcj = MAX( vmask(ji,jj-1,jk ) + vmask(ji,jj,jk-1) & 294 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk ) , zeps ) * e2t(ji,jj) 295 zai = ( zgru (ji-1,jj,jk ) + zgru (ji,jj,jk-1) & 296 & + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk ) ) / zci * wmask (ji,jj,jk) 297 zaj = ( zgrv (ji,jj-1,jk ) + zgrv (ji,jj,jk-1) & 298 & + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk ) ) / zcj * wmask (ji,jj,jk) 299 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 300 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 301 zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zai ) ) 302 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zaj ) ) 303 ! ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 304 zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) ) ! zfk=1 in the ML otherwise zfk=0 305 zck = ( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj) - fsdepw(ji,jj,mikt(ji,jj)), 10._wp ) 306 zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk ) * wmask(ji,jj,jk) 307 zww(ji,jj,jk) = ( zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk ) * wmask(ji,jj,jk) 320 308 321 309 !!gm modif to suppress omlmask.... (as in Griffies operator) 322 ! 323 ! 324 ! 325 ! 326 ! 310 ! ! ! jk must be >= ML level for zfk=1. otherwise zfk=0. 311 ! zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp ) 312 ! zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10. ) 313 ! zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk) 314 ! zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk) 327 315 !!gm end modif 328 END DO 329 END DO 330 END DO 331 CALL lbc_lnk( zwz, 'T', -1. ) ; CALL lbc_lnk( zww, 'T', -1. ) ! lateral boundary conditions 332 ! 333 ! !* horizontal Shapiro filter 334 DO jk = 2, jpkm1 335 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 336 DO ji = 2, jpim1 337 zcofw = tmask(ji,jj,jk) * z1_16 338 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 339 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 340 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 341 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 342 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 343 344 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 345 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 346 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 347 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 348 & + 4.* zww(ji ,jj ,jk) ) * zcofw 349 END DO 350 END DO 351 DO jj = 3, jpj-2 ! other rows 352 DO ji = fs_2, fs_jpim1 ! vector opt. 353 zcofw = tmask(ji,jj,jk) * z1_16 354 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 355 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 356 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 357 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 358 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 359 360 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 361 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 362 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 363 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 364 & + 4.* zww(ji ,jj ,jk) ) * zcofw 365 END DO 366 END DO 367 ! !* decrease along coastal boundaries 368 DO jj = 2, jpjm1 369 DO ji = fs_2, fs_jpim1 ! vector opt. 370 zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) & 371 & * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 372 wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck * wmask(ji,jj,jk) 373 wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck * wmask(ji,jj,jk) 374 END DO 375 END DO 376 END DO 377 378 ! III. Specific grid points 379 ! =========================== 380 ! 381 IF( cp_cfg == "orca" .AND. jp_cfg == 4 ) THEN ! ORCA_R4 configuration: horizontal diffusion in specific area 382 ! ! Gibraltar Strait 383 ij0 = 50 ; ij1 = 53 384 ii0 = 69 ; ii1 = 71 ; uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 385 ij0 = 51 ; ij1 = 53 386 ii0 = 68 ; ii1 = 71 ; vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 387 ii0 = 69 ; ii1 = 71 ; wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 388 ii0 = 69 ; ii1 = 71 ; wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 389 ! 390 ! ! Mediterrannean Sea 391 ij0 = 49 ; ij1 = 56 392 ii0 = 71 ; ii1 = 90 ; uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 393 ij0 = 50 ; ij1 = 56 394 ii0 = 70 ; ii1 = 90 ; vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 395 ii0 = 71 ; ii1 = 90 ; wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 396 ii0 = 71 ; ii1 = 90 ; wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 397 ENDIF 398 399 400 ! IV. Lateral boundary conditions 401 ! =============================== 402 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) 403 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 404 405 406 IF(ln_ctl) THEN 407 CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk) 408 CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk) 409 ENDIF 410 ! 411 412 ELSEIF ( lk_vvl ) THEN 413 414 IF(lwp) THEN 415 WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 416 ENDIF 417 418 ! geopotential diffusion in s-coordinates on tracers and/or momentum 419 ! The slopes of s-surfaces are computed at each time step due to vvl 420 ! The slopes for momentum diffusion are i- or j- averaged of those on tracers 421 422 ! set the slope of diffusion to the slope of s-surfaces 423 ! ( c a u t i o n : minus sign as fsdep has positive value ) 424 DO jj = 2, jpjm1 425 DO ji = fs_2, fs_jpim1 ! vector opt. 426 uslp(ji,jj,1) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,1) - fsdept_b(ji ,jj ,1) ) * umask(ji,jj,1) 427 vslp(ji,jj,1) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,1) - fsdept_b(ji ,jj ,1) ) * vmask(ji,jj,1) 428 wslpi(ji,jj,1) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,1) - fsdepw_b(ji-1,jj,1) ) * tmask(ji,jj,1) * 0.5 429 wslpj(ji,jj,1) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,1) - fsdepw_b(ji,jj-1,1) ) * tmask(ji,jj,1) * 0.5 430 END DO 431 END DO 432 433 DO jk = 2, jpk 434 DO jj = 2, jpjm1 435 DO ji = fs_2, fs_jpim1 ! vector opt. 436 uslp(ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk) 437 vslp(ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk) 438 wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) & 439 & * wmask(ji,jj,jk) * 0.5 440 wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) & 441 & * wmask(ji,jj,jk) * 0.5 442 END DO 443 END DO 444 END DO 445 446 ! Lateral boundary conditions on the slopes 447 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) 448 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 449 450 if( kt == nit000 ) then 451 IF(lwp) WRITE(numout,*) ' max slop: u',SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)), & 452 & ' wi', sqrt(MAXVAL(wslpi)), ' wj', sqrt(MAXVAL(wslpj)) 453 endif 454 455 IF(ln_ctl) THEN 456 CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk) 457 CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk) 458 ENDIF 459 316 END DO 317 END DO 318 END DO 319 CALL lbc_lnk( zwz, 'T', -1. ) ; CALL lbc_lnk( zww, 'T', -1. ) ! lateral boundary conditions 320 ! 321 ! !* horizontal Shapiro filter 322 DO jk = 2, jpkm1 323 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 324 DO ji = 2, jpim1 325 zcofw = wmask(ji,jj,jk) * z1_16 326 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 327 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 328 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 329 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 330 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 331 332 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 333 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 334 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 335 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 336 & + 4.* zww(ji ,jj ,jk) ) * zcofw 337 END DO 338 END DO 339 DO jj = 3, jpj-2 ! other rows 340 DO ji = fs_2, fs_jpim1 ! vector opt. 341 zcofw = wmask(ji,jj,jk) * z1_16 342 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 343 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 344 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 345 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 346 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 347 348 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 349 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 350 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 351 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 352 & + 4.* zww(ji ,jj ,jk) ) * zcofw 353 END DO 354 END DO 355 ! !* decrease in vicinity of topography 356 DO jj = 2, jpjm1 357 DO ji = fs_2, fs_jpim1 ! vector opt. 358 zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) & 359 & * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 360 wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck 361 wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck 362 END DO 363 END DO 364 END DO 365 366 ! IV. Lateral boundary conditions 367 ! =============================== 368 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) 369 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 370 371 IF(ln_ctl) THEN 372 CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk) 373 CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk) 460 374 ENDIF 461 375 ! 462 376 CALL wrk_dealloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv ) 463 CALL wrk_dealloc( jpi,jpj, zhmlpu, zhmlpv)377 CALL wrk_dealloc( jpi,jpj, zslpml_hmlpu, zslpml_hmlpv ) 464 378 ! 465 379 IF( nn_timing == 1 ) CALL timing_stop('ldf_slp') … … 468 382 469 383 470 SUBROUTINE ldf_slp_ grif( kt )471 !!---------------------------------------------------------------------- 472 !! *** ROUTINE ldf_slp_ grif***384 SUBROUTINE ldf_slp_triad ( kt ) 385 !!---------------------------------------------------------------------- 386 !! *** ROUTINE ldf_slp_triad *** 473 387 !! 474 388 !! ** Purpose : Compute the squared slopes of neutral surfaces (slope 475 !! of iso-pycnal surfaces referenced locally) (ln_traldf_ grif=T)389 !! of iso-pycnal surfaces referenced locally) (ln_traldf_triad=T) 476 390 !! at W-points using the Griffies quarter-cells. 477 391 !! … … 488 402 REAL(wp) :: zfacti, zfactj ! local scalars 489 403 REAL(wp) :: znot_thru_surface ! local scalars 490 REAL(wp) :: zdit, zdis, zdjt, zdjs, zdkt, zdks, zbu, zbv, zbti, zbtj 404 REAL(wp) :: zdit, zdis, zdkt, zbu, zbti, zisw 405 REAL(wp) :: zdjt, zdjs, zdks, zbv, zbtj, zjsw 491 406 REAL(wp) :: zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_g_raw, zti_g_lim 492 407 REAL(wp) :: zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_g_raw, ztj_g_lim 493 408 REAL(wp) :: zdzrho_raw 409 REAL(wp) :: zbeta0, ze3_e1, ze3_e2 494 410 REAL(wp), POINTER, DIMENSION(:,:) :: z1_mlbw 411 REAL(wp), POINTER, DIMENSION(:,:,:) :: zalbet 495 412 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zdxrho , zdyrho, zdzrho ! Horizontal and vertical density gradients 496 413 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zti_mlb, ztj_mlb ! for Griffies operator only 497 414 !!---------------------------------------------------------------------- 498 415 ! 499 IF( nn_timing == 1 ) CALL timing_start('ldf_slp_ grif')416 IF( nn_timing == 1 ) CALL timing_start('ldf_slp_triad') 500 417 ! 501 418 CALL wrk_alloc( jpi,jpj, z1_mlbw ) 419 CALL wrk_alloc( jpi,jpj,jpk, zalbet ) 502 420 CALL wrk_alloc( jpi,jpj,jpk,2, zdxrho , zdyrho, zdzrho, klstart = 0 ) 503 421 CALL wrk_alloc( jpi,jpj, 2,2, zti_mlb, ztj_mlb, kkstart = 0, klstart = 0 ) … … 517 435 zdjt = ( tsb(ji,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! j-gradient of T & S at v-point 518 436 zdjs = ( tsb(ji,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 519 zdxrho_raw = ( - rab_b(ji+ip,jj ,jk,jp_tem) * zdit + rab_b(ji+ip,jj ,jk,jp_sal) * zdis ) /e1u(ji,jj)520 zdyrho_raw = ( - rab_b(ji ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji ,jj+jp,jk,jp_sal) * zdjs ) /e2v(ji,jj)521 zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw) ! keep the sign522 zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw)437 zdxrho_raw = ( - rab_b(ji+ip,jj ,jk,jp_tem) * zdit + rab_b(ji+ip,jj ,jk,jp_sal) * zdis ) * r1_e1u(ji,jj) 438 zdyrho_raw = ( - rab_b(ji ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji ,jj+jp,jk,jp_sal) * zdjs ) * r1_e2v(ji,jj) 439 zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 440 zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 523 441 END DO 524 442 END DO … … 531 449 zdit = gtsu(ji,jj,jp_tem) ; zdjt = gtsv(ji,jj,jp_tem) ! i- & j-gradient of Temperature 532 450 zdis = gtsu(ji,jj,jp_sal) ; zdjs = gtsv(ji,jj,jp_sal) ! i- & j-gradient of Salinity 533 zdxrho_raw = ( - rab_b(ji+ip,jj ,iku,jp_tem) * zdit + rab_b(ji+ip,jj ,iku,jp_sal) * zdis ) /e1u(ji,jj)534 zdyrho_raw = ( - rab_b(ji ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji ,jj+jp,ikv,jp_sal) * zdjs ) /e2v(ji,jj)451 zdxrho_raw = ( - rab_b(ji+ip,jj ,iku,jp_tem) * zdit + rab_b(ji+ip,jj ,iku,jp_sal) * zdis ) * r1_e1u(ji,jj) 452 zdyrho_raw = ( - rab_b(ji ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji ,jj+jp,ikv,jp_sal) * zdjs ) * r1_e2v(ji,jj) 535 453 zdxrho(ji+ip,jj ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 536 454 zdyrho(ji ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) … … 552 470 zdks = 0._wp 553 471 ENDIF 554 zdzrho_raw = ( - rab_b(ji,jj,jk,jp_tem) * zdkt + rab_b(ji,jj,jk,jp_sal) *zdks ) / fse3w(ji,jj,jk+kp)555 zdzrho(ji,jj,jk,kp) = - MIN( - repsln , zdzrho_raw ) ! force zdzrho >= repsln472 zdzrho_raw = ( - zalbet(ji,jj,jk) * zdkt + zbeta0*zdks ) / fse3w(ji,jj,jk+kp) 473 zdzrho(ji,jj,jk,kp) = - MIN( - repsln , zdzrho_raw ) ! force zdzrho >= repsln 556 474 END DO 557 475 END DO … … 586 504 ! 587 505 jk = nmln(ji+ip,jj) + 1 588 IF( jk > mbkt(ji+ip,jj) ) THEN !ML reaches bottom 589 zti_mlb(ji+ip,jj ,1-ip,kp) = 0.0_wp 590 ELSE 591 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 592 zti_g_raw = ( zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp) & 593 & - ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) / e1u(ji,jj) ) * umask(ji,jj,jk) 594 zti_mlb(ji+ip,jj ,1-ip,kp) = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 506 IF( jk > mbkt(ji+ip,jj) ) THEN ! ML reaches bottom 507 zti_mlb(ji+ip,jj ,1-ip,kp) = 0.0_wp 508 ELSE 509 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 510 zti_g_raw = ( zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp) & 511 & - ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) * r1_e1u(ji,jj) ) * umask(ji,jj,jk) 512 ze3_e1 = fse3w(ji+ip,jj,jk-kp) * r1_e1u(ji,jj) 513 zti_mlb(ji+ip,jj ,1-ip,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1 , ABS( zti_g_raw ) ), zti_g_raw ) 595 514 ENDIF 596 515 ! 597 516 jk = nmln(ji,jj+jp) + 1 598 517 IF( jk > mbkt(ji,jj+jp) ) THEN !ML reaches bottom 599 ztj_mlb(ji ,jj+jp,1-jp,kp) = 0.0_wp518 ztj_mlb(ji ,jj+jp,1-jp,kp) = 0.0_wp 600 519 ELSE 601 ztj_g_raw = ( zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp) & 602 & - ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) 603 ztj_mlb(ji ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) 520 ztj_g_raw = ( zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp) & 521 & - ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) 522 ze3_e2 = fse3w(ji,jj+jp,jk-kp) / e2v(ji,jj) 523 ztj_mlb(ji ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2 , ABS( ztj_g_raw ) ), ztj_g_raw ) 604 524 ENDIF 605 525 END DO … … 630 550 zti_raw = zdxrho(ji+ip,jj ,jk,1-ip) / zdzrho(ji+ip,jj ,jk,kp) ! unmasked 631 551 ztj_raw = zdyrho(ji ,jj+jp,jk,1-jp) / zdzrho(ji ,jj+jp,jk,kp) 632 552 ! 633 553 ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface 634 zti_coord = znot_thru_surface * ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) /e1u(ji,jj)635 ztj_coord = znot_thru_surface * ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj)! unmasked554 zti_coord = znot_thru_surface * ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) * r1_e1u(ji,jj) 555 ztj_coord = znot_thru_surface * ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) * r1_e2v(ji,jj) ! unmasked 636 556 zti_g_raw = zti_raw - zti_coord ! ref to geopot surfaces 637 557 ztj_g_raw = ztj_raw - ztj_coord 638 zti_g_lim = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 639 ztj_g_lim = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) 558 ! additional limit required in bilaplacian case 559 ze3_e1 = fse3w(ji+ip,jj ,jk+kp) * r1_e1u(ji,jj) 560 ze3_e2 = fse3w(ji ,jj+jp,jk+kp) * r1_e2v(ji,jj) 561 ! NB: hard coded factor 5 (can be a namelist parameter...) 562 zti_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1, ABS( zti_g_raw ) ), zti_g_raw ) 563 ztj_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2, ABS( ztj_g_raw ) ), ztj_g_raw ) 640 564 ! 641 565 ! Below ML use limited zti_g as is & mask … … 666 590 ! 667 591 IF( ln_triad_iso ) THEN 668 zti_raw = zti_lim* *2/ zti_raw669 ztj_raw = ztj_lim* *2/ ztj_raw592 zti_raw = zti_lim*zti_lim / zti_raw 593 ztj_raw = ztj_lim*ztj_lim / ztj_raw 670 594 zti_raw = SIGN( MIN( ABS(zti_lim), ABS( zti_raw ) ), zti_raw ) 671 595 ztj_raw = SIGN( MIN( ABS(ztj_lim), ABS( ztj_raw ) ), ztj_raw ) 672 zti_lim = zfacti * zti_lim & 673 & + ( 1._wp - zfacti ) * zti_raw 674 ztj_lim = zfactj * ztj_lim & 675 & + ( 1._wp - zfactj ) * ztj_raw 596 zti_lim = zfacti * zti_lim + ( 1._wp - zfacti ) * zti_raw 597 ztj_lim = zfactj * ztj_lim + ( 1._wp - zfactj ) * ztj_raw 676 598 ENDIF 677 triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim 678 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim 679 ! 680 zbu = e1u(ji ,jj) * e2u(ji ,jj) * fse3u(ji ,jj,jk ) 681 zbv = e1v(ji ,jj) * e2v(ji ,jj) * fse3v(ji ,jj,jk ) 682 zbti = e1t(ji+ip,jj) * e2t(ji+ip,jj) * fse3w(ji+ip,jj,jk+kp) 683 zbtj = e1t(ji,jj+jp) * e2t(ji,jj+jp) * fse3w(ji,jj+jp,jk+kp) 684 ! 685 !!gm this may inhibit vectorization on Vect Computers, and even on scalar computers.... ==> to be checked 686 wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim**2 ! masked 687 wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_g_lim**2 599 ! ! switching triad scheme 600 zisw = (rn_sw_triad - 1._wp ) + rn_sw_triad & 601 & * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - ip ) * SIGN( 1._wp , zdxrho(ji+ip,jj,jk,1-ip) ) ) 602 zjsw = (rn_sw_triad - 1._wp ) + rn_sw_triad & 603 & * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - jp ) * SIGN( 1._wp , zdyrho(ji,jj+jp,jk,1-jp) ) ) 604 ! 605 triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim * zisw 606 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim * zjsw 607 ! 608 zbu = e1e2u(ji ,jj ) * fse3u(ji ,jj ,jk ) 609 zbv = e1e2v(ji ,jj ) * fse3v(ji ,jj ,jk ) 610 zbti = e1e2t(ji+ip,jj ) * fse3w(ji+ip,jj ,jk+kp) 611 zbtj = e1e2t(ji ,jj+jp) * fse3w(ji ,jj+jp,jk+kp) 612 ! 613 wslp2(ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim*zti_g_lim ! masked 614 wslp2(ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_g_lim*ztj_g_lim 688 615 END DO 689 616 END DO … … 697 624 ! 698 625 CALL wrk_dealloc( jpi,jpj, z1_mlbw ) 626 CALL wrk_dealloc( jpi,jpj,jpk, zalbet ) 699 627 CALL wrk_dealloc( jpi,jpj,jpk,2, zdxrho , zdyrho, zdzrho, klstart = 0 ) 700 628 CALL wrk_dealloc( jpi,jpj, 2,2, zti_mlb, ztj_mlb, kkstart = 0, klstart = 0 ) 701 629 ! 702 IF( nn_timing == 1 ) CALL timing_stop('ldf_slp_ grif')703 ! 704 END SUBROUTINE ldf_slp_ grif630 IF( nn_timing == 1 ) CALL timing_stop('ldf_slp_triad') 631 ! 632 END SUBROUTINE ldf_slp_triad 705 633 706 634 … … 728 656 INTEGER :: ji , jj , jk ! dummy loop indices 729 657 INTEGER :: iku, ikv, ik, ikm1 ! local integers 730 REAL(wp) :: zeps, zm1_g, zm1_2g 658 REAL(wp) :: zeps, zm1_g, zm1_2g, z1_slpmax ! local scalars 731 659 REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - - 732 660 REAL(wp) :: zcj, zfj, zav, zbv, zaj, zbj ! - - … … 739 667 zm1_g = -1.0_wp / grav 740 668 zm1_2g = -0.5_wp / grav 669 z1_slpmax = 1._wp / rn_slpmax 741 670 ! 742 671 uslpml (1,:) = 0._wp ; uslpml (jpi,:) = 0._wp … … 750 679 DO ji = 1, jpi 751 680 ik = nmln(ji,jj) - 1 752 IF( jk <= ik ) THEN ;omlmask(ji,jj,jk) = 1._wp753 ELSE ;omlmask(ji,jj,jk) = 0._wp681 IF( jk <= ik ) THEN ; omlmask(ji,jj,jk) = 1._wp 682 ELSE ; omlmask(ji,jj,jk) = 0._wp 754 683 ENDIF 755 684 END DO … … 773 702 ! 774 703 ! !- vertical density gradient for u- and v-slopes (from dzr at T-point) 775 iku = MIN( MAX( nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1 ) ! ML (MAX of T-pts, bound by jpkm1)776 ikv = MIN( MAX( nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) !704 iku = MIN( MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1 ) ! ML (MAX of T-pts, bound by jpkm1) 705 ikv = MIN( MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) ! 777 706 zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj ,iku) ) 778 707 zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji ,jj+1,ikv) ) 779 708 ! !- horizontal density gradient at u- & v-points 780 zau = p_gru(ji,jj,iku) /e1u(ji,jj)781 zav = p_grv(ji,jj,ikv) /e2v(ji,jj)709 zau = p_gru(ji,jj,iku) * r1_e1u(ji,jj) 710 zav = p_grv(ji,jj,ikv) * r1_e2v(ji,jj) 782 711 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 783 712 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 784 zbu = MIN( zbu , - 100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau ) )785 zbv = MIN( zbv , - 100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav ) )713 zbu = MIN( zbu , - z1_slpmax * ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau ) ) 714 zbv = MIN( zbv , - z1_slpmax * ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav ) ) 786 715 ! !- Slope at u- & v-points (uslpml, vslpml) 787 716 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku) … … 808 737 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zaj ) ) 809 738 ! !- i- & j-slope at w-points (wslpiml, wslpjml) 810 wslpiml(ji,jj) = zai / ( zbi - zeps ) * wmask (ji,jj,ik)811 wslpjml(ji,jj) = zaj / ( zbj - zeps ) * wmask (ji,jj,ik)739 wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) 740 wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik) 812 741 END DO 813 742 END DO … … 827 756 !! ** Purpose : Initialization for the isopycnal slopes computation 828 757 !! 829 !! ** Method : read the nammbf namelist and check the parameter 830 !! values called by tra_dmp at the first timestep (nit000) 758 !! ** Method : 831 759 !!---------------------------------------------------------------------- 832 760 INTEGER :: ji, jj, jk ! dummy loop indices … … 841 769 WRITE(numout,*) '~~~~~~~~~~~~' 842 770 ENDIF 843 844 IF( ln_traldf_grif ) THEN ! Griffies operator : triad of slopes 845 ALLOCATE( triadi_g(jpi,jpj,jpk,0:1,0:1) , triadj_g(jpi,jpj,jpk,0:1,0:1) , wslp2(jpi,jpj,jpk) , STAT=ierr ) 846 ALLOCATE( triadi (jpi,jpj,jpk,0:1,0:1) , triadj (jpi,jpj,jpk,0:1,0:1) , STAT=ierr ) 847 IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Griffies operator slope' ) 848 ! 771 ! 772 ALLOCATE( ah_wslp2(jpi,jpj,jpk) , akz(jpi,jpj,jpk) , STAT=ierr ) 773 IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate ah_slp2 or akz' ) 774 ! 775 IF( ln_traldf_triad ) THEN ! Griffies operator : triad of slopes 776 IF(lwp) WRITE(numout,*) ' Griffies (triad) operator initialisation' 777 ALLOCATE( triadi_g(jpi,jpj,jpk,0:1,0:1) , triadj_g(jpi,jpj,jpk,0:1,0:1) , & 778 & triadi (jpi,jpj,jpk,0:1,0:1) , triadj (jpi,jpj,jpk,0:1,0:1) , & 779 & wslp2 (jpi,jpj,jpk) , STAT=ierr ) 780 IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Griffies operator slope' ) 849 781 IF( ln_dynldf_iso ) CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' ) 850 782 ! 851 783 ELSE ! Madec operator : slopes at u-, v-, and w-points 852 ALLOCATE( uslp(jpi,jpj,jpk) , vslp(jpi,jpj,jpk) , wslpi(jpi,jpj,jpk) , wslpj(jpi,jpj,jpk) , & 853 & omlmask(jpi,jpj,jpk) , uslpml(jpi,jpj) , vslpml(jpi,jpj) , wslpiml(jpi,jpj) , wslpjml(jpi,jpj) , STAT=ierr ) 784 IF(lwp) WRITE(numout,*) ' Madec operator initialisation' 785 ALLOCATE( omlmask(jpi,jpj,jpk) , & 786 & uslp(jpi,jpj,jpk) , uslpml(jpi,jpj) , wslpi(jpi,jpj,jpk) , wslpiml(jpi,jpj) , & 787 & vslp(jpi,jpj,jpk) , vslpml(jpi,jpj) , wslpj(jpi,jpj,jpk) , wslpjml(jpi,jpj) , STAT=ierr ) 854 788 IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Madec operator slope ' ) 855 789 … … 861 795 wslpj(:,:,:) = 0._wp ; wslpjml(:,:) = 0._wp 862 796 863 IF(ln_sco .AND. (ln_traldf_hor .OR. ln_dynldf_hor )) THEN 864 IF(lwp) WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 865 866 ! geopotential diffusion in s-coordinates on tracers and/or momentum 867 ! The slopes of s-surfaces are computed once (no call to ldfslp in step) 868 ! The slopes for momentum diffusion are i- or j- averaged of those on tracers 869 870 ! set the slope of diffusion to the slope of s-surfaces 871 ! ( c a u t i o n : minus sign as fsdep has positive value ) 872 DO jk = 1, jpk 873 DO jj = 2, jpjm1 874 DO ji = fs_2, fs_jpim1 ! vector opt. 875 uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk) 876 vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk) 877 wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5 878 wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5 879 END DO 880 END DO 881 END DO 882 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) ! Lateral boundary conditions 883 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 884 ENDIF 797 !!gm I no longer understand this..... 798 !!gm IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN 799 ! IF(lwp) WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 800 ! 801 ! ! geopotential diffusion in s-coordinates on tracers and/or momentum 802 ! ! The slopes of s-surfaces are computed once (no call to ldfslp in step) 803 ! ! The slopes for momentum diffusion are i- or j- averaged of those on tracers 804 ! 805 ! ! set the slope of diffusion to the slope of s-surfaces 806 ! ! ( c a u t i o n : minus sign as fsdep has positive value ) 807 ! DO jk = 1, jpk 808 ! DO jj = 2, jpjm1 809 ! DO ji = fs_2, fs_jpim1 ! vector opt. 810 ! uslp (ji,jj,jk) = - ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 811 ! vslp (ji,jj,jk) = - ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 812 ! wslpi(ji,jj,jk) = - ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * r1_e1t(ji,jj) * wmask(ji,jj,jk) * 0.5 813 ! wslpj(ji,jj,jk) = - ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * r1_e2t(ji,jj) * wmask(ji,jj,jk) * 0.5 814 ! END DO 815 ! END DO 816 ! END DO 817 ! CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) ! Lateral boundary conditions 818 ! CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 819 !!gm ENDIF 885 820 ENDIF 886 821 ! … … 888 823 ! 889 824 END SUBROUTINE ldf_slp_init 890 891 #else892 !!------------------------------------------------------------------------893 !! Dummy module : NO Rotation of lateral mixing tensor894 !!------------------------------------------------------------------------895 LOGICAL, PUBLIC, PARAMETER :: lk_ldfslp = .FALSE. !: slopes flag896 CONTAINS897 SUBROUTINE ldf_slp( kt, prd, pn2 ) ! Dummy routine898 INTEGER, INTENT(in) :: kt899 REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2900 WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1)901 END SUBROUTINE ldf_slp902 SUBROUTINE ldf_slp_grif( kt ) ! Dummy routine903 INTEGER, INTENT(in) :: kt904 WRITE(*,*) 'ldf_slp_grif: You should not have seen this print! error?', kt905 END SUBROUTINE ldf_slp_grif906 SUBROUTINE ldf_slp_init ! Dummy routine907 END SUBROUTINE ldf_slp_init908 #endif909 825 910 826 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
