Changeset 5951 for branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
- Timestamp:
- 2015-11-30T12:48:01+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r5950 r5951 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, iku ! temporary integer 106 115 INTEGER :: ij0, ij1, ikv ! 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 ! - -112 120 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz, zww 113 121 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdzr 114 122 REAL(wp), POINTER, DIMENSION(:,:,:) :: zgru, zgrv 115 REAL(wp), POINTER, DIMENSION(:,: ) :: zhmlpu, zhmlpv116 123 !!---------------------------------------------------------------------- 117 124 ! … … 119 126 ! 120 127 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) / e1u(ji,jj) 211 zav = zgrv(ji,jj,jk) / 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 128 129 zeps = 1.e-20_wp !== Local constant initialization ==! 130 z1_16 = 1.0_wp / 16._wp 131 zm1_g = -1.0_wp / grav 132 zm1_2g = -0.5_wp / grav 133 z1_slpmax = 1._wp / rn_slpmax 134 ! 135 zww(:,:,:) = 0._wp 136 zwz(:,:,:) = 0._wp 137 ! 138 DO jk = 1, jpk !== i- & j-gradient of density ==! 139 DO jj = 1, jpjm1 140 DO ji = 1, fs_jpim1 ! vector opt. 141 zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj ,jk) - prd(ji,jj,jk) ) 142 zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji ,jj+1,jk) - prd(ji,jj,jk) ) 143 END DO 144 END DO 145 END DO 146 IF( ln_zps ) THEN ! partial steps correction at the bottom ocean level 147 DO jj = 1, jpjm1 148 DO ji = 1, jpim1 149 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 150 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 151 END DO 152 END DO 153 ENDIF 154 ! 155 zdzr(:,:,1) = 0._wp !== Local vertical density gradient at T-point == ! (evaluated from N^2) 156 DO jk = 2, jpkm1 157 ! ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point 158 ! ! trick: tmask(ik ) = 0 => all pn2 = 0 => zdzr = 0 159 ! ! else tmask(ik+1) = 0 => pn2(ik+1) = 0 => zdzr divides by 1 160 ! ! umask(ik+1) /= 0 => all pn2 /= 0 => zdzr divides by 2 161 ! ! NB: 1/(tmask+1) = (1-.5*tmask) substitute a / by a * ==> faster 162 zdzr(:,:,jk) = zm1_g * ( prd(:,:,jk) + 1._wp ) & 163 & * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask(:,:,jk+1) ) 164 END DO 165 ! 166 ! !== Slopes just below the mixed layer ==! 167 CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr ) ! output: uslpml, vslpml, wslpiml, wslpjml 168 169 170 ! I. slopes at u and v point | uslp = d/di( prd ) / d/dz( prd ) 171 ! =========================== | vslp = d/dj( prd ) / d/dz( prd ) 172 ! 173 DO jk = 2, jpkm1 !* Slopes at u and v points 174 DO jj = 2, jpjm1 175 DO ji = fs_2, fs_jpim1 ! vector opt. 176 ! ! horizontal and vertical density gradient at u- and v-points 177 zau = zgru(ji,jj,jk) * r1_e1u(ji,jj) 178 zav = zgrv(ji,jj,jk) * r1_e2v(ji,jj) 179 zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj ,jk) ) 180 zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji ,jj+1,jk) ) 181 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 182 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 183 zbu = MIN( zbu, - z1_slpmax * ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,jk)* ABS( zau ) ) 184 zbv = MIN( zbv, - z1_slpmax * ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav ) ) 185 ! ! uslp and vslp output in zwz and zww, resp. 186 zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 187 zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 188 zwz(ji,jj,jk) = ( ( 1. - zfi) * zau / ( zbu - zeps ) & 189 & + zfi * uslpml(ji,jj) & 190 & * 0.5_wp * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) ) & 191 & / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5._wp ) ) * umask(ji,jj,jk) 192 zww(ji,jj,jk) = ( ( 1. - zfj) * zav / ( zbv - zeps ) & 193 & + zfj * vslpml(ji,jj) & 194 & * 0.5_wp * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) ) & 195 & / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk) 235 196 !!gm modif to suppress omlmask.... (as in Griffies case) 236 ! 237 ! 238 ! 239 ! 240 ! 241 ! 242 ! 197 ! ! ! jk must be >= ML level for zf=1. otherwise zf=0. 198 ! zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp ) 199 ! zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp ) 200 ! zci = 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) ) 201 ! zcj = 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) ) 202 ! zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk) 203 ! zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk) 243 204 !!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) 205 END DO 206 END DO 207 END DO 208 CALL lbc_lnk( zwz, 'U', -1. ) ; CALL lbc_lnk( zww, 'V', -1. ) ! lateral boundary conditions 209 ! 210 ! !* horizontal Shapiro filter 211 DO jk = 2, jpkm1 212 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 213 DO ji = 2, jpim1 214 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 215 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 216 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 217 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 218 & + 4.* zwz(ji ,jj ,jk) ) 219 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 220 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 221 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 222 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 223 & + 4.* zww(ji,jj ,jk) ) 224 END DO 225 END DO 226 DO jj = 3, jpj-2 ! other rows 227 DO ji = fs_2, fs_jpim1 ! vector opt. 228 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 229 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 230 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 231 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 232 & + 4.* zwz(ji ,jj ,jk) ) 233 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 234 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 235 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 236 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 237 & + 4.* zww(ji,jj ,jk) ) 238 END DO 239 END DO 240 ! !* decrease along coastal boundaries 241 DO jj = 2, jpjm1 242 DO ji = fs_2, fs_jpim1 ! vector opt. 243 uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk ) ) * 0.5_wp & 244 & * ( umask(ji,jj ,jk) + umask(ji,jj ,jk+1) ) * 0.5_wp 245 vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk ) ) * 0.5_wp & 246 & * ( vmask(ji ,jj,jk) + vmask(ji ,jj,jk+1) ) * 0.5_wp 247 END DO 248 END DO 249 END DO 250 251 252 ! II. slopes at w point | wslpi = mij( d/di( prd ) / d/dz( prd ) 253 ! =========================== | wslpj = mij( d/dj( prd ) / d/dz( prd ) 254 ! 255 DO jk = 2, jpkm1 256 DO jj = 2, jpjm1 257 DO ji = fs_2, fs_jpim1 ! vector opt. 258 ! !* Local vertical density gradient evaluated from N^2 259 zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) 260 ! !* Slopes at w point 261 ! ! i- & j-gradient of density at w-points 262 zci = MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk ) & 263 & + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps ) * e1t(ji,jj) 264 zcj = MAX( vmask(ji,jj-1,jk ) + vmask(ji,jj,jk-1) & 265 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk ) , zeps ) * e2t(ji,jj) 266 zai = ( zgru (ji-1,jj,jk ) + zgru (ji,jj,jk-1) & 267 & + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk ) ) / zci * tmask (ji,jj,jk) 268 zaj = ( zgrv (ji,jj-1,jk ) + zgrv (ji,jj,jk-1) & 269 & + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk ) ) / zcj * tmask (ji,jj,jk) 270 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 271 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 272 zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zai ) ) 273 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zaj ) ) 274 ! ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 275 zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) ) ! zfk=1 in the ML otherwise zfk=0 276 zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10._wp ) 277 zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk ) * tmask(ji,jj,jk) 278 zww(ji,jj,jk) = ( zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk ) * tmask(ji,jj,jk) 322 279 323 280 !!gm modif to suppress omlmask.... (as in Griffies operator) 324 ! 325 ! 326 ! 327 ! 328 ! 281 ! ! ! jk must be >= ML level for zfk=1. otherwise zfk=0. 282 ! zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp ) 283 ! zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10. ) 284 ! zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk) 285 ! zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk) 329 286 !!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) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,1) - fsdept_b(ji ,jj ,1) ) * umask(ji,jj,1) 429 vslp(ji,jj,1) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,1) - fsdept_b(ji ,jj ,1) ) * vmask(ji,jj,1) 430 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 431 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 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) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk) 439 vslp(ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk) 440 wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) & 441 & * wmask(ji,jj,jk) * 0.5 442 wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) & 443 & * wmask(ji,jj,jk) * 0.5 444 END DO 445 END DO 446 END DO 447 448 ! Lateral boundary conditions on the slopes 449 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) 450 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 451 452 if( kt == nit000 ) then 453 IF(lwp) WRITE(numout,*) ' max slop: u',SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)), & 454 & ' wi', sqrt(MAXVAL(wslpi)), ' wj', sqrt(MAXVAL(wslpj)) 455 endif 456 457 IF(ln_ctl) THEN 458 CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk) 459 CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk) 460 ENDIF 461 287 END DO 288 END DO 289 END DO 290 CALL lbc_lnk( zwz, 'T', -1. ) ; CALL lbc_lnk( zww, 'T', -1. ) ! lateral boundary conditions 291 ! 292 ! !* horizontal Shapiro filter 293 DO jk = 2, jpkm1 294 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 295 DO ji = 2, jpim1 296 zcofw = wmask(ji,jj,jk) * z1_16 297 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 298 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 299 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 300 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 301 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 302 303 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 304 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 305 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 306 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 307 & + 4.* zww(ji ,jj ,jk) ) * zcofw 308 END DO 309 END DO 310 DO jj = 3, jpj-2 ! other rows 311 DO ji = fs_2, fs_jpim1 ! vector opt. 312 zcofw = wmask(ji,jj,jk) * z1_16 313 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 314 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 315 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 316 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 317 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 318 319 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 320 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 321 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 322 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 323 & + 4.* zww(ji ,jj ,jk) ) * zcofw 324 END DO 325 END DO 326 ! !* decrease in vicinity of topography 327 DO jj = 2, jpjm1 328 DO ji = fs_2, fs_jpim1 ! vector opt. 329 zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) & 330 & * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 331 wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck 332 wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck 333 END DO 334 END DO 335 END DO 336 337 ! IV. Lateral boundary conditions 338 ! =============================== 339 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) 340 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 341 342 343 IF(ln_ctl) THEN 344 CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk) 345 CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk) 462 346 ENDIF 463 347 ! 464 348 CALL wrk_dealloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv ) 465 CALL wrk_dealloc( jpi,jpj, zhmlpu, zhmlpv)466 349 ! 467 350 IF( nn_timing == 1 ) CALL timing_stop('ldf_slp') … … 470 353 471 354 472 SUBROUTINE ldf_slp_ grif( kt )473 !!---------------------------------------------------------------------- 474 !! *** ROUTINE ldf_slp_ grif***355 SUBROUTINE ldf_slp_triad ( kt ) 356 !!---------------------------------------------------------------------- 357 !! *** ROUTINE ldf_slp_triad *** 475 358 !! 476 359 !! ** Purpose : Compute the squared slopes of neutral surfaces (slope 477 !! of iso-pycnal surfaces referenced locally) (ln_traldf_ grif=T)360 !! of iso-pycnal surfaces referenced locally) (ln_traldf_triad=T) 478 361 !! at W-points using the Griffies quarter-cells. 479 362 !! … … 490 373 REAL(wp) :: zfacti, zfactj ! local scalars 491 374 REAL(wp) :: znot_thru_surface ! local scalars 492 REAL(wp) :: zdit, zdis, zdjt, zdjs, zdkt, zdks, zbu, zbv, zbti, zbtj 375 REAL(wp) :: zdit, zdis, zdkt, zbu, zbti, zisw 376 REAL(wp) :: zdjt, zdjs, zdks, zbv, zbtj, zjsw 493 377 REAL(wp) :: zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_g_raw, zti_g_lim 494 378 REAL(wp) :: zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_g_raw, ztj_g_lim 495 379 REAL(wp) :: zdzrho_raw 380 REAL(wp) :: zbeta0, ze3_e1, ze3_e2 496 381 REAL(wp), POINTER, DIMENSION(:,:) :: z1_mlbw 382 REAL(wp), POINTER, DIMENSION(:,:,:) :: zalbet 497 383 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zdxrho , zdyrho, zdzrho ! Horizontal and vertical density gradients 498 384 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zti_mlb, ztj_mlb ! for Griffies operator only 499 385 !!---------------------------------------------------------------------- 500 386 ! 501 IF( nn_timing == 1 ) CALL timing_start('ldf_slp_ grif')387 IF( nn_timing == 1 ) CALL timing_start('ldf_slp_triad') 502 388 ! 503 389 CALL wrk_alloc( jpi,jpj, z1_mlbw ) 390 CALL wrk_alloc( jpi,jpj,jpk, zalbet ) 504 391 CALL wrk_alloc( jpi,jpj,jpk,2, zdxrho , zdyrho, zdzrho, klstart = 0 ) 505 392 CALL wrk_alloc( jpi,jpj, 2,2, zti_mlb, ztj_mlb, kkstart = 0, klstart = 0 ) … … 519 406 zdjt = ( tsb(ji,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! j-gradient of T & S at v-point 520 407 zdjs = ( tsb(ji,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 521 zdxrho_raw = ( - rab_b(ji+ip,jj ,jk,jp_tem) * zdit + rab_b(ji+ip,jj ,jk,jp_sal) * zdis ) /e1u(ji,jj)522 zdyrho_raw = ( - rab_b(ji ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji ,jj+jp,jk,jp_sal) * zdjs ) /e2v(ji,jj)523 zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw) ! keep the sign524 zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw)408 zdxrho_raw = ( - rab_b(ji+ip,jj ,jk,jp_tem) * zdit + rab_b(ji+ip,jj ,jk,jp_sal) * zdis ) * r1_e1u(ji,jj) 409 zdyrho_raw = ( - rab_b(ji ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji ,jj+jp,jk,jp_sal) * zdjs ) * r1_e2v(ji,jj) 410 zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 411 zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 525 412 END DO 526 413 END DO … … 533 420 zdit = gtsu(ji,jj,jp_tem) ; zdjt = gtsv(ji,jj,jp_tem) ! i- & j-gradient of Temperature 534 421 zdis = gtsu(ji,jj,jp_sal) ; zdjs = gtsv(ji,jj,jp_sal) ! i- & j-gradient of Salinity 535 zdxrho_raw = ( - rab_b(ji+ip,jj ,iku,jp_tem) * zdit + rab_b(ji+ip,jj ,iku,jp_sal) * zdis ) /e1u(ji,jj)536 zdyrho_raw = ( - rab_b(ji ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji ,jj+jp,ikv,jp_sal) * zdjs ) /e2v(ji,jj)422 zdxrho_raw = ( - rab_b(ji+ip,jj ,iku,jp_tem) * zdit + rab_b(ji+ip,jj ,iku,jp_sal) * zdis ) * r1_e1u(ji,jj) 423 zdyrho_raw = ( - rab_b(ji ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji ,jj+jp,ikv,jp_sal) * zdjs ) * r1_e2v(ji,jj) 537 424 zdxrho(ji+ip,jj ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 538 425 zdyrho(ji ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) … … 554 441 zdks = 0._wp 555 442 ENDIF 556 zdzrho_raw = ( - rab_b(ji,jj,jk,jp_tem) * zdkt + rab_b(ji,jj,jk,jp_sal) *zdks ) / fse3w(ji,jj,jk+kp)557 zdzrho(ji,jj,jk,kp) = - MIN( - repsln , zdzrho_raw ) ! force zdzrho >= repsln443 zdzrho_raw = ( - zalbet(ji,jj,jk) * zdkt + zbeta0*zdks ) / fse3w(ji,jj,jk+kp) 444 zdzrho(ji,jj,jk,kp) = - MIN( - repsln , zdzrho_raw ) ! force zdzrho >= repsln 558 445 END DO 559 446 END DO … … 588 475 ! 589 476 jk = nmln(ji+ip,jj) + 1 590 IF( jk .GT. mbkt(ji+ip,jj) ) THEN !ML reaches bottom 591 zti_mlb(ji+ip,jj ,1-ip,kp) = 0.0_wp 592 ELSE 593 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 594 zti_g_raw = ( zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp) & 595 & - ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) / e1u(ji,jj) ) * umask(ji,jj,jk) 596 zti_mlb(ji+ip,jj ,1-ip,kp) = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 477 IF( jk > mbkt(ji+ip,jj) ) THEN ! ML reaches bottom 478 zti_mlb(ji+ip,jj ,1-ip,kp) = 0.0_wp 479 ELSE 480 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 481 zti_g_raw = ( zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp) & 482 & - ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) * r1_e1u(ji,jj) ) * umask(ji,jj,jk) 483 ze3_e1 = fse3w(ji+ip,jj,jk-kp) * r1_e1u(ji,jj) 484 zti_mlb(ji+ip,jj ,1-ip,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1 , ABS( zti_g_raw ) ), zti_g_raw ) 597 485 ENDIF 598 486 ! 599 487 jk = nmln(ji,jj+jp) + 1 600 488 IF( jk .GT. mbkt(ji,jj+jp) ) THEN !ML reaches bottom 601 ztj_mlb(ji ,jj+jp,1-jp,kp) = 0.0_wp489 ztj_mlb(ji ,jj+jp,1-jp,kp) = 0.0_wp 602 490 ELSE 603 ztj_g_raw = ( zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp) & 604 & - ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) 605 ztj_mlb(ji ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) 491 ztj_g_raw = ( zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp) & 492 & - ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) 493 ze3_e2 = fse3w(ji,jj+jp,jk-kp) / e2v(ji,jj) 494 ztj_mlb(ji ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2 , ABS( ztj_g_raw ) ), ztj_g_raw ) 606 495 ENDIF 607 496 END DO … … 632 521 zti_raw = zdxrho(ji+ip,jj ,jk,1-ip) / zdzrho(ji+ip,jj ,jk,kp) ! unmasked 633 522 ztj_raw = zdyrho(ji ,jj+jp,jk,1-jp) / zdzrho(ji ,jj+jp,jk,kp) 634 523 ! 635 524 ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface 636 zti_coord = znot_thru_surface * ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) /e1u(ji,jj)637 ztj_coord = znot_thru_surface * ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj)! unmasked525 zti_coord = znot_thru_surface * ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) * r1_e1u(ji,jj) 526 ztj_coord = znot_thru_surface * ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) * r1_e2v(ji,jj) ! unmasked 638 527 zti_g_raw = zti_raw - zti_coord ! ref to geopot surfaces 639 528 ztj_g_raw = ztj_raw - ztj_coord 640 zti_g_lim = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 641 ztj_g_lim = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) 529 ! additional limit required in bilaplacian case 530 ze3_e1 = fse3w(ji+ip,jj ,jk+kp) * r1_e1u(ji,jj) 531 ze3_e2 = fse3w(ji ,jj+jp,jk+kp) * r1_e2v(ji,jj) 532 ! NB: hard coded factor 5 (can be a namelist parameter...) 533 zti_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1, ABS( zti_g_raw ) ), zti_g_raw ) 534 ztj_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2, ABS( ztj_g_raw ) ), ztj_g_raw ) 642 535 ! 643 536 ! Below ML use limited zti_g as is & mask … … 668 561 ! 669 562 IF( ln_triad_iso ) THEN 670 zti_raw = zti_lim* *2/ zti_raw671 ztj_raw = ztj_lim* *2/ ztj_raw563 zti_raw = zti_lim*zti_lim / zti_raw 564 ztj_raw = ztj_lim*ztj_lim / ztj_raw 672 565 zti_raw = SIGN( MIN( ABS(zti_lim), ABS( zti_raw ) ), zti_raw ) 673 566 ztj_raw = SIGN( MIN( ABS(ztj_lim), ABS( ztj_raw ) ), ztj_raw ) 674 zti_lim = zfacti * zti_lim & 675 & + ( 1._wp - zfacti ) * zti_raw 676 ztj_lim = zfactj * ztj_lim & 677 & + ( 1._wp - zfactj ) * ztj_raw 567 zti_lim = zfacti * zti_lim + ( 1._wp - zfacti ) * zti_raw 568 ztj_lim = zfactj * ztj_lim + ( 1._wp - zfactj ) * ztj_raw 678 569 ENDIF 679 triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim 680 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim 681 ! 682 zbu = e1u(ji ,jj) * e2u(ji ,jj) * fse3u(ji ,jj,jk ) 683 zbv = e1v(ji ,jj) * e2v(ji ,jj) * fse3v(ji ,jj,jk ) 684 zbti = e1t(ji+ip,jj) * e2t(ji+ip,jj) * fse3w(ji+ip,jj,jk+kp) 685 zbtj = e1t(ji,jj+jp) * e2t(ji,jj+jp) * fse3w(ji,jj+jp,jk+kp) 686 ! 687 !!gm this may inhibit vectorization on Vect Computers, and even on scalar computers.... ==> to be checked 688 wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim**2 ! masked 689 wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_g_lim**2 570 ! ! switching triad scheme 571 zisw = (rn_sw_triad - 1._wp ) + rn_sw_triad & 572 & * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - ip ) * SIGN( 1._wp , zdxrho(ji+ip,jj,jk,1-ip) ) ) 573 zjsw = (rn_sw_triad - 1._wp ) + rn_sw_triad & 574 & * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - jp ) * SIGN( 1._wp , zdyrho(ji,jj+jp,jk,1-jp) ) ) 575 ! 576 triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim * zisw 577 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim * zjsw 578 ! 579 zbu = e1e2u(ji ,jj ) * fse3u(ji ,jj ,jk ) 580 zbv = e1e2v(ji ,jj ) * fse3v(ji ,jj ,jk ) 581 zbti = e1e2t(ji+ip,jj ) * fse3w(ji+ip,jj ,jk+kp) 582 zbtj = e1e2t(ji ,jj+jp) * fse3w(ji ,jj+jp,jk+kp) 583 ! 584 wslp2(ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim*zti_g_lim ! masked 585 wslp2(ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_g_lim*ztj_g_lim 690 586 END DO 691 587 END DO … … 699 595 ! 700 596 CALL wrk_dealloc( jpi,jpj, z1_mlbw ) 597 CALL wrk_dealloc( jpi,jpj,jpk, zalbet ) 701 598 CALL wrk_dealloc( jpi,jpj,jpk,2, zdxrho , zdyrho, zdzrho, klstart = 0 ) 702 599 CALL wrk_dealloc( jpi,jpj, 2,2, zti_mlb, ztj_mlb, kkstart = 0, klstart = 0 ) 703 600 ! 704 IF( nn_timing == 1 ) CALL timing_stop('ldf_slp_ grif')705 ! 706 END SUBROUTINE ldf_slp_ grif601 IF( nn_timing == 1 ) CALL timing_stop('ldf_slp_triad') 602 ! 603 END SUBROUTINE ldf_slp_triad 707 604 708 605 … … 730 627 INTEGER :: ji , jj , jk ! dummy loop indices 731 628 INTEGER :: iku, ikv, ik, ikm1 ! local integers 732 REAL(wp) :: zeps, zm1_g, zm1_2g 629 REAL(wp) :: zeps, zm1_g, zm1_2g, z1_slpmax ! local scalars 733 630 REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - - 734 631 REAL(wp) :: zcj, zfj, zav, zbv, zaj, zbj ! - - … … 741 638 zm1_g = -1.0_wp / grav 742 639 zm1_2g = -0.5_wp / grav 640 z1_slpmax = 1._wp / rn_slpmax 743 641 ! 744 642 uslpml (1,:) = 0._wp ; uslpml (jpi,:) = 0._wp … … 752 650 DO ji = 1, jpi 753 651 ik = nmln(ji,jj) - 1 754 IF( jk <= ik .AND. jk >= mikt(ji,jj) ) THEN 755 omlmask(ji,jj,jk) = 1._wp 756 ELSE 757 omlmask(ji,jj,jk) = 0._wp 652 IF( jk <= ik ) THEN ; omlmask(ji,jj,jk) = 1._wp 653 ELSE ; omlmask(ji,jj,jk) = 0._wp 758 654 ENDIF 759 655 END DO … … 777 673 ! 778 674 ! !- vertical density gradient for u- and v-slopes (from dzr at T-point) 779 iku = MIN( MAX( miku(ji,jj)+1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1 ) ! ML (MAX of T-pts, bound by jpkm1)780 ikv = MIN( MAX( mikv(ji,jj)+1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) !675 iku = MIN( MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1 ) ! ML (MAX of T-pts, bound by jpkm1) 676 ikv = MIN( MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) ! 781 677 zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj ,iku) ) 782 678 zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji ,jj+1,ikv) ) 783 679 ! !- horizontal density gradient at u- & v-points 784 zau = p_gru(ji,jj,iku) /e1u(ji,jj)785 zav = p_grv(ji,jj,ikv) /e2v(ji,jj)680 zau = p_gru(ji,jj,iku) * r1_e1u(ji,jj) 681 zav = p_grv(ji,jj,ikv) * r1_e2v(ji,jj) 786 682 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 787 683 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 788 zbu = MIN( zbu , - 100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau ) )789 zbv = MIN( zbv , - 100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav ) )684 zbu = MIN( zbu , - z1_slpmax * ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau ) ) 685 zbv = MIN( zbv , - z1_slpmax * ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav ) ) 790 686 ! !- Slope at u- & v-points (uslpml, vslpml) 791 687 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku) … … 812 708 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zaj ) ) 813 709 ! !- i- & j-slope at w-points (wslpiml, wslpjml) 814 wslpiml(ji,jj) = zai / ( zbi - zeps ) * wmask (ji,jj,ik)815 wslpjml(ji,jj) = zaj / ( zbj - zeps ) * wmask (ji,jj,ik)710 wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) 711 wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik) 816 712 END DO 817 713 END DO … … 831 727 !! ** Purpose : Initialization for the isopycnal slopes computation 832 728 !! 833 !! ** Method : read the nammbf namelist and check the parameter 834 !! values called by tra_dmp at the first timestep (nit000) 729 !! ** Method : 835 730 !!---------------------------------------------------------------------- 836 731 INTEGER :: ji, jj, jk ! dummy loop indices … … 845 740 WRITE(numout,*) '~~~~~~~~~~~~' 846 741 ENDIF 847 848 IF( ln_traldf_grif ) THEN ! Griffies operator : triad of slopes 849 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 ) 850 ALLOCATE( triadi (jpi,jpj,jpk,0:1,0:1) , triadj (jpi,jpj,jpk,0:1,0:1) , STAT=ierr ) 851 IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Griffies operator slope' ) 852 ! 742 ! 743 ALLOCATE( ah_wslp2(jpi,jpj,jpk) , akz(jpi,jpj,jpk) , STAT=ierr ) 744 IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate ah_slp2 or akz' ) 745 ! 746 IF( ln_traldf_triad ) THEN ! Griffies operator : triad of slopes 747 IF(lwp) WRITE(numout,*) ' Griffies (triad) operator initialisation' 748 ALLOCATE( triadi_g(jpi,jpj,jpk,0:1,0:1) , triadj_g(jpi,jpj,jpk,0:1,0:1) , & 749 & triadi (jpi,jpj,jpk,0:1,0:1) , triadj (jpi,jpj,jpk,0:1,0:1) , & 750 & wslp2 (jpi,jpj,jpk) , STAT=ierr ) 751 IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Griffies operator slope' ) 853 752 IF( ln_dynldf_iso ) CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' ) 854 753 ! 855 754 ELSE ! Madec operator : slopes at u-, v-, and w-points 856 ALLOCATE( uslp(jpi,jpj,jpk) , vslp(jpi,jpj,jpk) , wslpi(jpi,jpj,jpk) , wslpj(jpi,jpj,jpk) , & 857 & omlmask(jpi,jpj,jpk) , uslpml(jpi,jpj) , vslpml(jpi,jpj) , wslpiml(jpi,jpj) , wslpjml(jpi,jpj) , STAT=ierr ) 755 IF(lwp) WRITE(numout,*) ' Madec operator initialisation' 756 ALLOCATE( omlmask(jpi,jpj,jpk) , & 757 & uslp(jpi,jpj,jpk) , uslpml(jpi,jpj) , wslpi(jpi,jpj,jpk) , wslpiml(jpi,jpj) , & 758 & vslp(jpi,jpj,jpk) , vslpml(jpi,jpj) , wslpj(jpi,jpj,jpk) , wslpjml(jpi,jpj) , STAT=ierr ) 858 759 IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Madec operator slope ' ) 859 760 … … 865 766 wslpj(:,:,:) = 0._wp ; wslpjml(:,:) = 0._wp 866 767 867 IF(ln_sco .AND. (ln_traldf_hor .OR. ln_dynldf_hor )) THEN 868 IF(lwp) WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 869 870 ! geopotential diffusion in s-coordinates on tracers and/or momentum 871 ! The slopes of s-surfaces are computed once (no call to ldfslp in step) 872 ! The slopes for momentum diffusion are i- or j- averaged of those on tracers 873 874 ! set the slope of diffusion to the slope of s-surfaces 875 ! ( c a u t i o n : minus sign as fsdep has positive value ) 876 DO jk = 1, jpk 877 DO jj = 2, jpjm1 878 DO ji = fs_2, fs_jpim1 ! vector opt. 879 uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk) 880 vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk) 881 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 882 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 883 END DO 884 END DO 885 END DO 886 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) ! Lateral boundary conditions 887 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 888 ENDIF 768 !!gm I no longer understand this..... 769 !!gm IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN 770 ! IF(lwp) WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 771 ! 772 ! ! geopotential diffusion in s-coordinates on tracers and/or momentum 773 ! ! The slopes of s-surfaces are computed once (no call to ldfslp in step) 774 ! ! The slopes for momentum diffusion are i- or j- averaged of those on tracers 775 ! 776 ! ! set the slope of diffusion to the slope of s-surfaces 777 ! ! ( c a u t i o n : minus sign as fsdep has positive value ) 778 ! DO jk = 1, jpk 779 ! DO jj = 2, jpjm1 780 ! DO ji = fs_2, fs_jpim1 ! vector opt. 781 ! uslp (ji,jj,jk) = - ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 782 ! vslp (ji,jj,jk) = - ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 783 ! wslpi(ji,jj,jk) = - ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * r1_e1t(ji,jj) * wmask(ji,jj,jk) * 0.5 784 ! wslpj(ji,jj,jk) = - ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * r1_e2t(ji,jj) * wmask(ji,jj,jk) * 0.5 785 ! END DO 786 ! END DO 787 ! END DO 788 ! CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) ! Lateral boundary conditions 789 ! CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 790 !!gm ENDIF 889 791 ENDIF 890 792 ! … … 892 794 ! 893 795 END SUBROUTINE ldf_slp_init 894 895 #else896 !!------------------------------------------------------------------------897 !! Dummy module : NO Rotation of lateral mixing tensor898 !!------------------------------------------------------------------------899 LOGICAL, PUBLIC, PARAMETER :: lk_ldfslp = .FALSE. !: slopes flag900 CONTAINS901 SUBROUTINE ldf_slp( kt, prd, pn2 ) ! Dummy routine902 INTEGER, INTENT(in) :: kt903 REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2904 WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1)905 END SUBROUTINE ldf_slp906 SUBROUTINE ldf_slp_grif( kt ) ! Dummy routine907 INTEGER, INTENT(in) :: kt908 WRITE(*,*) 'ldf_slp_grif: You should not have seen this print! error?', kt909 END SUBROUTINE ldf_slp_grif910 SUBROUTINE ldf_slp_init ! Dummy routine911 END SUBROUTINE ldf_slp_init912 #endif913 796 914 797 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.