Changeset 2528 for trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90
- Timestamp:
- 2010-12-27T18:33:53+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90
- Property svn:eol-style deleted
r1152 r2528 2 2 !!============================================================================== 3 3 !! *** MODULE traldf_bilapg *** 4 !! Ocean active tracers: horizontal component of the lateral tracer mixing trend 4 !! Ocean tracers: horizontal component of the lateral tracer mixing trend 5 !!============================================================================== 6 !! History : 8.0 ! 1997-07 (G. Madec) Original code 7 !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module 8 !! 3.3 ! 2010-06 (C. Ethe, G. Madec) Merge TRA-TRC 5 9 !!============================================================================== 6 10 #if defined key_ldfslp || defined key_esopa … … 12 16 !! ldfght : ??? 13 17 !!---------------------------------------------------------------------- 14 !! * Modules used15 18 USE oce ! ocean dynamics and tracers variables 16 19 USE dom_oce ! ocean space and time domain variables 17 20 USE ldftra_oce ! ocean active tracers: lateral physics 18 USE trdmod ! ocean active tracers trends19 USE trdmod_oce ! ocean variables trends20 21 USE in_out_manager ! I/O manager 21 22 USE ldfslp ! iso-neutral slopes available 22 23 USE lbclnk ! ocean lateral boundary condition (or mpp link) 23 24 USE diaptr ! poleward transport diagnostics 24 USE prtctl ! Print control25 USE trc_oce ! share passive tracers/Ocean variables 25 26 26 27 IMPLICIT NONE 27 28 PRIVATE 28 29 29 !! * Routine accessibility 30 PUBLIC tra_ldf_bilapg ! routine called by step.F90 30 PUBLIC tra_ldf_bilapg ! routine called by step.F90 31 31 32 32 !! * Substitutions … … 35 35 # include "ldfeiv_substitute.h90" 36 36 !!---------------------------------------------------------------------- 37 !! OPA 9.0 , LOCEAN-IPSL (2005) 38 !! $Id$ 39 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 40 !!---------------------------------------------------------------------- 41 37 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 38 !! $Id$ 39 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 40 !!---------------------------------------------------------------------- 42 41 CONTAINS 43 42 44 SUBROUTINE tra_ldf_bilapg( kt )43 SUBROUTINE tra_ldf_bilapg( kt, cdtype, ptb, pta, kjpt ) 45 44 !!---------------------------------------------------------------------- 46 45 !! *** ROUTINE tra_ldf_bilapg *** 47 46 !! 48 !! ** Purpose : Compute the before horizontal tracer (t & s)diffusive47 !! ** Purpose : Compute the before horizontal tracer diffusive 49 48 !! trend and add it to the general trend of tracer equation. 50 49 !! … … 54 53 !! computed in routine inildf. 55 54 !! -1- compute the geopotential harmonic operator applied to 56 !! (tb,sb)and multiply it by the eddy diffusivity coefficient57 !! (done by a call to ldfght routine, result in (wk1,wk2)arrays).55 !! ptb and multiply it by the eddy diffusivity coefficient 56 !! (done by a call to ldfght routine, result in wk1 arrays). 58 57 !! Applied the domain lateral boundary conditions by call to lbc_lnk 59 58 !! -2- compute the geopotential harmonic operator applied to 60 !! (wk1,wk2) by a second call to ldfght routine (result in (wk3,wk4)59 !! wk1 by a second call to ldfght routine (result in wk2) 61 60 !! arrays). 62 !! -3- Add this trend to the general trend (ta,sa):63 !! (ta,sa) = (ta,sa) + (wk3,wk4)64 !! 65 !! ** Action : - Update (ta,sa)arrays with the before geopotential61 !! -3- Add this trend to the general trend 62 !! pta = pta + wk2 63 !! 64 !! ** Action : - Update pta arrays with the before geopotential 66 65 !! biharmonic mixing trend. 67 !! 68 !! History : 69 !! 8.0 ! 97-07 (G. Madec) Original code 70 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 71 !! 9.0 ! 04-08 (C. Talandier) New trends organization 72 !!---------------------------------------------------------------------- 73 !! * Modules used 74 USE oce , wk1 => ua, & ! use ua as workspace 75 & wk2 => va ! use va as workspace 76 77 !! * Arguments 78 INTEGER, INTENT( in ) :: kt ! ocean time-step index 79 80 !! * Local declarations 81 INTEGER :: ji, jj, jk ! dummy loop indices 82 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 83 wk3, wk4 ! work array used for rotated biharmonic 84 ! ! operator on tracers and/or momentum 85 !!---------------------------------------------------------------------- 86 87 IF( kt == nit000 ) THEN 66 !!---------------------------------------------------------------------- 67 INTEGER , INTENT(in ) :: kt ! ocean time-step index 68 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 69 INTEGER , INTENT(in ) :: kjpt ! number of tracers 70 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 71 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 72 !! 73 INTEGER :: ji, jj, jk, jn ! dummy loop indices 74 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt) :: wk1, wk2 ! 4D workspace 75 !!---------------------------------------------------------------------- 76 77 IF( kt == nit000 ) THEN 88 78 IF(lwp) WRITE(numout,*) 89 IF(lwp) WRITE(numout,*) 'tra_ldf_bilapg : horizontal biharmonic operator in s-coordinate '79 IF(lwp) WRITE(numout,*) 'tra_ldf_bilapg : horizontal biharmonic operator in s-coordinate on ', cdtype 90 80 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~' 91 81 ENDIF 92 82 93 ! 1. Laplacian of (tb,sb)* aht83 ! 1. Laplacian of ptb * aht 94 84 ! ----------------------------- 95 ! rotated harmonic operator applied to (tb,sb) 96 ! and multiply by aht (output in (wk1,wk2) ) 97 98 CALL ldfght ( kt, tb, sb, wk1, wk2, 1 ) 99 100 101 ! Lateral boundary conditions on (wk1,wk2) (unchanged sign) 102 CALL lbc_lnk( wk1, 'T', 1. ) ; CALL lbc_lnk( wk2, 'T', 1. ) 103 104 ! 2. Bilaplacian of (tb,sb) 85 CALL ldfght( kt, cdtype, ptb, wk1, kjpt, 1 ) ! rotated harmonic operator applied to ptb and multiply by aht 86 ! ! output in wk1 87 ! 88 DO jn = 1, kjpt 89 CALL lbc_lnk( wk1(:,:,:,jn) , 'T', 1. ) ! Lateral boundary conditions on wk1 (unchanged sign) 90 END DO 91 92 ! 2. Bilaplacian of ptb 105 93 ! ------------------------- 106 ! rotated harmonic operator applied to (wk1,wk2) 107 ! (output in (wk3,wk4) ) 108 109 CALL ldfght ( kt, wk1, wk2, wk3, wk4, 2 ) 94 CALL ldfght( kt, cdtype, wk1, wk2, kjpt, 2 ) ! rotated harmonic operator applied to wk1 ; output in wk2 110 95 111 96 112 97 ! 3. Update the tracer trends (j-slab : 2, jpj-1) 113 98 ! --------------------------- 114 ! ! =============== 115 DO jj = 2, jpjm1 ! Vertical slab 116 ! ! =============== 117 DO jk = 1, jpkm1 118 DO ji = 2, jpim1 119 ! add it to the general tracer trends 120 ta(ji,jj,jk) = ta(ji,jj,jk) + wk3(ji,jj,jk) 121 sa(ji,jj,jk) = sa(ji,jj,jk) + wk4(ji,jj,jk) 99 DO jn = 1, kjpt 100 DO jj = 2, jpjm1 101 DO jk = 1, jpkm1 102 DO ji = 2, jpim1 103 ! add it to the general tracer trends 104 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + wk2(ji,jj,jk,jn) 105 END DO 122 106 END DO 123 107 END DO 124 ! ! =============== 125 END DO ! End of slab 126 ! ! =============== 127 108 END DO 109 ! 128 110 END SUBROUTINE tra_ldf_bilapg 129 111 130 112 131 SUBROUTINE ldfght ( kt, pt, ps, plt, pls, kaht )113 SUBROUTINE ldfght ( kt, cdtype, pt, plt, kjpt, kaht ) 132 114 !!---------------------------------------------------------------------- 133 115 !! *** ROUTINE ldfght *** 134 116 !! 135 !! ** Purpose : Apply a geopotential harmonic operator to (pt ,ps) and117 !! ** Purpose : Apply a geopotential harmonic operator to (pt) and 136 118 !! multiply it by the eddy diffusivity coefficient (if kaht=1). 137 119 !! Routine only used in s-coordinates (l_sco=T) with bilaplacian … … 140 122 !! 141 123 !! ** Method : The harmonic operator rotated along geopotential 142 !! surfaces is applied to (pt ,ps) using the slopes of geopotential124 !! surfaces is applied to (pt) using the slopes of geopotential 143 125 !! surfaces computed in inildf routine. The result is provided in 144 126 !! (plt,pls) arrays. It is computed in 2 steps: … … 166 148 !! plt = 1 / (e1t*e2t*e3t) { plt + dk[ zftw ] } 167 149 !! 168 !! * Action : 169 !! 'key_trdtra' defined: the trend is saved for diagnostics. 170 !! 171 !! History : 172 !! 8.0 ! 97-07 (G. Madec) Original code 173 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 174 !!---------------------------------------------------------------------- 175 !! * Arguments 176 INTEGER, INTENT( in ) :: kt ! ocean time-step index 177 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: & 178 pt, ps ! tracer fields (before t and s for 1st call 179 ! ! and laplacian of these fields for 2nd call. 180 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) :: & 181 plt, pls ! partial harmonic operator applied to 182 ! ! pt & ps components except 183 ! ! second order vertical derivative term) 184 INTEGER, INTENT( in ) :: & 185 kaht ! =1 multiply the laplacian by the eddy diffusivity coeff. 186 ! ! =2 no multiplication 187 188 !! * Local declarations 189 INTEGER :: ji, jj, jk ! dummy loop indices 190 REAL(wp) :: & 191 zabe1, zabe2, zmku, zmkv, & ! temporary scalars 192 zbtr, ztah, zsah, ztav, zsav, & 193 zcof0, zcof1, zcof2, & 194 zcof3, zcof4 195 REAL(wp), DIMENSION(jpi,jpj) :: & 196 zftu, zfsu, & ! workspace 197 zdkt, zdk1t, & 198 zdks, zdk1s 199 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 200 zftv, zfsv ! workspace (only v components for ptr) 201 REAL(wp), DIMENSION(jpi,jpk) :: & 202 zftw, zfsw, & ! workspace 203 zdit, zdjt, zdj1t, & 204 zdis, zdjs, zdj1s 205 !!---------------------------------------------------------------------- 206 207 ! ! ********** ! ! =============== 208 DO jk = 1, jpkm1 ! First step ! ! Horizontal slab 209 ! ! ********** ! ! =============== 210 211 ! I.1 Vertical gradient of pt and ps at level jk and jk+1 212 ! ------------------------------------------------------- 213 ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 214 215 zdk1t(:,:) = ( pt(:,:,jk) - pt(:,:,jk+1) ) * tmask(:,:,jk+1) 216 zdk1s(:,:) = ( ps(:,:,jk) - ps(:,:,jk+1) ) * tmask(:,:,jk+1) 217 218 IF( jk == 1 ) THEN 219 zdkt(:,:) = zdk1t(:,:) 220 zdks(:,:) = zdk1s(:,:) 221 ELSE 222 zdkt(:,:) = ( pt(:,:,jk-1) - pt(:,:,jk) ) * tmask(:,:,jk) 223 zdks(:,:) = ( ps(:,:,jk-1) - ps(:,:,jk) ) * tmask(:,:,jk) 150 !!---------------------------------------------------------------------- 151 USE oce , zftv => ua ! use ua as workspace 152 !! 153 INTEGER , INTENT(in ) :: kt ! ocean time-step index 154 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 155 INTEGER , INTENT(in ) :: kjpt !: dimension of 156 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk,kjpt) :: pt ! tracer fields ( before for 1st call 157 ! ! and laplacian of these fields for 2nd call. 158 REAL(wp) , INTENT(out), DIMENSION(jpi,jpj,jpk,kjpt) :: plt !: partial harmonic operator applied to pt components except 159 ! !: second order vertical derivative term 160 INTEGER , INTENT(in ) :: kaht !: =1 multiply the laplacian by the eddy diffusivity coeff. 161 ! !: =2 no multiplication 162 !! 163 INTEGER :: ji, jj, jk,jn ! dummy loop indices 164 ! ! temporary scalars 165 REAL(wp) :: zabe1, zabe2, zmku, zmkv 166 REAL(wp) :: zbtr, ztah, ztav 167 REAL(wp) :: zcof0, zcof1, zcof2, zcof3, zcof4 168 REAL(wp), DIMENSION(jpi,jpj) :: zftu, zdkt, zdk1t ! workspace 169 REAL(wp), DIMENSION(jpi,jpk) :: zftw, zdit, zdjt, zdj1t ! 170 !!---------------------------------------------------------------------- 171 172 ! 173 DO jn = 1, kjpt 174 ! ! ********** ! ! =============== 175 DO jk = 1, jpkm1 ! First step ! ! Horizontal slab 176 ! ! ********** ! ! =============== 177 178 ! I.1 Vertical gradient of pt and ps at level jk and jk+1 179 ! ------------------------------------------------------- 180 ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 181 182 zdk1t(:,:) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1) 183 IF( jk == 1 ) THEN 184 zdkt(:,:) = zdk1t(:,:) 185 ELSE 186 zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk) 187 ENDIF 188 189 190 ! I.2 Horizontal fluxes 191 ! --------------------- 192 193 DO jj = 1, jpjm1 194 DO ji = 1, jpim1 195 zabe1 = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) 196 zabe2 = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) 197 198 zmku = 1./MAX( tmask(ji+1,jj,jk )+tmask(ji,jj,jk+1) & 199 & +tmask(ji+1,jj,jk+1)+tmask(ji,jj,jk ),1. ) 200 zmkv = 1./MAX( tmask(ji,jj+1,jk )+tmask(ji,jj,jk+1) & 201 & +tmask(ji,jj+1,jk+1)+tmask(ji,jj,jk ),1. ) 202 203 zcof1 = -e2u(ji,jj) * uslp(ji,jj,jk) * zmku 204 zcof2 = -e1v(ji,jj) * vslp(ji,jj,jk) * zmkv 205 206 zftu(ji,jj)= umask(ji,jj,jk) * & 207 & ( zabe1 *( pt (ji+1,jj,jk,jn) - pt(ji,jj,jk,jn) ) & 208 & + zcof1 *( zdkt (ji+1,jj) + zdk1t(ji,jj) & 209 & +zdk1t(ji+1,jj) + zdkt (ji,jj) ) ) 210 211 zftv(ji,jj,jk)= vmask(ji,jj,jk) * & 212 & ( zabe2 *( pt(ji,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) & 213 & + zcof2 *( zdkt (ji,jj+1) + zdk1t(ji,jj) & 214 & +zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) 215 END DO 216 END DO 217 218 219 ! I.3 Second derivative (divergence) (not divided by the volume) 220 ! --------------------- 221 222 DO jj = 2 , jpjm1 223 DO ji = 2 , jpim1 224 ztah = zftu(ji,jj) - zftu(ji-1,jj) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) 225 plt(ji,jj,jk,jn) = ztah 226 END DO 227 END DO 228 ! ! =============== 229 END DO ! End of slab 230 ! ! =============== 231 ! "Poleward" diffusive heat or salt transport 232 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( kaht == 2 ) .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 233 IF( jn == jp_tem) htr_ldf(:) = ptr_vj( zftv(:,:,:) ) 234 IF( jn == jp_sal) str_ldf(:) = ptr_vj( zftv(:,:,:) ) 224 235 ENDIF 225 236 226 227 ! I.2 Horizontal fluxes 228 ! --------------------- 229 230 DO jj = 1, jpjm1 231 DO ji = 1, jpim1 232 zabe1 = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) 233 zabe2 = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) 234 235 zmku=1./MAX( tmask(ji+1,jj,jk )+tmask(ji,jj,jk+1) & 236 +tmask(ji+1,jj,jk+1)+tmask(ji,jj,jk ),1. ) 237 zmkv=1./MAX( tmask(ji,jj+1,jk )+tmask(ji,jj,jk+1) & 238 +tmask(ji,jj+1,jk+1)+tmask(ji,jj,jk ),1. ) 239 240 zcof1= -e2u(ji,jj) * uslp(ji,jj,jk) * zmku 241 zcof2= -e1v(ji,jj) * vslp(ji,jj,jk) * zmkv 242 243 zftu(ji,jj)= umask(ji,jj,jk) * & 244 ( zabe1 *( pt(ji+1,jj,jk) - pt(ji,jj,jk) ) & 245 + zcof1 *( zdkt (ji+1,jj) + zdk1t(ji,jj) & 246 +zdk1t(ji+1,jj) + zdkt (ji,jj) ) ) 247 248 zftv(ji,jj,jk)= vmask(ji,jj,jk) * & 249 ( zabe2 *( pt(ji,jj+1,jk) - pt(ji,jj,jk) ) & 250 + zcof2 *( zdkt (ji,jj+1) + zdk1t(ji,jj) & 251 +zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) 252 253 zfsu(ji,jj)= umask(ji,jj,jk) * & 254 ( zabe1 *( ps(ji+1,jj,jk) - ps(ji,jj,jk) ) & 255 + zcof1 *( zdks (ji+1,jj) + zdk1s(ji,jj) & 256 +zdk1s(ji+1,jj) + zdks (ji,jj) ) ) 257 258 zfsv(ji,jj,jk)= vmask(ji,jj,jk) * & 259 ( zabe2 *( ps(ji,jj+1,jk) - ps(ji,jj,jk) ) & 260 + zcof2 *( zdks (ji,jj+1) + zdk1s(ji,jj) & 261 +zdk1s(ji,jj+1) + zdks (ji,jj) ) ) 262 END DO 263 END DO 264 265 266 ! I.3 Second derivative (divergence) (not divided by the volume) 267 ! --------------------- 268 269 DO jj = 2 , jpjm1 270 DO ji = 2 , jpim1 271 ztah = zftu(ji,jj) - zftu(ji-1,jj) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) 272 zsah = zfsu(ji,jj) - zfsu(ji-1,jj) + zfsv(ji,jj,jk) - zfsv(ji,jj-1,jk) 273 plt(ji,jj,jk) = ztah 274 pls(ji,jj,jk) = zsah 275 END DO 276 END DO 277 ! ! =============== 278 END DO ! End of slab 279 ! ! =============== 280 281 !!but this should be done somewhere after 282 ! "zonal" mean diffusive heat and salt transport 283 IF( ln_diaptr .AND. ( kaht == 2 ) .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 284 pht_ldf(:) = ptr_vj( zftv(:,:,:) ) 285 pst_ldf(:) = ptr_vj( zfsv(:,:,:) ) 286 ENDIF 287 288 ! ! ************ ! ! =============== 289 DO jj = 2, jpjm1 ! Second step ! ! Horizontal slab 290 ! ! ************ ! ! =============== 291 292 ! II.1 horizontal tracer gradient 293 ! ------------------------------- 294 295 DO jk = 1, jpk 296 DO ji = 1, jpim1 297 zdit (ji,jk) = ( pt(ji+1,jj ,jk) - pt(ji,jj ,jk) ) * umask(ji,jj ,jk) 298 zdis (ji,jk) = ( ps(ji+1,jj ,jk) - ps(ji,jj ,jk) ) * umask(ji,jj ,jk) 299 zdjt (ji,jk) = ( pt(ji ,jj+1,jk) - pt(ji,jj ,jk) ) * vmask(ji,jj ,jk) 300 zdjs (ji,jk) = ( ps(ji ,jj+1,jk) - ps(ji,jj ,jk) ) * vmask(ji,jj ,jk) 301 zdj1t(ji,jk) = ( pt(ji ,jj ,jk) - pt(ji,jj-1,jk) ) * vmask(ji,jj-1,jk) 302 zdj1s(ji,jk) = ( ps(ji ,jj ,jk) - ps(ji,jj-1,jk) ) * vmask(ji,jj-1,jk) 303 END DO 304 END DO 305 306 307 ! II.2 Vertical fluxes 308 ! -------------------- 309 310 ! Surface and bottom vertical fluxes set to zero 311 zftw(:, 1 ) = 0.e0 312 zfsw(:, 1 ) = 0.e0 313 zftw(:,jpk) = 0.e0 314 zfsw(:,jpk) = 0.e0 315 316 ! interior (2=<jk=<jpk-1) 317 DO jk = 2, jpkm1 318 DO ji = 2, jpim1 319 zcof0 = e1t(ji,jj) * e2t(ji,jj) / fse3w(ji,jj,jk) & 320 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 321 + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) 322 323 zmku =1./MAX( umask(ji ,jj,jk-1)+umask(ji-1,jj,jk) & 324 +umask(ji-1,jj,jk-1)+umask(ji ,jj,jk), 1. ) 325 326 zmkv =1./MAX( vmask(ji,jj ,jk-1)+vmask(ji,jj-1,jk) & 327 +vmask(ji,jj-1,jk-1)+vmask(ji,jj ,jk), 1. ) 328 329 zcof3 = - e2t(ji,jj) * wslpi (ji,jj,jk) * zmku 330 zcof4 = - e1t(ji,jj) * wslpj (ji,jj,jk) * zmkv 331 332 zftw(ji,jk) = tmask(ji,jj,jk) * & 333 ( zcof0 * ( pt (ji,jj,jk-1) - pt (ji,jj,jk) ) & 334 + zcof3 * ( zdit (ji ,jk-1) + zdit (ji-1,jk) & 335 +zdit (ji-1,jk-1) + zdit (ji ,jk) ) & 336 + zcof4 * ( zdjt (ji ,jk-1) + zdj1t(ji ,jk) & 337 +zdj1t(ji ,jk-1) + zdjt (ji ,jk) ) ) 338 339 zfsw(ji,jk) = tmask(ji,jj,jk) * & 340 ( zcof0 * ( ps (ji,jj,jk-1) - ps (ji,jj,jk) ) & 341 + zcof3 * ( zdis (ji ,jk-1) + zdis (ji-1,jk) & 342 +zdis (ji-1,jk-1) + zdis (ji ,jk) ) & 343 + zcof4 * ( zdjs (ji ,jk-1) + zdj1s(ji ,jk) & 344 +zdj1s(ji ,jk-1) + zdjs (ji ,jk) ) ) 345 END DO 346 END DO 347 348 349 ! II.3 Divergence of vertical fluxes added to the horizontal divergence 350 ! --------------------------------------------------------------------- 351 352 IF( kaht == 1 ) THEN 353 ! multiply the laplacian by the eddy diffusivity coefficient 354 DO jk = 1, jpkm1 237 ! ! ************ ! ! =============== 238 DO jj = 2, jpjm1 ! Second step ! ! Horizontal slab 239 ! ! ************ ! ! =============== 240 241 ! II.1 horizontal tracer gradient 242 ! ------------------------------- 243 244 DO jk = 1, jpk 245 DO ji = 1, jpim1 246 zdit (ji,jk) = ( pt(ji+1,jj ,jk,jn) - pt(ji,jj ,jk,jn) ) * umask(ji,jj ,jk) 247 zdjt (ji,jk) = ( pt(ji ,jj+1,jk,jn) - pt(ji,jj ,jk,jn) ) * vmask(ji,jj ,jk) 248 zdj1t(ji,jk) = ( pt(ji ,jj ,jk,jn) - pt(ji,jj-1,jk,jn) ) * vmask(ji,jj-1,jk) 249 END DO 250 END DO 251 252 253 ! II.2 Vertical fluxes 254 ! -------------------- 255 256 ! Surface and bottom vertical fluxes set to zero 257 zftw(:, 1 ) = 0.e0 258 zftw(:,jpk) = 0.e0 259 260 ! interior (2=<jk=<jpk-1) 261 DO jk = 2, jpkm1 355 262 DO ji = 2, jpim1 356 ! eddy coef. divided by the volume element 357 zbtr = fsahtt(ji,jj,jk) / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) ) 358 ! vertical divergence 359 ztav = zftw(ji,jk) - zftw(ji,jk+1) 360 zsav = zfsw(ji,jk) - zfsw(ji,jk+1) 361 ! harmonic operator applied to (pt,ps) and multiply by aht 362 plt(ji,jj,jk) = ( plt(ji,jj,jk) + ztav ) * zbtr 363 pls(ji,jj,jk) = ( pls(ji,jj,jk) + zsav ) * zbtr 364 END DO 365 END DO 366 ELSEIF( kaht == 2 ) THEN 367 ! second call, no multiplication 368 DO jk = 1, jpkm1 369 DO ji = 2, jpim1 370 ! inverse of the volume element 371 zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) ) 372 ! vertical divergence 373 ztav = zftw(ji,jk) - zftw(ji,jk+1) 374 zsav = zfsw(ji,jk) - zfsw(ji,jk+1) 375 ! harmonic operator applied to (pt,ps) 376 plt(ji,jj,jk) = ( plt(ji,jj,jk) + ztav ) * zbtr 377 pls(ji,jj,jk) = ( pls(ji,jj,jk) + zsav ) * zbtr 378 END DO 379 END DO 380 ELSE 381 IF(lwp) WRITE(numout,*) ' ldfght: kaht= 1 or 2, here =', kaht 382 IF(lwp) WRITE(numout,*) ' We stop' 383 STOP 'ldfght' 384 ENDIF 385 ! ! =============== 386 END DO ! End of slab 387 ! ! =============== 263 zcof0 = e1t(ji,jj) * e2t(ji,jj) / fse3w(ji,jj,jk) & 264 & * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 265 & + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) 266 267 zmku = 1./MAX( umask(ji ,jj,jk-1)+umask(ji-1,jj,jk) & 268 & +umask(ji-1,jj,jk-1)+umask(ji ,jj,jk), 1. ) 269 270 zmkv = 1./MAX( vmask(ji,jj ,jk-1)+vmask(ji,jj-1,jk) & 271 & +vmask(ji,jj-1,jk-1)+vmask(ji,jj ,jk), 1. ) 272 273 zcof3 = - e2t(ji,jj) * wslpi (ji,jj,jk) * zmku 274 zcof4 = - e1t(ji,jj) * wslpj (ji,jj,jk) * zmkv 275 276 zftw(ji,jk) = tmask(ji,jj,jk) * & 277 & ( zcof0 * ( pt (ji,jj,jk-1,jn) - pt (ji ,jj,jk,jn) ) & 278 & + zcof3 * ( zdit (ji ,jk-1 ) + zdit (ji-1,jk ) & 279 & +zdit (ji-1 ,jk-1 ) + zdit (ji ,jk ) ) & 280 & + zcof4 * ( zdjt (ji ,jk-1 ) + zdj1t(ji ,jk) & 281 & +zdj1t(ji ,jk-1 ) + zdjt (ji ,jk ) ) ) 282 END DO 283 END DO 284 285 286 ! II.3 Divergence of vertical fluxes added to the horizontal divergence 287 ! --------------------------------------------------------------------- 288 289 IF( kaht == 1 ) THEN 290 ! multiply the laplacian by the eddy diffusivity coefficient 291 DO jk = 1, jpkm1 292 DO ji = 2, jpim1 293 ! eddy coef. divided by the volume element 294 zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 295 ! vertical divergence 296 ztav = fsahtt(ji,jj,jk) * ( zftw(ji,jk) - zftw(ji,jk+1) ) 297 ! harmonic operator applied to (pt,ps) and multiply by aht 298 plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr 299 END DO 300 END DO 301 ELSEIF( kaht == 2 ) THEN 302 ! second call, no multiplication 303 DO jk = 1, jpkm1 304 DO ji = 2, jpim1 305 ! inverse of the volume element 306 zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 307 ! vertical divergence 308 ztav = zftw(ji,jk) - zftw(ji,jk+1) 309 ! harmonic operator applied to (pt,ps) 310 plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr 311 END DO 312 END DO 313 ELSE 314 IF(lwp) WRITE(numout,*) ' ldfght: kaht= 1 or 2, here =', kaht 315 IF(lwp) WRITE(numout,*) ' We stop' 316 STOP 'ldfght' 317 ENDIF 318 ! ! =============== 319 END DO ! End of slab 320 ! ! =============== 321 END DO 322 ! 388 323 END SUBROUTINE ldfght 389 324 … … 393 328 !!---------------------------------------------------------------------- 394 329 CONTAINS 395 SUBROUTINE tra_ldf_bilapg( kt ) ! Dummy routine 396 WRITE(*,*) 'tra_ldf_bilapg: You should not have seen this print! error?', kt 330 SUBROUTINE tra_ldf_bilapg( kt, cdtype, ptb, pta, kjpt ) ! Empty routine 331 CHARACTER(len=3) :: cdtype 332 REAL, DIMENSION(:,:,:,:) :: ptb, pta 333 WRITE(*,*) 'tra_ldf_iso: You should not have seen this print! error?', kt, cdtype, ptb(1,1,1,1), pta(1,1,1,1), kjpt 397 334 END SUBROUTINE tra_ldf_bilapg 398 335 #endif
Note: See TracChangeset
for help on using the changeset viewer.
