- Timestamp:
- 2015-09-24T08:31:40+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r5737 r5758 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 !!gm 25 ! USE ldfdyn ! lateral diffusion: eddy viscosity coef. 26 !!gm to be removed 27 USE ldfdyn_oce ! lateral diffusion: eddy viscosity coef. 28 !!gm 27 29 USE phycst ! physical constants 28 30 USE zdfmxl ! mixed layer depth … … 30 32 ! 31 33 USE in_out_manager ! I/O manager 34 USE prtctl ! Print control 32 35 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 33 USE prtctl ! Print control 36 USE lib_mpp ! distribued memory computing library 37 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 34 38 USE wrk_nemo ! work arrays 35 39 USE timing ! Timing 36 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)37 40 38 41 IMPLICIT NONE 39 42 PRIVATE 40 43 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 44 PUBLIC ldf_slp ! routine called by step.F90 45 PUBLIC ldf_slp_triad ! routine called by step.F90 46 PUBLIC ldf_slp_init ! routine called by nemogcm.F90 47 48 LOGICAL , PUBLIC :: l_ldfslp = .FALSE. !: slopes flag 49 50 LOGICAL , PUBLIC :: ln_traldf_iso = .TRUE. !: iso-neutral direction 51 LOGICAL , PUBLIC :: ln_traldf_triad = .FALSE. !: griffies triad scheme 52 53 LOGICAL , PUBLIC :: ln_triad_iso = .FALSE. !: pure horizontal mixing in ML 54 LOGICAL , PUBLIC :: ln_botmix_triad = .FALSE. !: mixing on bottom 55 REAL(wp), PUBLIC :: rn_sw_triad = 1._wp !: =1 switching triads ; =0 all four triads used 56 REAL(wp), PUBLIC :: rn_slpmax = 0.01_wp !: slope limit 57 58 LOGICAL , PUBLIC :: l_grad_zps = .FALSE. !: special treatment for Horz Tgradients w partial steps (triad operator) 59 60 ! !! Classic operator (Madec) 48 61 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: uslp, wslpi !: i_slope at U- and W-points 49 62 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: vslp, wslpj !: j-slope at V- and W-points 50 ! !! Griffies operator63 ! !! triad operator (Griffies) 51 64 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wslp2 !: wslp**2 from Griffies quarter cells 52 65 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) :: triadi_g, triadj_g !: skew flux slopes relative to geopotentials 53 66 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 67 ! !! both operators 68 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ah_wslp2 !: ah * slope^2 at w-point 69 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akz !: stabilizing vertical diffusivity 70 71 ! !! Madec operator 57 72 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: omlmask ! mask of the surface mixed layer at T-pt 58 73 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: uslpml, wslpiml ! i_slope at U- and W-points just below the mixed layer … … 63 78 !! * Substitutions 64 79 # include "domzgr_substitute.h90" 65 # include "ldftra_substitute.h90"66 # include "ldfeiv_substitute.h90"67 80 # include "vectopt_loop_substitute.h90" 68 81 !!---------------------------------------------------------------------- 69 !! NEMO/OPA 4.0 , NEMO Consortium (201 1)82 !! NEMO/OPA 4.0 , NEMO Consortium (2014) 70 83 !! $Id$ 71 84 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 105 118 INTEGER :: ii0, ii1, iku ! temporary integer 106 119 INTEGER :: ij0, ij1, ikv ! temporary integer 107 REAL(wp) :: zeps, zm1_g, zm1_2g, z1_16, zcofw ! local scalars120 REAL(wp) :: zeps, zm1_g, zm1_2g, z1_16, zcofw, z1_slpmax ! local scalars 108 121 REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - - 109 122 REAL(wp) :: zcj, zfj, zav, zbv, zaj, zbj ! - - 110 123 REAL(wp) :: zck, zfk, zbw ! - - 111 REAL(wp) :: zdepv, zdepu ! - -112 124 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz, zww 113 125 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdzr 114 126 REAL(wp), POINTER, DIMENSION(:,:,:) :: zgru, zgrv 115 REAL(wp), POINTER, DIMENSION(:,: ) :: zhmlpu, zhmlpv116 127 !!---------------------------------------------------------------------- 117 128 ! … … 119 130 ! 120 131 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 152 IF ( miku(ji,jj) > 1 ) zgru(ji,jj,miku(ji,jj)) = grui(ji,jj) 153 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 ! ! umask(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 IF (miku(ji,jj) .GT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji ,jj ), 5._wp) 191 IF (miku(ji,jj) .LT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji+1,jj ), 5._wp) 192 IF (miku(ji,jj) .EQ. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji ,jj ), hmlpt(ji+1,jj ), 5._wp) 193 IF (mikv(ji,jj) .GT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji ,jj ), 5._wp) 194 IF (mikv(ji,jj) .LT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji ,jj+1), 5._wp) 195 IF (mikv(ji,jj) .EQ. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji ,jj ), hmlpt(ji ,jj+1), 5._wp) 196 ENDDO 197 ENDDO 198 ELSE 199 DO jj = 2, jpjm1 200 DO ji = fs_2, fs_jpim1 ! vector opt. 201 zhmlpu(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji+1,jj ), 5._wp) 202 zhmlpv(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji ,jj+1), 5._wp) 203 ENDDO 204 ENDDO 205 END IF 206 DO jk = 2, jpkm1 !* Slopes at u and v points 207 DO jj = 2, jpjm1 208 DO ji = fs_2, fs_jpim1 ! vector opt. 209 ! ! horizontal and vertical density gradient at u- and v-points 210 zau = zgru(ji,jj,jk) * r1_e1u(ji,jj) 211 zav = zgrv(ji,jj,jk) * r1_e2v(ji,jj) 212 zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj ,jk) ) 213 zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji ,jj+1,jk) ) 214 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 215 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 216 zbu = MIN( zbu, -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,jk)* ABS( zau ) ) 217 zbv = MIN( zbv, -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav ) ) 218 ! ! uslp and vslp output in zwz and zww, resp. 219 zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj ,jk) ) 220 zfj = MAX( omlmask(ji,jj,jk), omlmask(ji ,jj+1,jk) ) 221 ! thickness of water column between surface and level k at u/v point 222 zdepu = 0.5_wp * ( ( fsdept(ji,jj,jk) + fsdept(ji+1,jj ,jk) ) & 223 - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj ) ) - fse3u(ji,jj,miku(ji,jj)) ) 224 zdepv = 0.5_wp * ( ( fsdept(ji,jj,jk) + fsdept(ji ,jj+1,jk) ) & 225 - 2 * MAX( risfdep(ji,jj), risfdep(ji ,jj+1) ) - fse3v(ji,jj,mikv(ji,jj)) ) 226 ! 227 zwz(ji,jj,jk) = ( 1. - zfi) * zau / ( zbu - zeps ) & 228 & + zfi * uslpml(ji,jj) * zdepu / zhmlpu(ji,jj) 229 zwz(ji,jj,jk) = zwz(ji,jj,jk) * wumask(ji,jj,jk) 230 zww(ji,jj,jk) = ( 1. - zfj) * zav / ( zbv - zeps ) & 231 & + zfj * vslpml(ji,jj) * zdepv / zhmlpv(ji,jj) 232 zww(ji,jj,jk) = zww(ji,jj,jk) * wvmask(ji,jj,jk) 233 234 132 133 zeps = 1.e-20_wp !== Local constant initialization ==! 134 z1_16 = 1.0_wp / 16._wp 135 zm1_g = -1.0_wp / grav 136 zm1_2g = -0.5_wp / grav 137 z1_slpmax = 1._wp / rn_slpmax 138 ! 139 zww(:,:,:) = 0._wp 140 zwz(:,:,:) = 0._wp 141 ! 142 DO jk = 1, jpk !== i- & j-gradient of density ==! 143 DO jj = 1, jpjm1 144 DO ji = 1, fs_jpim1 ! vector opt. 145 zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj ,jk) - prd(ji,jj,jk) ) 146 zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji ,jj+1,jk) - prd(ji,jj,jk) ) 147 END DO 148 END DO 149 END DO 150 IF( ln_zps ) THEN ! partial steps correction at the bottom ocean level 151 DO jj = 1, jpjm1 152 DO ji = 1, jpim1 153 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 154 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 155 END DO 156 END DO 157 ENDIF 158 ! 159 zdzr(:,:,1) = 0._wp !== Local vertical density gradient at T-point == ! (evaluated from N^2) 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 ! ! umask(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 ! 170 ! !== Slopes just below the mixed layer ==! 171 CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr ) ! output: uslpml, vslpml, wslpiml, wslpjml 172 173 174 ! I. slopes at u and v point | uslp = d/di( prd ) / d/dz( prd ) 175 ! =========================== | vslp = d/dj( prd ) / d/dz( prd ) 176 ! 177 DO jk = 2, jpkm1 !* Slopes at u and v points 178 DO jj = 2, jpjm1 179 DO ji = fs_2, fs_jpim1 ! vector opt. 180 ! ! horizontal and vertical density gradient at u- and v-points 181 zau = zgru(ji,jj,jk) * r1_e1u(ji,jj) 182 zav = zgrv(ji,jj,jk) * r1_e2v(ji,jj) 183 zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj ,jk) ) 184 zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji ,jj+1,jk) ) 185 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 186 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 187 zbu = MIN( zbu, - z1_slpmax * ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,jk)* ABS( zau ) ) 188 zbv = MIN( zbv, - z1_slpmax * ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav ) ) 189 ! ! uslp and vslp output in zwz and zww, resp. 190 zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 191 zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 192 zwz(ji,jj,jk) = ( ( 1. - zfi) * zau / ( zbu - zeps ) & 193 & + zfi * uslpml(ji,jj) & 194 & * 0.5_wp * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) ) & 195 & / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5._wp ) ) * umask(ji,jj,jk) 196 zww(ji,jj,jk) = ( ( 1. - zfj) * zav / ( zbv - zeps ) & 197 & + zfj * vslpml(ji,jj) & 198 & * 0.5_wp * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) ) & 199 & / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk) 235 200 !!gm modif to suppress omlmask.... (as in Griffies case) 236 ! 237 ! 238 ! 239 ! 240 ! 241 ! 242 ! 201 ! ! ! jk must be >= ML level for zf=1. otherwise zf=0. 202 ! zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp ) 203 ! zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp ) 204 ! zci = 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) ) 205 ! zcj = 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) ) 206 ! zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk) 207 ! zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk) 243 208 !!gm end modif 244 END DO 245 END DO 246 END DO 247 CALL lbc_lnk( zwz, 'U', -1. ) ; CALL lbc_lnk( zww, 'V', -1. ) ! lateral boundary conditions 248 ! 249 ! !* horizontal Shapiro filter 250 DO jk = 2, jpkm1 251 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 252 DO ji = 2, jpim1 253 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 254 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 255 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 256 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 257 & + 4.* zwz(ji ,jj ,jk) ) 258 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 259 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 260 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 261 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 262 & + 4.* zww(ji,jj ,jk) ) 263 END DO 264 END DO 265 DO jj = 3, jpj-2 ! other rows 266 DO ji = fs_2, fs_jpim1 ! vector opt. 267 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 268 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 269 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 270 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 271 & + 4.* zwz(ji ,jj ,jk) ) 272 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 273 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 274 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 275 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 276 & + 4.* zww(ji,jj ,jk) ) 277 END DO 278 END DO 279 ! !* decrease along coastal boundaries 280 DO jj = 2, jpjm1 281 DO ji = fs_2, fs_jpim1 ! vector opt. 282 uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk ) ) * 0.5_wp & 283 & * ( umask(ji,jj ,jk) + umask(ji,jj ,jk+1) ) * 0.5_wp & 284 & * umask(ji,jj,jk-1) 285 vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk ) ) * 0.5_wp & 286 & * ( vmask(ji ,jj,jk) + vmask(ji ,jj,jk+1) ) * 0.5_wp & 287 & * vmask(ji,jj,jk-1) 288 END DO 289 END DO 290 END DO 291 292 293 ! II. slopes at w point | wslpi = mij( d/di( prd ) / d/dz( prd ) 294 ! =========================== | wslpj = mij( d/dj( prd ) / d/dz( prd ) 295 ! 296 DO jk = 2, jpkm1 297 DO jj = 2, jpjm1 298 DO ji = fs_2, fs_jpim1 ! vector opt. 299 ! !* Local vertical density gradient evaluated from N^2 300 zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) * wmask(ji,jj,jk) 301 ! !* Slopes at w point 302 ! ! i- & j-gradient of density at w-points 303 zci = MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk ) & 304 & + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps ) * e1t(ji,jj) 305 zcj = MAX( vmask(ji,jj-1,jk ) + vmask(ji,jj,jk-1) & 306 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk ) , zeps ) * e2t(ji,jj) 307 zai = ( zgru (ji-1,jj,jk ) + zgru (ji,jj,jk-1) & 308 & + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk ) ) / zci 309 zaj = ( zgrv (ji,jj-1,jk ) + zgrv (ji,jj,jk-1) & 310 & + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk ) ) / zcj 311 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 312 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 313 zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zai ) ) 314 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zaj ) ) 315 ! ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 316 zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) ) ! zfk=1 in the ML otherwise zfk=0 317 zck = ( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj), 10._wp ) 318 zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) & 319 & + zck * wslpiml(ji,jj) * zfk ) * wmask(ji,jj,jk) 320 zww(ji,jj,jk) = ( zaj / ( zbj - zeps ) * ( 1._wp - zfk ) & 321 & + zck * wslpjml(ji,jj) * zfk ) * wmask(ji,jj,jk) 209 END DO 210 END DO 211 END DO 212 CALL lbc_lnk( zwz, 'U', -1. ) ; CALL lbc_lnk( zww, 'V', -1. ) ! lateral boundary conditions 213 ! 214 ! !* horizontal Shapiro filter 215 DO jk = 2, jpkm1 216 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 217 DO ji = 2, jpim1 218 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 219 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 220 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 221 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 222 & + 4.* zwz(ji ,jj ,jk) ) 223 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 224 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 225 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 226 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 227 & + 4.* zww(ji,jj ,jk) ) 228 END DO 229 END DO 230 DO jj = 3, jpj-2 ! other rows 231 DO ji = fs_2, fs_jpim1 ! vector opt. 232 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 233 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 234 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 235 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 236 & + 4.* zwz(ji ,jj ,jk) ) 237 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 238 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 239 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 240 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 241 & + 4.* zww(ji,jj ,jk) ) 242 END DO 243 END DO 244 ! !* decrease along coastal boundaries 245 DO jj = 2, jpjm1 246 DO ji = fs_2, fs_jpim1 ! vector opt. 247 uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk ) ) * 0.5_wp & 248 & * ( umask(ji,jj ,jk) + umask(ji,jj ,jk+1) ) * 0.5_wp 249 vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk ) ) * 0.5_wp & 250 & * ( vmask(ji ,jj,jk) + vmask(ji ,jj,jk+1) ) * 0.5_wp 251 END DO 252 END DO 253 END DO 254 255 256 ! II. slopes at w point | wslpi = mij( d/di( prd ) / d/dz( prd ) 257 ! =========================== | wslpj = mij( d/dj( prd ) / d/dz( prd ) 258 ! 259 DO jk = 2, jpkm1 260 DO jj = 2, jpjm1 261 DO ji = fs_2, fs_jpim1 ! vector opt. 262 ! !* Local vertical density gradient evaluated from N^2 263 zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) 264 ! !* Slopes at w point 265 ! ! i- & j-gradient of density at w-points 266 zci = MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk ) & 267 & + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps ) * e1t(ji,jj) 268 zcj = MAX( vmask(ji,jj-1,jk ) + vmask(ji,jj,jk-1) & 269 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk ) , zeps ) * e2t(ji,jj) 270 zai = ( zgru (ji-1,jj,jk ) + zgru (ji,jj,jk-1) & 271 & + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk ) ) / zci * tmask (ji,jj,jk) 272 zaj = ( zgrv (ji,jj-1,jk ) + zgrv (ji,jj,jk-1) & 273 & + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk ) ) / zcj * tmask (ji,jj,jk) 274 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 275 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 276 zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zai ) ) 277 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zaj ) ) 278 ! ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 279 zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) ) ! zfk=1 in the ML otherwise zfk=0 280 zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10._wp ) 281 zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk ) * tmask(ji,jj,jk) 282 zww(ji,jj,jk) = ( zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk ) * tmask(ji,jj,jk) 322 283 323 284 !!gm modif to suppress omlmask.... (as in Griffies operator) 324 ! 325 ! 326 ! 327 ! 328 ! 285 ! ! ! jk must be >= ML level for zfk=1. otherwise zfk=0. 286 ! zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp ) 287 ! zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10. ) 288 ! zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk) 289 ! zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk) 329 290 !!gm end modif 330 END DO 331 END DO 332 END DO 333 CALL lbc_lnk( zwz, 'T', -1. ) ; CALL lbc_lnk( zww, 'T', -1. ) ! lateral boundary conditions 334 ! 335 ! !* horizontal Shapiro filter 336 DO jk = 2, jpkm1 337 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 338 DO ji = 2, jpim1 339 zcofw = tmask(ji,jj,jk) * z1_16 340 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 341 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 342 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 343 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 344 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 345 346 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 347 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 348 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 349 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 350 & + 4.* zww(ji ,jj ,jk) ) * zcofw 351 END DO 352 END DO 353 DO jj = 3, jpj-2 ! other rows 354 DO ji = fs_2, fs_jpim1 ! vector opt. 355 zcofw = tmask(ji,jj,jk) * z1_16 356 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 357 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 358 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 359 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 360 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 361 362 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 363 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 364 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 365 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 366 & + 4.* zww(ji ,jj ,jk) ) * zcofw 367 END DO 368 END DO 369 ! !* decrease along coastal boundaries 370 DO jj = 2, jpjm1 371 DO ji = fs_2, fs_jpim1 ! vector opt. 372 zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) & 373 & * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 374 wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck * wmask(ji,jj,jk) 375 wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck * wmask(ji,jj,jk) 376 END DO 377 END DO 378 END DO 379 380 ! III. Specific grid points 381 ! =========================== 382 ! 383 IF( cp_cfg == "orca" .AND. jp_cfg == 4 ) THEN ! ORCA_R4 configuration: horizontal diffusion in specific area 384 ! ! Gibraltar Strait 385 ij0 = 50 ; ij1 = 53 386 ii0 = 69 ; ii1 = 71 ; uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 387 ij0 = 51 ; ij1 = 53 388 ii0 = 68 ; ii1 = 71 ; vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 389 ii0 = 69 ; ii1 = 71 ; wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 390 ii0 = 69 ; ii1 = 71 ; wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 391 ! 392 ! ! Mediterrannean Sea 393 ij0 = 49 ; ij1 = 56 394 ii0 = 71 ; ii1 = 90 ; uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 395 ij0 = 50 ; ij1 = 56 396 ii0 = 70 ; ii1 = 90 ; vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 397 ii0 = 71 ; ii1 = 90 ; wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 398 ii0 = 71 ; ii1 = 90 ; wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 399 ENDIF 400 401 402 ! IV. Lateral boundary conditions 403 ! =============================== 404 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) 405 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 406 407 408 IF(ln_ctl) THEN 409 CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk) 410 CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk) 411 ENDIF 412 ! 413 414 ELSEIF ( lk_vvl ) THEN 415 416 IF(lwp) THEN 417 WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 418 ENDIF 419 420 ! geopotential diffusion in s-coordinates on tracers and/or momentum 421 ! The slopes of s-surfaces are computed at each time step due to vvl 422 ! The slopes for momentum diffusion are i- or j- averaged of those on tracers 423 424 ! set the slope of diffusion to the slope of s-surfaces 425 ! ( c a u t i o n : minus sign as fsdep has positive value ) 426 DO jj = 2, jpjm1 427 DO ji = fs_2, fs_jpim1 ! vector opt. 428 uslp (ji,jj,1) = -r1_e1u(ji,jj) * ( fsdept_b(ji+1,jj,1) - fsdept_b(ji ,jj ,1) ) * umask(ji,jj,1) 429 vslp (ji,jj,1) = -r1_e2v(ji,jj) * ( fsdept_b(ji,jj+1,1) - fsdept_b(ji ,jj ,1) ) * vmask(ji,jj,1) 430 wslpi(ji,jj,1) = -r1_e1t(ji,jj) * ( fsdepw_b(ji+1,jj,1) - fsdepw_b(ji-1,jj,1) ) * tmask(ji,jj,1) * 0.5 431 wslpj(ji,jj,1) = -r1_e2t(ji,jj) * ( fsdepw_b(ji,jj+1,1) - fsdepw_b(ji,jj-1,1) ) * tmask(ji,jj,1) * 0.5 432 END DO 433 END DO 434 435 DO jk = 2, jpk 436 DO jj = 2, jpjm1 437 DO ji = fs_2, fs_jpim1 ! vector opt. 438 uslp (ji,jj,jk) = -r1_e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk) 439 vslp (ji,jj,jk) = -r1_e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk) 440 wslpi(ji,jj,jk) = -r1_e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * wmask(ji,jj,jk) * 0.5 441 wslpj(ji,jj,jk) = -r1_e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * 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 291 END DO 292 END DO 293 END DO 294 CALL lbc_lnk( zwz, 'T', -1. ) ; CALL lbc_lnk( zww, 'T', -1. ) ! lateral boundary conditions 295 ! 296 ! !* horizontal Shapiro filter 297 DO jk = 2, jpkm1 298 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 299 DO ji = 2, jpim1 300 zcofw = wmask(ji,jj,jk) * z1_16 301 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 302 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 303 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 304 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 305 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 306 307 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 308 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 309 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 310 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 311 & + 4.* zww(ji ,jj ,jk) ) * zcofw 312 END DO 313 END DO 314 DO jj = 3, jpj-2 ! other rows 315 DO ji = fs_2, fs_jpim1 ! vector opt. 316 zcofw = wmask(ji,jj,jk) * z1_16 317 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 318 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 319 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 320 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 321 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 322 323 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 324 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 325 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 326 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 327 & + 4.* zww(ji ,jj ,jk) ) * zcofw 328 END DO 329 END DO 330 ! !* decrease in vicinity of topography 331 DO jj = 2, jpjm1 332 DO ji = fs_2, fs_jpim1 ! vector opt. 333 zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) & 334 & * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 335 wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck 336 wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck 337 END DO 338 END DO 339 END DO 340 341 ! IV. Lateral boundary conditions 342 ! =============================== 343 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) 344 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 345 346 347 IF(ln_ctl) THEN 348 CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk) 349 CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk) 460 350 ENDIF 461 351 ! 462 352 CALL wrk_dealloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv ) 463 CALL wrk_dealloc( jpi,jpj, zhmlpu, zhmlpv)464 353 ! 465 354 IF( nn_timing == 1 ) CALL timing_stop('ldf_slp') … … 468 357 469 358 470 SUBROUTINE ldf_slp_ grif( kt )471 !!---------------------------------------------------------------------- 472 !! *** ROUTINE ldf_slp_ grif***359 SUBROUTINE ldf_slp_triad ( kt ) 360 !!---------------------------------------------------------------------- 361 !! *** ROUTINE ldf_slp_triad *** 473 362 !! 474 363 !! ** Purpose : Compute the squared slopes of neutral surfaces (slope 475 !! of iso-pycnal surfaces referenced locally) (ln_traldf_ grif=T)364 !! of iso-pycnal surfaces referenced locally) (ln_traldf_triad=T) 476 365 !! at W-points using the Griffies quarter-cells. 477 366 !! … … 488 377 REAL(wp) :: zfacti, zfactj ! local scalars 489 378 REAL(wp) :: znot_thru_surface ! local scalars 490 REAL(wp) :: zdit, zdis, zdjt, zdjs, zdkt, zdks, zbu, zbv, zbti, zbtj 379 REAL(wp) :: zdit, zdis, zdkt, zbu, zbti, zisw 380 REAL(wp) :: zdjt, zdjs, zdks, zbv, zbtj, zjsw 491 381 REAL(wp) :: zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_g_raw, zti_g_lim 492 382 REAL(wp) :: zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_g_raw, ztj_g_lim 493 383 REAL(wp) :: zdzrho_raw 384 REAL(wp) :: zbeta0, ze3_e1, ze3_e2 494 385 REAL(wp), POINTER, DIMENSION(:,:) :: z1_mlbw 386 REAL(wp), POINTER, DIMENSION(:,:,:) :: zalbet 495 387 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zdxrho , zdyrho, zdzrho ! Horizontal and vertical density gradients 496 388 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zti_mlb, ztj_mlb ! for Griffies operator only 497 389 !!---------------------------------------------------------------------- 498 390 ! 499 IF( nn_timing == 1 ) CALL timing_start('ldf_slp_ grif')391 IF( nn_timing == 1 ) CALL timing_start('ldf_slp_triad') 500 392 ! 501 393 CALL wrk_alloc( jpi,jpj, z1_mlbw ) 394 CALL wrk_alloc( jpi,jpj,jpk, zalbet ) 502 395 CALL wrk_alloc( jpi,jpj,jpk,2, zdxrho , zdyrho, zdzrho, klstart = 0 ) 503 396 CALL wrk_alloc( jpi,jpj, 2,2, zti_mlb, ztj_mlb, kkstart = 0, klstart = 0 ) … … 519 412 zdxrho_raw = ( - rab_b(ji+ip,jj ,jk,jp_tem) * zdit + rab_b(ji+ip,jj ,jk,jp_sal) * zdis ) * r1_e1u(ji,jj) 520 413 zdyrho_raw = ( - rab_b(ji ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji ,jj+jp,jk,jp_sal) * zdjs ) * r1_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)414 zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 415 zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 523 416 END DO 524 417 END DO … … 552 445 zdks = 0._wp 553 446 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 >= repsln447 zdzrho_raw = ( - zalbet(ji,jj,jk) * zdkt + zbeta0*zdks ) / fse3w(ji,jj,jk+kp) 448 zdzrho(ji,jj,jk,kp) = - MIN( - repsln , zdzrho_raw ) ! force zdzrho >= repsln 556 449 END DO 557 450 END DO … … 586 479 ! 587 480 jk = nmln(ji+ip,jj) + 1 588 IF( jk .GT. 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) ) * r1_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 ) 481 IF( jk > mbkt(ji+ip,jj) ) THEN ! ML reaches bottom 482 zti_mlb(ji+ip,jj ,1-ip,kp) = 0.0_wp 483 ELSE 484 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 485 zti_g_raw = ( zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp) & 486 & - ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) * r1_e1u(ji,jj) ) * umask(ji,jj,jk) 487 ze3_e1 = fse3w(ji+ip,jj,jk-kp) * r1_e1u(ji,jj) 488 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 489 ENDIF 596 490 ! 597 491 jk = nmln(ji,jj+jp) + 1 598 492 IF( jk .GT. mbkt(ji,jj+jp) ) THEN !ML reaches bottom 599 ztj_mlb(ji ,jj+jp,1-jp,kp) = 0.0_wp493 ztj_mlb(ji ,jj+jp,1-jp,kp) = 0.0_wp 600 494 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) ) * r1_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 ) 495 ztj_g_raw = ( zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp) & 496 & - ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) 497 ze3_e2 = fse3w(ji,jj+jp,jk-kp) / e2v(ji,jj) 498 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 499 ENDIF 605 500 END DO … … 628 523 ! raw slopes: unmasked unbounded slopes (relative to geopotential (zti_g) and model surface (zti) 629 524 ! 630 zti_raw = zdxrho(ji+ip,jj ,jk,1-ip) / zdzrho(ji+ip,jj ,jk,kp) 525 zti_raw = zdxrho(ji+ip,jj ,jk,1-ip) / zdzrho(ji+ip,jj ,jk,kp) ! unmasked 631 526 ztj_raw = zdyrho(ji ,jj+jp,jk,1-jp) / zdzrho(ji ,jj+jp,jk,kp) 632 527 ! 633 528 ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface 634 529 zti_coord = znot_thru_surface * ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) * r1_e1u(ji,jj) … … 636 531 zti_g_raw = zti_raw - zti_coord ! ref to geopot surfaces 637 532 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 ) 533 ! additional limit required in bilaplacian case 534 ze3_e1 = fse3w(ji+ip,jj ,jk+kp) * r1_e1u(ji,jj) 535 ze3_e2 = fse3w(ji ,jj+jp,jk+kp) * r1_e2v(ji,jj) 536 ! NB: hard coded factor 5 (can be a namelist parameter...) 537 zti_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1, ABS( zti_g_raw ) ), zti_g_raw ) 538 ztj_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2, ABS( ztj_g_raw ) ), ztj_g_raw ) 640 539 ! 641 540 ! Below ML use limited zti_g as is & mask … … 666 565 ! 667 566 IF( ln_triad_iso ) THEN 668 zti_raw = zti_lim* *2/ zti_raw669 ztj_raw = ztj_lim* *2/ ztj_raw567 zti_raw = zti_lim*zti_lim / zti_raw 568 ztj_raw = ztj_lim*ztj_lim / ztj_raw 670 569 zti_raw = SIGN( MIN( ABS(zti_lim), ABS( zti_raw ) ), zti_raw ) 671 570 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 571 zti_lim = zfacti * zti_lim + ( 1._wp - zfacti ) * zti_raw 572 ztj_lim = zfactj * ztj_lim + ( 1._wp - zfactj ) * ztj_raw 676 573 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 = e1e2u(ji ,jj) * fse3u(ji ,jj,jk ) 681 zbv = e1e2v(ji ,jj) * fse3v(ji ,jj,jk ) 682 zbti = e1e2t(ji+ip,jj) * fse3w(ji+ip,jj,jk+kp) 683 zbtj = e1e2t(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 574 ! ! switching triad scheme 575 zisw = (rn_sw_triad - 1._wp ) + rn_sw_triad & 576 & * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - ip ) * SIGN( 1._wp , zdxrho(ji+ip,jj,jk,1-ip) ) ) 577 zjsw = (rn_sw_triad - 1._wp ) + rn_sw_triad & 578 & * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - jp ) * SIGN( 1._wp , zdyrho(ji,jj+jp,jk,1-jp) ) ) 579 ! 580 triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim * zisw 581 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim * zjsw 582 ! 583 zbu = e1e2u(ji ,jj ) * fse3u(ji ,jj ,jk ) 584 zbv = e1e2v(ji ,jj ) * fse3v(ji ,jj ,jk ) 585 zbti = e1e2t(ji+ip,jj ) * fse3w(ji+ip,jj ,jk+kp) 586 zbtj = e1e2t(ji ,jj+jp) * fse3w(ji ,jj+jp,jk+kp) 587 ! 588 wslp2(ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim*zti_g_lim ! masked 589 wslp2(ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_g_lim*ztj_g_lim 688 590 END DO 689 591 END DO … … 697 599 ! 698 600 CALL wrk_dealloc( jpi,jpj, z1_mlbw ) 601 CALL wrk_dealloc( jpi,jpj,jpk, zalbet ) 699 602 CALL wrk_dealloc( jpi,jpj,jpk,2, zdxrho , zdyrho, zdzrho, klstart = 0 ) 700 603 CALL wrk_dealloc( jpi,jpj, 2,2, zti_mlb, ztj_mlb, kkstart = 0, klstart = 0 ) 701 604 ! 702 IF( nn_timing == 1 ) CALL timing_stop('ldf_slp_ grif')703 ! 704 END SUBROUTINE ldf_slp_ grif605 IF( nn_timing == 1 ) CALL timing_stop('ldf_slp_triad') 606 ! 607 END SUBROUTINE ldf_slp_triad 705 608 706 609 … … 728 631 INTEGER :: ji , jj , jk ! dummy loop indices 729 632 INTEGER :: iku, ikv, ik, ikm1 ! local integers 730 REAL(wp) :: zeps, zm1_g, zm1_2g 633 REAL(wp) :: zeps, zm1_g, zm1_2g, z1_slpmax ! local scalars 731 634 REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - - 732 635 REAL(wp) :: zcj, zfj, zav, zbv, zaj, zbj ! - - … … 739 642 zm1_g = -1.0_wp / grav 740 643 zm1_2g = -0.5_wp / grav 644 z1_slpmax = 1._wp / rn_slpmax 741 645 ! 742 646 uslpml (1,:) = 0._wp ; uslpml (jpi,:) = 0._wp … … 750 654 DO ji = 1, jpi 751 655 ik = nmln(ji,jj) - 1 752 IF( jk <= ik .AND. jk >= mikt(ji,jj) ) THEN 753 omlmask(ji,jj,jk) = 1._wp 754 ELSE 755 omlmask(ji,jj,jk) = 0._wp 656 IF( jk <= ik ) THEN ; omlmask(ji,jj,jk) = 1._wp 657 ELSE ; omlmask(ji,jj,jk) = 0._wp 756 658 ENDIF 757 659 END DO … … 775 677 ! 776 678 ! !- vertical density gradient for u- and v-slopes (from dzr at T-point) 777 iku = MIN( MAX( miku(ji,jj)+1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1 ) ! ML (MAX of T-pts, bound by jpkm1)778 ikv = MIN( MAX( mikv(ji,jj)+1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) !679 iku = MIN( MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1 ) ! ML (MAX of T-pts, bound by jpkm1) 680 ikv = MIN( MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) ! 779 681 zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj ,iku) ) 780 682 zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji ,jj+1,ikv) ) … … 784 686 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 785 687 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 786 zbu = MIN( zbu , - 100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau ) )787 zbv = MIN( zbv , - 100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav ) )688 zbu = MIN( zbu , - z1_slpmax * ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau ) ) 689 zbv = MIN( zbv , - z1_slpmax * ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav ) ) 788 690 ! !- Slope at u- & v-points (uslpml, vslpml) 789 691 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku) … … 810 712 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zaj ) ) 811 713 ! !- i- & j-slope at w-points (wslpiml, wslpjml) 812 wslpiml(ji,jj) = zai / ( zbi - zeps ) * wmask (ji,jj,ik)813 wslpjml(ji,jj) = zaj / ( zbj - zeps ) * wmask (ji,jj,ik)714 wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) 715 wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik) 814 716 END DO 815 717 END DO … … 829 731 !! ** Purpose : Initialization for the isopycnal slopes computation 830 732 !! 831 !! ** Method : read the nammbf namelist and check the parameter 832 !! values called by tra_dmp at the first timestep (nit000) 733 !! ** Method : 833 734 !!---------------------------------------------------------------------- 834 735 INTEGER :: ji, jj, jk ! dummy loop indices … … 843 744 WRITE(numout,*) '~~~~~~~~~~~~' 844 745 ENDIF 845 846 IF( ln_traldf_grif ) THEN ! Griffies operator : triad of slopes 847 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 ) 848 ALLOCATE( triadi (jpi,jpj,jpk,0:1,0:1) , triadj (jpi,jpj,jpk,0:1,0:1) , STAT=ierr ) 849 IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Griffies operator slope' ) 850 ! 746 ! 747 ALLOCATE( ah_wslp2(jpi,jpj,jpk) , akz(jpi,jpj,jpk) , STAT=ierr ) 748 IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate ah_slp2 or akz' ) 749 ! 750 IF( ln_traldf_triad ) THEN ! Griffies operator : triad of slopes 751 IF(lwp) WRITE(numout,*) ' Griffies (triad) operator initialisation' 752 ALLOCATE( triadi_g(jpi,jpj,jpk,0:1,0:1) , triadj_g(jpi,jpj,jpk,0:1,0:1) , & 753 & triadi (jpi,jpj,jpk,0:1,0:1) , triadj (jpi,jpj,jpk,0:1,0:1) , & 754 & wslp2 (jpi,jpj,jpk) , STAT=ierr ) 755 IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Griffies operator slope' ) 851 756 IF( ln_dynldf_iso ) CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' ) 852 757 ! 853 758 ELSE ! Madec operator : slopes at u-, v-, and w-points 854 ALLOCATE( uslp(jpi,jpj,jpk) , vslp(jpi,jpj,jpk) , wslpi(jpi,jpj,jpk) , wslpj(jpi,jpj,jpk) , & 855 & omlmask(jpi,jpj,jpk) , uslpml(jpi,jpj) , vslpml(jpi,jpj) , wslpiml(jpi,jpj) , wslpjml(jpi,jpj) , STAT=ierr ) 759 IF(lwp) WRITE(numout,*) ' Madec operator initialisation' 760 ALLOCATE( omlmask(jpi,jpj,jpk) , & 761 & uslp(jpi,jpj,jpk) , uslpml(jpi,jpj) , wslpi(jpi,jpj,jpk) , wslpiml(jpi,jpj) , & 762 & vslp(jpi,jpj,jpk) , vslpml(jpi,jpj) , wslpj(jpi,jpj,jpk) , wslpjml(jpi,jpj) , STAT=ierr ) 856 763 IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Madec operator slope ' ) 857 764 … … 863 770 wslpj(:,:,:) = 0._wp ; wslpjml(:,:) = 0._wp 864 771 865 IF(ln_sco .AND. (ln_traldf_hor .OR. ln_dynldf_hor )) THEN 866 IF(lwp) WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 867 868 ! geopotential diffusion in s-coordinates on tracers and/or momentum 869 ! The slopes of s-surfaces are computed once (no call to ldfslp in step) 870 ! The slopes for momentum diffusion are i- or j- averaged of those on tracers 871 872 ! set the slope of diffusion to the slope of s-surfaces 873 ! ( c a u t i o n : minus sign as fsdep has positive value ) 874 DO jk = 1, jpk 875 DO jj = 2, jpjm1 876 DO ji = fs_2, fs_jpim1 ! vector opt. 877 uslp (ji,jj,jk) = -r1_e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk) 878 vslp (ji,jj,jk) = -r1_e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk) 879 wslpi(ji,jj,jk) = -r1_e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * wmask(ji,jj,jk) * 0.5 880 wslpj(ji,jj,jk) = -r1_e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * wmask(ji,jj,jk) * 0.5 881 END DO 882 END DO 883 END DO 884 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) ! Lateral boundary conditions 885 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 886 ENDIF 772 !!gm I no longer understand this..... 773 !!gm IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN 774 ! IF(lwp) WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 775 ! 776 ! ! geopotential diffusion in s-coordinates on tracers and/or momentum 777 ! ! The slopes of s-surfaces are computed once (no call to ldfslp in step) 778 ! ! The slopes for momentum diffusion are i- or j- averaged of those on tracers 779 ! 780 ! ! set the slope of diffusion to the slope of s-surfaces 781 ! ! ( c a u t i o n : minus sign as fsdep has positive value ) 782 ! DO jk = 1, jpk 783 ! DO jj = 2, jpjm1 784 ! DO ji = fs_2, fs_jpim1 ! vector opt. 785 ! uslp (ji,jj,jk) = - ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 786 ! vslp (ji,jj,jk) = - ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 787 ! wslpi(ji,jj,jk) = - ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * r1_e1t(ji,jj) * wmask(ji,jj,jk) * 0.5 788 ! wslpj(ji,jj,jk) = - ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * r1_e2t(ji,jj) * wmask(ji,jj,jk) * 0.5 789 ! END DO 790 ! END DO 791 ! END DO 792 ! CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) ! Lateral boundary conditions 793 ! CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 794 !!gm ENDIF 887 795 ENDIF 888 796 ! … … 890 798 ! 891 799 END SUBROUTINE ldf_slp_init 892 893 #else894 !!------------------------------------------------------------------------895 !! Dummy module : NO Rotation of lateral mixing tensor896 !!------------------------------------------------------------------------897 LOGICAL, PUBLIC, PARAMETER :: lk_ldfslp = .FALSE. !: slopes flag898 CONTAINS899 SUBROUTINE ldf_slp( kt, prd, pn2 ) ! Dummy routine900 INTEGER, INTENT(in) :: kt901 REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2902 WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1)903 END SUBROUTINE ldf_slp904 SUBROUTINE ldf_slp_grif( kt ) ! Dummy routine905 INTEGER, INTENT(in) :: kt906 WRITE(*,*) 'ldf_slp_grif: You should not have seen this print! error?', kt907 END SUBROUTINE ldf_slp_grif908 SUBROUTINE ldf_slp_init ! Dummy routine909 END SUBROUTINE ldf_slp_init910 #endif911 800 912 801 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
