Changeset 2528 for trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
- Timestamp:
- 2010-12-27T18:33:53+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
- Property svn:eol-style deleted
r2389 r2528 7 7 !! 8.0 ! 1997-06 (G. Madec) optimization, lbc 8 8 !! 8.1 ! 1999-10 (A. Jouzeau) NEW profile in the mixed layer 9 !! NEMO 0.5 ! 2002-10 (G. Madec) Free form, F90 10 !! 1.0 ! 2005-10 (A. Beckmann) correction for s-coordinates 11 !! 3.2 ! 2010-11 (F. Dupond, G. Madec) bug correction in slopes just below the ML 9 !! NEMO 1.0 ! 2002-10 (G. Madec) Free form, F90 10 !! - ! 2005-10 (A. Beckmann) correction for s-coordinates 11 !! 3.3 ! 2010-10 (G. Nurser, C. Harris, G. Madec) add Griffies operator 12 !! - ! 2010-11 (F. Dupond, G. Madec) bug correction in slopes just below the ML 12 13 !!---------------------------------------------------------------------- 13 14 #if defined key_ldfslp || defined key_esopa … … 15 16 !! 'key_ldfslp' Rotation of lateral mixing tensor 16 17 !!---------------------------------------------------------------------- 17 !! ldf_slp : compute the slopes of iso-neutral surface 18 !! ldf_slp_mxl : compute the slopes of iso-neutral surface just below the Mixed Layer 18 !! ldf_slp_grif : calculates the triads of isoneutral slopes (Griffies operator) 19 !! ldf_slp : calculates the slopes of neutral surface (Madec operator) 20 !! ldf_slp_mxl : calculates the slopes at the base of the mixed layer (Madec operator) 19 21 !! ldf_slp_init : initialization of the slopes computation 20 22 !!---------------------------------------------------------------------- 21 23 USE oce ! ocean dynamics and tracers 22 24 USE dom_oce ! ocean space and time domain 23 USE ldftra_oce 24 USE ldfdyn_oce 25 USE ldftra_oce ! lateral diffusion: traceur 26 USE ldfdyn_oce ! lateral diffusion: dynamics 25 27 USE phycst ! physical constants 26 28 USE zdfmxl ! mixed layer depth 29 USE eosbn2 ! equation of states 27 30 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 28 31 USE in_out_manager ! I/O manager … … 32 35 PRIVATE 33 36 34 PUBLIC ldf_slp ! routine called by step.F90 35 36 LOGICAL , PUBLIC, PARAMETER :: lk_ldfslp = .TRUE. !: slopes flag 37 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: uslp, wslpi !: i_slope at U- and W-points 38 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: vslp, wslpj !: j-slope at V- and W-points 39 40 REAL(wp), DIMENSION(jpi,jpj,jpk) :: omlmask ! mask of the surface mixed layer at T-pt 41 REAL(wp), DIMENSION(jpi,jpj) :: uslpml, wslpiml ! i_slope at U- and W-points just below the mixed layer 42 REAL(wp), DIMENSION(jpi,jpj) :: vslpml, wslpjml ! j_slope at V- and W-points just below the mixed layer 37 PUBLIC ldf_slp ! routine called by step.F90 38 PUBLIC ldf_slp_grif ! routine called by step.F90 39 PUBLIC ldf_slp_init ! routine called by opa.F90 40 41 LOGICAL , PUBLIC, PARAMETER :: lk_ldfslp = .TRUE. !: slopes flag 42 ! !! Madec operator 43 REAL(wp), PUBLIC, DIMENSION(:,:,:) , ALLOCATABLE :: uslp, wslpi !: i_slope at U- and W-points 44 REAL(wp), PUBLIC, DIMENSION(:,:,:) , ALLOCATABLE :: vslp, wslpj !: j-slope at V- and W-points 45 ! !! Griffies operator 46 REAL(wp), PUBLIC, DIMENSION(:,:,:) , ALLOCATABLE :: wslp2 !: wslp**2 from Griffies quarter cells 47 REAL(wp), PUBLIC, DIMENSION(:,:,:,:,:), ALLOCATABLE :: triadi_g, triadj_g !: skew flux slopes relative to geopotentials 48 REAL(wp), PUBLIC, DIMENSION(:,:,:,:,:), ALLOCATABLE :: triadi , triadj !: isoneutral slopes relative to model-coordinate 49 50 ! !! Madec operator 51 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: omlmask ! mask of the surface mixed layer at T-pt 52 REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: uslpml, wslpiml ! i_slope at U- and W-points just below the mixed layer 53 REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: vslpml, wslpjml ! j_slope at V- and W-points just below the mixed layer 54 55 REAL(wp) :: repsln = 1.e-25_wp ! tiny value used as minium of di(rho), dj(rho) and dk(rho) 43 56 44 57 !! * Substitutions 45 58 # include "domzgr_substitute.h90" 59 # include "ldftra_substitute.h90" 60 # include "ldfeiv_substitute.h90" 46 61 # include "vectopt_loop_substitute.h90" 47 62 !!---------------------------------------------------------------------- 48 !! NEMO/OPA 3. 2 , LOCEAN-IPSL (2009)63 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 49 64 !! $Id$ 50 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)65 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 51 66 !!---------------------------------------------------------------------- 52 53 67 CONTAINS 54 68 … … 58 72 !! 59 73 !! ** Purpose : Compute the slopes of neutral surface (slope of isopycnal 60 !! surfaces referenced locally) ( 'key_traldfiso').74 !! surfaces referenced locally) (ln_traldf_iso=T). 61 75 !! 62 76 !! ** Method : The slope in the i-direction is computed at U- and … … 80 94 USE oce , zgru => ua ! use ua as workspace 81 95 USE oce , zgrv => va ! use va as workspace 82 USE oce , zw y=> ta ! use ta as workspace96 USE oce , zww => ta ! use ta as workspace 83 97 USE oce , zwz => sa ! use sa as workspace 84 98 !! … … 90 104 INTEGER :: ii0, ii1, iku ! temporary integer 91 105 INTEGER :: ij0, ij1, ikv ! temporary integer 92 REAL(wp) :: zeps, zmg, zm05g, zalpha ! temporary scalars 93 REAL(wp) :: zcoef1, zcoef2, zcoef3 ! - - 94 REAL(wp) :: zcofu , zcofv , zcofw ! - - 95 REAL(wp) :: zau, zbu, zai, zbi, z1u, z1wu ! - - 96 REAL(wp) :: zav, zbv, zaj, zbj, z1v, z1wv ! 97 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zww ! 3D workspace 106 REAL(wp) :: zeps, zm1_g, zm1_2g, z1_16 ! local scalars 107 REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - - 108 REAL(wp) :: zcj, zfj, zav, zbv, zaj, zbj ! - - 109 REAL(wp) :: zck, zfk, zbw ! - - 110 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdzr ! 3D workspace 98 111 !!---------------------------------------------------------------------- 99 112 100 IF( kt == nit000 ) CALL ldf_slp_init ! initialization (first time-step only) 101 102 zeps = 1.e-20 ! Local constant initialization 103 zmg = -1.0 / grav 104 zm05g = -0.5 / grav 105 ! 106 zww(:,:,:) = 0.e0 107 zwz(:,:,:) = 0.e0 108 ! ! horizontal density gradient computation 109 DO jk = 1, jpk 113 zeps = 1.e-20_wp !== Local constant initialization ==! 114 z1_16 = 1.0_wp / 16._wp 115 zm1_g = -1.0_wp / grav 116 zm1_2g = -0.5_wp / grav 117 ! 118 zww(:,:,:) = 0._wp 119 zwz(:,:,:) = 0._wp 120 ! 121 DO jk = 1, jpk !== i- & j-gradient of density ==! 110 122 DO jj = 1, jpjm1 111 123 DO ji = 1, fs_jpim1 ! vector opt. … … 123 135 DO ji = 1, jpim1 124 136 # endif 125 iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj) ) - 1 ! last ocean level 126 ikv = MIN ( mbathy(ji,jj), mbathy(ji,jj+1) ) - 1 127 zgru(ji,jj,iku) = gru(ji,jj) 128 zgrv(ji,jj,ikv) = grv(ji,jj) 137 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 138 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 129 139 END DO 130 140 END DO 131 141 ENDIF 132 133 134 ! !* Local vertical density gradient evaluated from N^2 135 DO jk = 2, jpkm1 ! zwy = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point 136 zwy(:,:,jk) = zmg * ( prd(:,:,jk) + 1. ) * ( pn2 (:,:,jk) + pn2 (:,:,jk+1) ) & 137 & / MAX( tmask(:,:,jk) + tmask(:,:,jk+1), 1. ) 138 END DO 139 zwy(:,:,1) = 0.e0 ! surface value set to zero 140 141 142 CALL ldf_slp_mxl( prd, pn2 ) !* Slopes of isopycnal surfaces just below the mixed layer 143 142 ! 143 zdzr(:,:,1) = 0._wp !== Local vertical density gradient at T-point == ! (evaluated from N^2) 144 DO jk = 2, jpkm1 145 ! ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point 146 ! ! trick: tmask(ik ) = 0 => all pn2 = 0 => zdzr = 0 147 ! ! else tmask(ik+1) = 0 => pn2(ik+1) = 0 => zdzr divides by 1 148 ! ! umask(ik+1) /= 0 => all pn2 /= 0 => zdzr divides by 2 149 ! ! NB: 1/(tmask+1) = (1-.5*tmask) substitute a / by a * ==> faster 150 zdzr(:,:,jk) = zm1_g * ( prd(:,:,jk) + 1._wp ) & 151 & * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask(:,:,jk+1) ) 152 END DO 153 ! 154 ! !== Slopes just below the mixed layer ==! 155 CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr ) ! output: uslpml, vslpml, wslpiml, wslpjml 156 144 157 145 158 ! I. slopes at u and v point | uslp = d/di( prd ) / d/dz( prd ) … … 149 162 DO jj = 2, jpjm1 150 163 DO ji = fs_2, fs_jpim1 ! vector opt. 151 ! horizontal and vertical density gradient at u- and v-points 152 zau = 1. / e1u(ji,jj) * zgru(ji,jj,jk) 153 zav = 1. / e2v(ji,jj) * zgrv(ji,jj,jk) 154 zbu = 0.5 * ( zwy(ji,jj,jk) + zwy(ji+1,jj ,jk) ) 155 zbv = 0.5 * ( zwy(ji,jj,jk) + zwy(ji ,jj+1,jk) ) 156 ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 157 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 158 zbu = MIN( zbu, -100.*ABS( zau ), -7.e+3/fse3u(ji,jj,jk)*ABS( zau ) ) 159 zbv = MIN( zbv, -100.*ABS( zav ), -7.e+3/fse3v(ji,jj,jk)*ABS( zav ) ) 160 ! uslp and vslp output in zwz and zww, resp. 161 zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 162 zwz (ji,jj,jk) = ( ( 1. - zalpha) * zau / ( zbu - zeps ) & 163 & + zalpha * uslpml(ji,jj) & 164 & * 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) ) & 165 & / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5. ) ) * umask(ji,jj,jk) 166 zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 167 zww (ji,jj,jk) = ( ( 1. - zalpha) * zav / ( zbv - zeps ) & 168 & + zalpha * vslpml(ji,jj) & 169 & * 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) ) & 170 & / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk) 164 ! ! horizontal and vertical density gradient at u- and v-points 165 zau = zgru(ji,jj,jk) / e1u(ji,jj) 166 zav = zgrv(ji,jj,jk) / e2v(ji,jj) 167 zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj ,jk) ) 168 zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji ,jj+1,jk) ) 169 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 170 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 171 zbu = MIN( zbu, -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,jk)* ABS( zau ) ) 172 zbv = MIN( zbv, -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav ) ) 173 ! ! uslp and vslp output in zwz and zww, resp. 174 zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 175 zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 176 zwz(ji,jj,jk) = ( ( 1. - zfi) * zau / ( zbu - zeps ) & 177 & + zfi * uslpml(ji,jj) & 178 & * 0.5_wp * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) ) & 179 & / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5._wp ) ) * umask(ji,jj,jk) 180 zww(ji,jj,jk) = ( ( 1. - zfj) * zav / ( zbv - zeps ) & 181 & + zfj * vslpml(ji,jj) & 182 & * 0.5_wp * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) ) & 183 & / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk) 184 !!gm modif to suppress omlmask.... (as in Griffies case) 185 ! ! ! jk must be >= ML level for zf=1. otherwise zf=0. 186 ! zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp ) 187 ! zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp ) 188 ! zci = 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) ) 189 ! zcj = 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) ) 190 ! zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk) 191 ! zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk) 192 !!gm end modif 171 193 END DO 172 194 END DO … … 174 196 CALL lbc_lnk( zwz, 'U', -1. ) ; CALL lbc_lnk( zww, 'V', -1. ) ! lateral boundary conditions 175 197 ! 176 zcofu = 1. / 16. !* horizontal Shapiro filter 177 zcofv = 1. / 16. 198 ! !* horizontal Shapiro filter 178 199 DO jk = 2, jpkm1 179 DO jj = 2, jpjm1, jpj-3! rows jj=2 and =jpjm1 only200 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 180 201 DO ji = 2, jpim1 181 uslp(ji,jj,jk) = z cofu* ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) &202 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 182 203 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 183 204 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 184 205 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 185 206 & + 4.* zwz(ji ,jj ,jk) ) 186 vslp(ji,jj,jk) = z cofv* ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) &207 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 187 208 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 188 209 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & … … 193 214 DO jj = 3, jpj-2 ! other rows 194 215 DO ji = fs_2, fs_jpim1 ! vector opt. 195 uslp(ji,jj,jk) = z cofu* ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) &216 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 196 217 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 197 218 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 198 219 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 199 220 & + 4.* zwz(ji ,jj ,jk) ) 200 vslp(ji,jj,jk) = z cofv* ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) &221 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 201 222 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 202 223 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & … … 208 229 DO jj = 2, jpjm1 209 230 DO ji = fs_2, fs_jpim1 ! vector opt. 210 z1u = ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk) )*.5 211 z1v = ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk) )*.5 212 z1wu = ( umask(ji,jj,jk) + umask(ji,jj,jk+1) )*.5 213 z1wv = ( vmask(ji,jj,jk) + vmask(ji,jj,jk+1) )*.5 214 uslp(ji,jj,jk) = uslp(ji,jj,jk) * z1u * z1wu 215 vslp(ji,jj,jk) = vslp(ji,jj,jk) * z1v * z1wv 231 uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk ) ) * 0.5_wp & 232 & * ( umask(ji,jj ,jk) + umask(ji,jj ,jk+1) ) * 0.5_wp 233 vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk ) ) * 0.5_wp & 234 & * ( vmask(ji ,jj,jk) + vmask(ji ,jj,jk+1) ) * 0.5_wp 216 235 END DO 217 236 END DO … … 222 241 ! =========================== | wslpj = mij( d/dj( prd ) / d/dz( prd ) 223 242 ! 224 ! !* Local vertical density gradient evaluated from N^2 225 DO jk = 2, jpkm1 ! zwy = d/dz(prd)= - mk ( prd ) / grav * pn2 -- at w point 226 DO jj = 1, jpj 227 DO ji = 1, jpi 228 zwy(ji,jj,jk) = zm05g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) 229 END DO 230 END DO 231 END DO 232 DO jk = 2, jpkm1 !* Slopes at w point 243 DO jk = 2, jpkm1 233 244 DO jj = 2, jpjm1 234 245 DO ji = fs_2, fs_jpim1 ! vector opt. 235 ! ! horizontal density i-gradient at w-points236 z coef1 = MAX( zeps, umask(ji-1,jj,jk )+umask(ji,jj,jk ) &237 & +umask(ji-1,jj,jk-1)+umask(ji,jj,jk-1) )238 zcoef1 = 1. / ( zcoef1 * e1t (ji,jj) )239 z ai = zcoef1 * ( zgru(ji ,jj,jk ) + zgru(ji ,jj,jk-1)&240 & + zgru(ji-1,jj,jk-1) + zgru(ji-1,jj,jk ) ) * tmask (ji,jj,jk)241 ! ! horizontal density j-gradient at w-points242 zcoef2 = MAX( zeps, vmask(ji,jj-1,jk )+vmask(ji,jj,jk-1) &243 & +vmask(ji,jj-1,jk-1)+vmask(ji,jj,jk ) )244 zcoef2 = 1.0 / ( zcoef2 * e2t (ji,jj))245 zaj = zcoef2 * ( zgrv(ji,jj ,jk ) + zgrv(ji,jj ,jk-1)&246 & + zgrv(ji,jj-1,jk-1) + zgrv(ji,jj-1,jk ) )* tmask (ji,jj,jk)246 ! !* Local vertical density gradient evaluated from N^2 247 zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) 248 ! !* Slopes at w point 249 ! ! i- & j-gradient of density at w-points 250 zci = MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk ) & 251 & + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps ) * e1t(ji,jj) 252 zcj = MAX( vmask(ji,jj-1,jk ) + vmask(ji,jj,jk-1) & 253 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk ) , zeps ) * e2t(ji,jj) 254 zai = ( zgru (ji-1,jj,jk ) + zgru (ji,jj,jk-1) & 255 & + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk ) ) / zci * tmask (ji,jj,jk) 256 zaj = ( zgrv (ji,jj-1,jk ) + zgrv (ji,jj,jk-1) & 257 & + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk ) ) / zcj * tmask (ji,jj,jk) 247 258 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 248 ! ! static instability: kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 249 zbi = MIN( zwy (ji,jj,jk),- 100.*ABS(zai), -7.e+3/fse3w(ji,jj,jk)*ABS(zai) ) 250 zbj = MIN( zwy (ji,jj,jk), -100.*ABS(zaj), -7.e+3/fse3w(ji,jj,jk)*ABS(zaj) ) 251 ! ! wslpi and wslpj output in zwz and zww, resp. 252 zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) ) 253 zcoef3 = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10. ) 254 zwz(ji,jj,jk) = ( zai / ( zbi - zeps) * ( 1. - zalpha ) & 255 & + zcoef3 * wslpiml(ji,jj) * zalpha ) * tmask (ji,jj,jk) 256 zww(ji,jj,jk) = ( zaj / ( zbj - zeps) * ( 1. - zalpha ) & 257 & + zcoef3 * wslpjml(ji,jj) * zalpha ) * tmask (ji,jj,jk) 258 END DO 259 END DO 260 END DO 261 CALL lbc_lnk( zwz, 'T', -1. ) ; CALL lbc_lnk( zww, 'T', -1. ) ! lateral boundary conditions on zwz and zww 259 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 260 zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zai ) ) 261 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zaj ) ) 262 ! ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 263 zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) ) ! zfk=1 in the ML otherwise zfk=0 264 zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10._wp ) 265 zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk ) * tmask(ji,jj,jk) 266 zww(ji,jj,jk) = ( zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk ) * tmask(ji,jj,jk) 267 268 !!gm modif to suppress omlmask.... (as in Griffies operator) 269 ! ! ! jk must be >= ML level for zfk=1. otherwise zfk=0. 270 ! zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp ) 271 ! zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10. ) 272 ! zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk) 273 ! zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk) 274 !!gm end modif 275 END DO 276 END DO 277 END DO 278 CALL lbc_lnk( zwz, 'T', -1. ) ; CALL lbc_lnk( zww, 'T', -1. ) ! lateral boundary conditions 262 279 ! 263 280 ! !* horizontal Shapiro filter 264 281 DO jk = 2, jpkm1 265 DO jj = 2, jpjm1, jpj-3 ! rows jj=2 and =jpjm1282 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 266 283 DO ji = 2, jpim1 267 zcofw = tmask(ji,jj,jk) / 16.268 284 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 269 285 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 270 286 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 271 287 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 272 & + 4.* zwz(ji ,jj ,jk) ) * zcofw288 & + 4.* zwz(ji ,jj ,jk) ) * z1_16 * tmask(ji,jj,jk) 273 289 274 290 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & … … 276 292 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 277 293 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 278 & + 4.* zww(ji ,jj ,jk) ) * zcofw294 & + 4.* zww(ji ,jj ,jk) ) * z1_16 * tmask(ji,jj,jk) 279 295 END DO 280 296 END DO 281 297 DO jj = 3, jpj-2 ! other rows 282 298 DO ji = fs_2, fs_jpim1 ! vector opt. 283 zcofw = tmask(ji,jj,jk) / 16.284 299 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 285 300 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 286 301 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 287 302 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 288 & + 4.* zwz(ji ,jj ,jk) ) * zcofw303 & + 4.* zwz(ji ,jj ,jk) ) * z1_16 * tmask(ji,jj,jk) 289 304 290 305 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & … … 292 307 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 293 308 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 294 & + 4.* zww(ji ,jj ,jk) ) * zcofw309 & + 4.* zww(ji ,jj ,jk) ) * z1_16 * tmask(ji,jj,jk) 295 310 END DO 296 311 END DO … … 298 313 DO jj = 2, jpjm1 299 314 DO ji = fs_2, fs_jpim1 ! vector opt. 300 z 1u = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) *.5301 z1v = ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) *.5302 wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * z 1u * z1v303 wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * z 1u * z1v315 zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) & 316 & * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 317 wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck 318 wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck 304 319 END DO 305 320 END DO 306 321 END DO 307 322 308 309 323 ! III. Specific grid points 310 324 ! =========================== … … 313 327 ! ! Gibraltar Strait 314 328 ij0 = 50 ; ij1 = 53 315 ii0 = 69 ; ii1 = 71 ; uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0. e0329 ii0 = 69 ; ii1 = 71 ; uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 316 330 ij0 = 51 ; ij1 = 53 317 ii0 = 68 ; ii1 = 71 ; vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0. e0318 ii0 = 69 ; ii1 = 71 ; wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0. e0319 ii0 = 69 ; ii1 = 71 ; wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0. e0331 ii0 = 68 ; ii1 = 71 ; vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 332 ii0 = 69 ; ii1 = 71 ; wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 333 ii0 = 69 ; ii1 = 71 ; wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 320 334 ! 321 335 ! ! Mediterrannean Sea 322 336 ij0 = 49 ; ij1 = 56 323 ii0 = 71 ; ii1 = 90 ; uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0. e0337 ii0 = 71 ; ii1 = 90 ; uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 324 338 ij0 = 50 ; ij1 = 56 325 ii0 = 70 ; ii1 = 90 ; vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0. e0326 ii0 = 71 ; ii1 = 90 ; wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0. e0327 ii0 = 71 ; ii1 = 90 ; wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0. e0339 ii0 = 70 ; ii1 = 90 ; vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 340 ii0 = 71 ; ii1 = 90 ; wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 341 ii0 = 71 ; ii1 = 90 ; wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 328 342 ENDIF 329 343 … … 341 355 ! 342 356 END SUBROUTINE ldf_slp 343 344 345 SUBROUTINE ldf_slp_mxl( prd, pn2 ) 357 358 359 SUBROUTINE ldf_slp_grif ( kt ) 360 !!---------------------------------------------------------------------- 361 !! *** ROUTINE ldf_slp_grif *** 362 !! 363 !! ** Purpose : Compute the squared slopes of neutral surfaces (slope 364 !! of iso-pycnal surfaces referenced locally) (ln_traldf_grif=T) 365 !! at W-points using the Griffies quarter-cells. 366 !! 367 !! ** Method : calculates alpha and beta at T-points 368 !! 369 !! ** Action : - triadi_g, triadj_g T-pts i- and j-slope triads relative to geopot. (used for eiv) 370 !! - triadi , triadj T-pts i- and j-slope triads relative to model-coordinate 371 !! - wslp2 squared slope of neutral surfaces at w-points. 372 !!---------------------------------------------------------------------- 373 USE oce, zdit => ua ! use ua as workspace 374 USE oce, zdis => va ! use va as workspace 375 USE oce, zdjt => ta ! use ta as workspace 376 USE oce, zdjs => sa ! use sa as workspace 377 !! 378 INTEGER, INTENT( in ) :: kt ! ocean time-step index 379 !! 380 INTEGER :: ji, jj, jk, jl, ip, jp, kp ! dummy loop indices 381 INTEGER :: iku, ikv ! temporary integer 382 REAL(wp) :: zfacti, zfactj, zatempw,zatempu,zatempv ! local scalars 383 REAL(wp) :: zbu, zbv, zbti, zbtj 384 REAL(wp) :: zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_lim2, zti_g_raw, zti_g_lim 385 REAL(wp) :: zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_lim2, ztj_g_raw, ztj_g_lim 386 REAL(wp) :: zdzrho_raw 387 REAL(wp), DIMENSION(jpi,jpj,jpk,0:1) :: zdzrho, zdyrho, zdxrho ! Horizontal and vertical density gradients 388 REAL(wp), DIMENSION(jpi,jpj,0:1,0:1) :: zti_mlb, ztj_mlb 389 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdkt, zdks 390 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zalpha, zbeta ! alpha, beta at T points, at depth fsgdept 391 REAL(wp), DIMENSION(jpi,jpj) :: z1_mlbw 392 !!---------------------------------------------------------------------- 393 394 !--------------------------------! 395 ! Some preliminary calculation ! 396 !--------------------------------! 397 ! 398 CALL eos_alpbet( tsb, zalpha, zbeta ) !== before thermal and haline expension coeff. at T-points ==! 399 ! 400 DO jk = 1, jpkm1 !== before lateral T & S gradients at T-level jk ==! 401 DO jj = 1, jpjm1 402 DO ji = 1, fs_jpim1 ! vector opt. 403 zdit(ji,jj,jk) = ( tb(ji+1,jj,jk) - tb(ji,jj,jk) ) * umask(ji,jj,jk) ! i-gradient of T and S at jj 404 zdis(ji,jj,jk) = ( sb(ji+1,jj,jk) - sb(ji,jj,jk) ) * umask(ji,jj,jk) 405 zdjt(ji,jj,jk) = ( tb(ji,jj+1,jk) - tb(ji,jj,jk) ) * vmask(ji,jj,jk) ! j-gradient of T and S at jj 406 zdjs(ji,jj,jk) = ( sb(ji,jj+1,jk) - sb(ji,jj,jk) ) * vmask(ji,jj,jk) 407 END DO 408 END DO 409 END DO 410 IF( ln_zps ) THEN ! partial steps: correction at the last level 411 # if defined key_vectopt_loop 412 DO jj = 1, 1 413 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 414 # else 415 DO jj = 1, jpjm1 416 DO ji = 1, jpim1 417 # endif 418 zdit(ji,jj,mbku(ji,jj)) = gtsu(ji,jj,jp_tem) ! i-gradient of T and S 419 zdis(ji,jj,mbku(ji,jj)) = gtsu(ji,jj,jp_sal) 420 zdjt(ji,jj,mbkv(ji,jj)) = gtsv(ji,jj,jp_tem) ! j-gradient of T and S 421 zdjs(ji,jj,mbkv(ji,jj)) = gtsv(ji,jj,jp_sal) 422 END DO 423 END DO 424 ENDIF 425 ! 426 zdkt(:,:,1) = 0._wp !== before vertical T & S gradient at w-level ==! 427 zdks(:,:,1) = 0._wp 428 DO jk = 2, jpk 429 zdkt(:,:,jk) = ( tb(:,:,jk-1) - tb(:,:,jk) ) * tmask(:,:,jk) 430 zdks(:,:,jk) = ( sb(:,:,jk-1) - sb(:,:,jk) ) * tmask(:,:,jk) 431 END DO 432 ! 433 ! 434 DO jl = 0, 1 !== density i-, j-, and k-gradients ==! 435 ip = jl ; jp = jl ! guaranteed nonzero gradients ( absolute value larger than repsln) 436 DO jk = 1, jpkm1 ! done each pair of triad 437 DO jj = 1, jpjm1 ! NB: not masked due to the minimum value set 438 DO ji = 1, fs_jpim1 ! vector opt. 439 zdxrho_raw = ( zalpha(ji+ip,jj ,jk) * zdit(ji,jj,jk) + zbeta(ji+ip,jj ,jk) * zdis(ji,jj,jk) ) / e1u(ji,jj) 440 zdyrho_raw = ( zalpha(ji ,jj+jp,jk) * zdjt(ji,jj,jk) + zbeta(ji ,jj+jp,jk) * zdjs(ji,jj,jk) ) / e2v(ji,jj) 441 zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 442 zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 443 END DO 444 END DO 445 END DO 446 END DO 447 DO kp = 0, 1 !== density i-, j-, and k-gradients ==! 448 DO jk = 1, jpkm1 ! done each pair of triad 449 DO jj = 1, jpj ! NB: not masked due to the minimum value set 450 DO ji = 1, jpi ! vector opt. 451 zdzrho_raw = ( zalpha(ji,jj,jk) * zdkt(ji,jj,jk+kp) + zbeta(ji,jj,jk) * zdks(ji,jj,jk+kp) ) & 452 & / fse3w(ji,jj,jk+kp) 453 zdzrho(ji ,jj ,jk, kp) = - MIN( - repsln, zdzrho_raw ) ! force zdzrho >= repsln 454 END DO 455 END DO 456 END DO 457 END DO 458 ! 459 DO jj = 1, jpj !== Reciprocal depth of the w-point below ML base ==! 460 DO ji = 1, jpi 461 jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1 ! MIN in case ML depth is the ocean depth 462 z1_mlbw(ji,jj) = 1._wp / fsdepw(ji,jj,jk) 463 END DO 464 END DO 465 ! 466 ! !== intialisations to zero ==! 467 ! 468 wslp2 (:,:,:) = 0._wp ! wslp2 will be cumulated 3D field set to zero 469 triadi_g(:,:,1,:,:) = 0._wp ; triadi_g(:,:,jpk,:,:) = 0._wp ! set surface and bottom slope to zero 470 triadj_g(:,:,1,:,:) = 0._wp ; triadj_g(:,:,jpk,:,:) = 0._wp 471 !!gm _iso set to zero missing 472 triadi (:,:,1,:,:) = 0._wp ; triadj (:,:,jpk,:,:) = 0._wp ! set surface and bottom slope to zero 473 triadj (:,:,1,:,:) = 0._wp ; triadj (:,:,jpk,:,:) = 0._wp 474 475 !-------------------------------------! 476 ! Triads just below the Mixed Layer ! 477 !-------------------------------------! 478 ! 479 DO jl = 0, 1 ! calculate slope of the 4 triads immediately ONE level below mixed-layer base 480 DO kp = 0, 1 ! with only the slope-max limit and MASKED 481 DO jj = 1, jpjm1 482 DO ji = 1, fs_jpim1 483 ip = jl ; jp = jl 484 jk = MIN( nmln(ji+ip,jj) , mbkt(ji+ip,jj) ) + 1 ! ML level+1 (MIN in case ML depth is the ocean 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) ) / e1u(ji,jj) ) * umask(ji,jj,jk) 487 jk = MIN( nmln(ji,jj+jp) , mbkt(ji,jj+jp) ) + 1 488 ztj_g_raw = ( zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp) & 489 & + ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) 490 zti_mlb(ji+ip,jj ,1-ip,kp) = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 491 ztj_mlb(ji ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) 492 END DO 493 END DO 494 END DO 495 END DO 496 497 !-------------------------------------! 498 ! Triads with surface limits ! 499 !-------------------------------------! 500 ! 501 DO kp = 0, 1 ! k-index of triads 502 DO jl = 0, 1 503 ip = jl ; jp = jl ! i- and j-indices of triads (i-k and j-k planes) 504 DO jk = 1, jpkm1 505 DO jj = 1, jpjm1 506 DO ji = 1, fs_jpim1 ! vector opt. 507 ! 508 ! Calculate slope relative to geopotentials used for GM skew fluxes 509 ! For s-coordinate, subtract slope at t-points (equivalent to *adding* gradient of depth) 510 ! Limit by slope *relative to geopotentials* by rn_slpmax, and mask by psi-point 511 ! masked by umask taken at the level of dz(rho) 512 ! 513 ! raw slopes: unmasked unbounded slopes (relative to geopotential (zti_g) and model surface (zti) 514 ! 515 zti_raw = zdxrho(ji+ip,jj ,jk,1-ip) / zdzrho(ji+ip,jj ,jk,kp) ! unmasked 516 ztj_raw = zdyrho(ji ,jj+jp,jk,1-jp) / zdzrho(ji ,jj+jp,jk,kp) 517 zti_coord = ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj) 518 ztj_coord = ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) 519 ! unmasked 520 zti_g_raw = zti_raw + zti_coord ! ref to geopot surfaces 521 ztj_g_raw = ztj_raw + ztj_coord 522 zti_g_lim = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 523 ztj_g_lim = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) 524 ! 525 ! Below ML use limited zti_g as is 526 ! Inside ML replace by linearly reducing sx_mlb towards surface 527 ! 528 zfacti = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji+ip,jj)), wp ) ! k index of uppermost point(s) of triad is jk+kp-1 529 zfactj = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp)), wp ) ! must be .ge. nmln(ji,jj) for zfact=1 530 ! ! otherwise zfact=0 531 zti_g_lim = zfacti * zti_g_lim & 532 & + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp) & 533 & * fsdepw(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj) 534 ztj_g_lim = zfactj * ztj_g_lim & 535 & + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp) & 536 & * fsdepw(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) ! masked 537 ! 538 triadi_g(ji+ip,jj ,jk,1-ip,kp) = zti_g_lim * umask(ji,jj,jk+kp) 539 triadj_g(ji ,jj+jp,jk,1-jp,kp) = ztj_g_lim * vmask(ji,jj,jk+kp) 540 ! 541 ! Get coefficients of isoneutral diffusion tensor 542 ! 1. Utilise gradients *relative* to s-coordinate, so add t-point slopes (*subtract* depth gradients) 543 ! 2. We require that isoneutral diffusion gives no vertical buoyancy flux 544 ! i.e. 33 term = (real slope* 31, 13 terms) 545 ! To do this, retain limited sx**2 in vertical flux, but divide by real slope for 13/31 terms 546 ! Equivalent to tapering A_iso = sx_limited**2/(real slope)**2 547 ! 548 zti_lim = zti_g_lim - zti_coord ! remove the coordinate slope ==> relative to coordinate surfaces 549 ztj_lim = ztj_g_lim - ztj_coord 550 zti_lim2 = zti_lim * zti_lim * umask(ji,jj,jk+kp) ! square of limited slopes ! masked <<== 551 ztj_lim2 = ztj_lim * ztj_lim * vmask(ji,jj,jk+kp) 552 ! 553 zbu = e1u(ji ,jj) * e2u(ji ,jj) * fse3u(ji ,jj,jk ) 554 zbv = e1v(ji ,jj) * e2v(ji ,jj) * fse3v(ji ,jj,jk ) 555 zbti = e1t(ji+ip,jj) * e2t(ji+ip,jj) * fse3w(ji+ip,jj,jk+kp) 556 zbtj = e1t(ji,jj+jp) * e2t(ji,jj+jp) * fse3w(ji,jj+jp,jk+kp) 557 ! 558 triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim2 / zti_raw ! masked 559 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim2 / ztj_raw 560 ! 561 !!gm this may inhibit vectorization on Vect Computers, and even on scalar computers.... ==> to be checked 562 wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_lim2 ! masked 563 wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_lim2 564 END DO 565 END DO 566 END DO 567 END DO 568 END DO 569 ! 570 wslp2(:,:,1) = 0._wp ! force the surface wslp to zero 571 572 CALL lbc_lnk( wslp2, 'W', 1. ) ! lateral boundary confition on wslp2 only ==>>> gm : necessary ? to be checked 573 ! 574 END SUBROUTINE ldf_slp_grif 575 576 577 SUBROUTINE ldf_slp_mxl( prd, pn2, p_gru, p_grv, p_dzr ) 346 578 !!---------------------------------------------------------------------- 347 579 !! *** ROUTINE ldf_slp_mxl *** … … 350 582 !! the mixed layer. 351 583 !! 352 !! ** Method : caution: use zgru, zgrv and zwy computed in ldf_slp 584 !! ** Method : The slope in the i-direction is computed at u- & w-points 585 !! (uslpml, wslpiml) and the slope in the j-direction is computed 586 !! at v- and w-points (vslpml, wslpjml) with the same bounds as 587 !! in ldf_slp. 353 588 !! 354 589 !! ** Action : uslpml, wslpiml : i- & j-slopes of neutral surfaces … … 356 591 !! omlmask : mixed layer mask 357 592 !!---------------------------------------------------------------------- 358 USE oce , zgru => ua ! zgru, zgrv (ua, va used as workspace) 359 USE oce , zgrv => va ! set to i- & j-density gradient in ldf_slp 360 USE oce , zwy => ta ! set to vertical density gradient at T-point in ldfslp 361 !! 362 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: prd ! in situ density 363 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: pn2 ! Brunt-Vaisala frequency (locally ref.) 364 !! 365 INTEGER :: ji, jj, jk ! dummy loop indices 366 INTEGER :: ik, ikm1 ! temporary integers 367 REAL(wp) :: zeps, zmg, zm05g ! temporary scalars 368 REAL(wp) :: zcoef1, zcoef2 ! - - 369 REAL(wp) :: zau, zbu, zai, zbi ! - - 370 REAL(wp) :: zav, zbv, zaj, zbj ! - - 371 REAL(wp) :: zbw ! - - 372 !!---------------------------------------------------------------------- 373 374 zeps = 1.e-20 ! Local constant initialization 375 zmg = -1.0 / grav 376 zm05g = -0.5 / grav 377 ! 378 uslpml (1,:) = 0.e0 ; uslpml (jpi,:) = 0.e0 379 vslpml (1,:) = 0.e0 ; vslpml (jpi,:) = 0.e0 380 wslpiml(1,:) = 0.e0 ; wslpiml(jpi,:) = 0.e0 381 wslpjml(1,:) = 0.e0 ; wslpjml(jpi,:) = 0.e0 382 383 ! ! surface mixed layer mask 384 DO jk = 1, jpk ! =1 inside the mixed layer, =0 otherwise 593 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: prd ! in situ density 594 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: pn2 ! Brunt-Vaisala frequency (locally ref.) 595 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: p_gru, p_grv ! i- & j-gradient of density (u- & v-pts) 596 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: p_dzr ! z-gradient of density (T-point) 597 !! 598 INTEGER :: ji , jj , jk ! dummy loop indices 599 INTEGER :: iku, ikv, ik, ikm1 ! local integers 600 REAL(wp) :: zeps, zm1_g, zm1_2g ! local scalars 601 REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - - 602 REAL(wp) :: zcj, zfj, zav, zbv, zaj, zbj ! - - 603 REAL(wp) :: zck, zfk, zbw ! - - 604 !!---------------------------------------------------------------------- 605 606 zeps = 1.e-20_wp !== Local constant initialization ==! 607 zm1_g = -1.0_wp / grav 608 zm1_2g = -0.5_wp / grav 609 ! 610 uslpml (1,:) = 0._wp ; uslpml (jpi,:) = 0._wp 611 vslpml (1,:) = 0._wp ; vslpml (jpi,:) = 0._wp 612 wslpiml(1,:) = 0._wp ; wslpiml(jpi,:) = 0._wp 613 wslpjml(1,:) = 0._wp ; wslpjml(jpi,:) = 0._wp 614 ! 615 ! !== surface mixed layer mask ! 616 DO jk = 1, jpk ! =1 inside the mixed layer, =0 otherwise 385 617 # if defined key_vectopt_loop 386 618 DO jj = 1, 1 … … 391 623 # endif 392 624 ik = nmln(ji,jj) - 1 393 IF( jk <= ik ) THEN ; omlmask(ji,jj,jk) = 1. e0394 ELSE ; omlmask(ji,jj,jk) = 0. e0625 IF( jk <= ik ) THEN ; omlmask(ji,jj,jk) = 1._wp 626 ELSE ; omlmask(ji,jj,jk) = 0._wp 395 627 ENDIF 396 628 END DO … … 409 641 !----------------------------------------------------------------------- 410 642 ! 411 ! !* Slope at u points412 643 # if defined key_vectopt_loop 413 644 DO jj = 1, 1 … … 417 648 DO ji = 2, jpim1 418 649 # endif 419 ! horizontal and vertical density gradient at u-points 420 ik = MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) 421 ik = MIN( ik, jpkm1 ) 422 zau = 1./ e1u(ji,jj) * zgru(ji,jj,ik) 423 zbu = 0.5*( zwy(ji,jj,ik) + zwy(ji+1,jj,ik) ) 424 ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 425 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 426 zbu = MIN( zbu, -100.*ABS(zau), -7.e+3/fse3u(ji,jj,ik)*ABS(zau) ) 427 ! uslpml 428 uslpml (ji,jj) = zau / ( zbu - zeps ) * umask (ji,jj,ik) 429 END DO 430 END DO 431 CALL lbc_lnk( uslpml, 'U', -1. ) ! lateral boundary conditions (i-gradient => sign change) 432 433 ! !* Slope at v points 434 # if defined key_vectopt_loop 435 DO jj = 1, 1 436 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling) 437 # else 438 DO jj = 2, jpjm1 439 DO ji = 2, jpim1 440 # endif 441 ! horizontal and vertical density gradient at v-points 442 ik = MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) 443 ik = MIN( ik,jpkm1 ) 444 zav = 1./ e2v(ji,jj) * zgrv(ji,jj,ik) 445 zbv = 0.5*( zwy(ji,jj,ik) + zwy(ji,jj+1,ik) ) 446 ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 447 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 448 zbv = MIN( zbv, -100.*ABS(zav), -7.e+3/fse3v(ji,jj,ik)*ABS( zav ) ) 449 ! vslpml 450 vslpml (ji,jj) = zav / ( zbv - zeps ) * vmask (ji,jj,ik) 451 END DO 452 END DO 453 CALL lbc_lnk( vslpml, 'V', -1. ) ! lateral boundary conditions (j-gradient => sign change) 454 455 ! !* Slopes at w points 456 # if defined key_vectopt_loop 457 DO jj = 1, 1 458 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling) 459 # else 460 DO jj = 2, jpjm1 461 DO ji = 2, jpim1 462 # endif 463 ik = nmln(ji,jj) + 1 464 ik = MIN( ik, jpk ) 465 ikm1 = MAX ( 1, ik-1 ) 466 ! horizontal density i-gradient at w-points 467 zcoef1 = MAX( zeps, umask(ji-1,jj,ik )+umask(ji,jj,ik ) & 468 & +umask(ji-1,jj,ikm1)+umask(ji,jj,ikm1) ) 469 zcoef1 = 1. / ( zcoef1 * e1t (ji,jj) ) 470 zai = zcoef1 * ( zgru(ji ,jj,ik ) + zgru(ji ,jj,ikm1) & 471 & + zgru(ji-1,jj,ikm1) + zgru(ji-1,jj,ik ) ) * tmask (ji,jj,ik) 472 ! horizontal density j-gradient at w-points 473 zcoef2 = MAX( zeps, vmask(ji,jj-1,ik )+vmask(ji,jj,ikm1) & 474 & +vmask(ji,jj-1,ikm1)+vmask(ji,jj,ik ) ) 475 zcoef2 = 1.0 / ( zcoef2 * e2t (ji,jj) ) 476 zaj = zcoef2 * ( zgrv(ji,jj ,ik ) + zgrv(ji,jj ,ikm1) & 477 & + zgrv(ji,jj-1,ikm1) + zgrv(ji,jj-1,ik ) ) * tmask (ji,jj,ik) 478 ! vertical density gradient at w-point (from N^2) 479 zbw = zm05g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. ) 480 ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 481 ! static instability: 482 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 483 zbi = MIN ( zbw,- 100.*ABS(zai), -7.e+3/fse3w(ji,jj,ik)*ABS(zai) ) 484 zbj = MIN ( zbw, -100.*ABS(zaj), -7.e+3/fse3w(ji,jj,ik)*ABS(zaj) ) 485 ! wslpiml and wslpjml 486 wslpiml (ji,jj) = zai / ( zbi - zeps) * tmask (ji,jj,ik) 487 wslpjml (ji,jj) = zaj / ( zbj - zeps) * tmask (ji,jj,ik) 488 END DO 489 END DO 490 CALL lbc_lnk( wslpiml, 'W', -1. ) ; CALL lbc_lnk( wslpjml, 'W', -1. ) ! lateral boundary conditions 650 ! !== Slope at u- & v-points just below the Mixed Layer ==! 651 ! 652 ! !- vertical density gradient for u- and v-slopes (from dzr at T-point) 653 iku = MIN( MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1 ) ! ML (MAX of T-pts, bound by jpkm1) 654 ikv = MIN( MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) ! 655 zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj ,iku) ) 656 zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji ,jj+1,ikv) ) 657 ! !- horizontal density gradient at u- & v-points 658 zau = p_gru(ji,jj,iku) / e1u(ji,jj) 659 zav = p_grv(ji,jj,ikv) / e2v(ji,jj) 660 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 661 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 662 zbu = MIN( zbu , -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,ik)* ABS( zau ) ) 663 zbv = MIN( zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ik)* ABS( zav ) ) 664 ! !- Slope at u- & v-points (uslpml, vslpml) 665 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,ik) 666 vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ik) 667 ! 668 ! !== i- & j-slopes at w-points just below the Mixed Layer ==! 669 ! 670 ik = MIN( nmln(ji,jj) + 1, jpk ) 671 ikm1 = MAX( 1, ik-1 ) 672 ! !- vertical density gradient for w-slope (from N^2) 673 zbw = zm1_2g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. ) 674 ! !- horizontal density i- & j-gradient at w-points 675 zci = MAX( umask(ji-1,jj,ik ) + umask(ji,jj,ik ) & 676 & + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps ) * e1t(ji,jj) 677 zcj = MAX( vmask(ji,jj-1,ik ) + vmask(ji,jj,ik ) & 678 & + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps ) * e2t(ji,jj) 679 zai = ( p_gru(ji-1,jj,ik ) + p_gru(ji,jj,ik) & 680 & + p_gru(ji-1,jj,ikm1) + p_gru(ji,jj,ikm1 ) ) / zci * tmask(ji,jj,ik) 681 zaj = ( p_grv(ji,jj-1,ik ) + p_grv(ji,jj,ik ) & 682 & + p_grv(ji,jj-1,ikm1) + p_grv(ji,jj,ikm1) ) / zcj * tmask(ji,jj,ik) 683 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0. 684 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 685 zbi = MIN( zbw , -100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zai ) ) 686 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zaj ) ) 687 ! !- i- & j-slope at w-points (wslpiml, wslpjml) 688 wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) 689 wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik) 690 END DO 691 END DO 692 !!gm this lbc_lnk should be useless.... 693 CALL lbc_lnk( uslpml , 'U', -1. ) ; CALL lbc_lnk( vslpml , 'V', -1. ) ! lateral boundary cond. (sign change) 694 CALL lbc_lnk( wslpiml, 'W', -1. ) ; CALL lbc_lnk( wslpjml, 'W', -1. ) ! lateral boundary conditions 491 695 ! 492 696 END SUBROUTINE ldf_slp_mxl … … 503 707 !!---------------------------------------------------------------------- 504 708 INTEGER :: ji, jj, jk ! dummy loop indices 709 INTEGER :: ierr ! local integer 505 710 !!---------------------------------------------------------------------- 506 711 507 712 IF(lwp) THEN 508 713 WRITE(numout,*) 509 WRITE(numout,*) 'ldf_slp : direction of lateral mixing'510 WRITE(numout,*) '~~~~~~~ '714 WRITE(numout,*) 'ldf_slp_init : direction of lateral mixing' 715 WRITE(numout,*) '~~~~~~~~~~~~' 511 716 ENDIF 512 513 ! Direction of lateral diffusion (tracers and/or momentum) 514 ! ------------------------------ 515 ! set the slope to zero (even in s-coordinates) 516 517 uslp (:,:,:) = 0.e0 518 vslp (:,:,:) = 0.e0 519 wslpi(:,:,:) = 0.e0 520 wslpj(:,:,:) = 0.e0 521 522 uslpml (:,:) = 0.e0 523 vslpml (:,:) = 0.e0 524 wslpiml(:,:) = 0.e0 525 wslpjml(:,:) = 0.e0 526 527 IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN 528 IF(lwp) THEN 529 WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 717 718 IF( ln_traldf_grif ) THEN ! Griffies operator : triad of slopes 719 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 ) 720 ALLOCATE( triadi (jpi,jpj,jpk,0:1,0:1) , triadj (jpi,jpj,jpk,0:1,0:1) , STAT=ierr ) 721 IF( ierr > 0 ) THEN 722 CALL ctl_stop( 'ldf_slp_init : unable to allocate Griffies operator slope ' ) ; RETURN 530 723 ENDIF 531 532 ! geopotential diffusion in s-coordinates on tracers and/or momentum 533 ! The slopes of s-surfaces are computed once (no call to ldfslp in step) 534 ! The slopes for momentum diffusion are i- or j- averaged of those on tracers 535 536 ! set the slope of diffusion to the slope of s-surfaces 537 ! ( c a u t i o n : minus sign as fsdep has positive value ) 538 DO jk = 1, jpk 539 DO jj = 2, jpjm1 540 DO ji = fs_2, fs_jpim1 ! vector opt. 541 uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk) 542 vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk) 543 wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5 544 wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5 724 ! 725 IF( ln_dynldf_iso ) CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' ) 726 ! 727 IF( ( ln_traldf_hor .AND. ln_dynldf_hor ) .AND. ln_sco ) & 728 & CALL ctl_stop( 'ldf_slp_init: horizontal Griffies operator ', & 729 & 'in s-coordinate not supported' ) 730 ! 731 ELSE ! Madec operator : slopes at u-, v-, and w-points 732 ALLOCATE( uslp(jpi,jpj,jpk) , vslp(jpi,jpj,jpk) , wslpi(jpi,jpj,jpk) , wslpj(jpi,jpj,jpk) , & 733 & omlmask(jpi,jpj,jpk) , uslpml(jpi,jpj) , vslpml(jpi,jpj) , wslpiml(jpi,jpj) , wslpjml(jpi,jpj) , STAT=ierr ) 734 IF( ierr > 0 ) THEN 735 CALL ctl_stop( 'ldf_slp_init : unable to allocate Madec operator slope ' ) ; RETURN 736 ENDIF 737 738 ! Direction of lateral diffusion (tracers and/or momentum) 739 ! ------------------------------ 740 uslp (:,:,:) = 0._wp ; uslpml (:,:) = 0._wp ! set the slope to zero (even in s-coordinates) 741 vslp (:,:,:) = 0._wp ; vslpml (:,:) = 0._wp 742 wslpi(:,:,:) = 0._wp ; wslpiml(:,:) = 0._wp 743 wslpj(:,:,:) = 0._wp ; wslpjml(:,:) = 0._wp 744 745 !!gm I no longer understand this..... 746 IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN 747 IF(lwp) THEN 748 WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 749 ENDIF 750 751 ! geopotential diffusion in s-coordinates on tracers and/or momentum 752 ! The slopes of s-surfaces are computed once (no call to ldfslp in step) 753 ! The slopes for momentum diffusion are i- or j- averaged of those on tracers 754 755 ! set the slope of diffusion to the slope of s-surfaces 756 ! ( c a u t i o n : minus sign as fsdep has positive value ) 757 DO jk = 1, jpk 758 DO jj = 2, jpjm1 759 DO ji = fs_2, fs_jpim1 ! vector opt. 760 uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk) 761 vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk) 762 wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5 763 wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5 764 END DO 545 765 END DO 546 766 END DO 547 END DO 548 ! Lateral boundary conditions on the slopes 549 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) 550 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 551 ENDIF 552 ! 767 ! Lateral boundary conditions on the slopes 768 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) 769 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 770 ENDIF 771 ENDIF ! 553 772 END SUBROUTINE ldf_slp_init 554 773 … … 564 783 WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1) 565 784 END SUBROUTINE ldf_slp 785 SUBROUTINE ldf_slp_init ! Dummy routine 786 END SUBROUTINE ldf_slp_init 566 787 #endif 567 788
Note: See TracChangeset
for help on using the changeset viewer.
