- Timestamp:
- 2010-10-04T15:53:42+02:00 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA/tranxt.F90
r2120 r2148 15 15 !! 3.0 ! 2008-06 (G. Madec) time stepping always done in trazdf 16 16 !! 3.1 ! 2009-02 (G. Madec, R. Benshila) re-introduce the vvl option 17 !! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA 17 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) semi-implicit hpg with asselin filter + modified LF-RA 18 !! - ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA 18 19 !!---------------------------------------------------------------------- 19 20 20 21 !!---------------------------------------------------------------------- 21 !! tra_nxt : time stepping on t emperature and salinity22 !! tra_nxt_fix : time stepping on t emperature and salinity: fixed volume case23 !! tra_nxt_vvl : time stepping on t emperature and salinity: variable volume case22 !! tra_nxt : time stepping on tracers 23 !! tra_nxt_fix : time stepping on tracers : fixed volume case 24 !! tra_nxt_vvl : time stepping on tracers : variable volume case 24 25 !!---------------------------------------------------------------------- 25 26 USE oce ! ocean dynamics and tracers variables 26 27 USE dom_oce ! ocean space and time domain variables 28 USE sbc_oce ! surface boundary condition: ocean 27 29 USE zdf_oce ! ??? 28 30 USE domvvl ! variable volume … … 38 40 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 39 41 USE prtctl ! Print control 42 USE traqsr ! penetrative solar radiation (needed for nksr) 40 43 USE traswp ! swap array 41 44 USE agrif_opa_update … … 49 52 PUBLIC tra_nxt_vvl ! to be used in trcnxt 50 53 54 REAL(wp) :: rbcp ! Brown & Campana parameters for semi-implicit hpg 51 55 REAL(wp), DIMENSION(jpk) :: r2dt ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler) 52 56 … … 96 100 IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap' 97 101 IF(lwp) WRITE(numout,*) '~~~~~~~' 102 ! 103 rbcp = 0.25 * (1. + atfp) * (1. + atfp) * ( 1. - atfp) ! Brown & Campana parameter for semi-implicit hpg 98 104 ENDIF 99 105 … … 107 113 #endif 108 114 109 #if defined key_obc 115 #if defined key_obc 110 116 IF( lk_obc ) CALL obc_tra( kt ) ! OBC open boundaries 111 117 #endif … … 114 120 #endif 115 121 #if defined key_agrif 116 CALL Agrif_tra 122 CALL Agrif_tra ! AGRIF zoom boundaries 117 123 #endif 118 124 … … 133 139 134 140 ! Leap-Frog + Asselin filter time stepping 135 IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt, tsb, tsn, tsa, jpts ) ! variable volume level (vvl)136 ELSE ; CALL tra_nxt_fix( kt, tsb, tsn, tsa, jpts ) ! fixed volume level141 IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt, 'TRA', tsb, tsn, tsa, jpts ) ! variable volume level (vvl) 142 ELSE ; CALL tra_nxt_fix( kt, tsb, tsn, tsa, jpts ) ! fixed volume level 137 143 ENDIF 138 144 … … 176 182 !! - swap tracer fields to prepare the next time_step. 177 183 !! This can be summurized for tempearture as: 178 !! ztm = (ta+2tn+tb)/4 ln_dynhpg_imp = T 179 !! ztm = 0 otherwise 184 !! ztm = tn + rbcp * [ta -2 tn + tb ] ln_dynhpg_imp = T 185 !! ztm = 0 otherwise 186 !! with rbcp=1/4 * (1-atfp^4) / (1-atfp) 180 187 !! tb = tn + atfp*[ tb - 2 tn + ta ] 181 !! tn = ta 188 !! tn = ta 182 189 !! ta = ztm (NB: reset to 0 after eos_bn2 call) 183 190 !! … … 185 192 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 186 193 !!---------------------------------------------------------------------- 187 INTEGER , INTENT(in ) :: kt 188 INTEGER , INTENT(in ) :: kjpt 189 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields190 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields191 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta 194 INTEGER , INTENT(in ) :: kt ! ocean time-step index 195 INTEGER , INTENT(in ) :: kjpt ! number of tracers 196 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields 197 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields 198 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend 192 199 !! 193 200 INTEGER :: ji, jj, jk, jn ! dummy loop indices 194 REAL(wp) :: zt m, ztf! temporary scalars201 REAL(wp) :: ztd, ztm ! temporary scalars 195 202 !!---------------------------------------------------------------------- 196 203 … … 205 212 ! ! (only swap) 206 213 DO jn = 1, kjpt 207 DO jk = 1, jpkm1 214 DO jk = 1, jpkm1 208 215 ptn(:,:,jk,jn) = pta(:,:,jk,jn) ! ptb <-- ptn 209 216 END DO … … 219 226 DO jj = 1, jpj 220 227 DO ji = 1, jpi 221 ztm = 0.25 * ( pta(ji,jj,jk,jn) + 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) ) ! mean pt 222 ztf = atfp * ( pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) ) ! Asselin filter on pt 223 ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + ztf ! ptb <-- filtered ptn 224 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 225 pta(ji,jj,jk,jn) = ztm ! pta <-- mean pt 228 ztd = pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) ! time laplacian on tracers 229 ztm = ptn(ji,jj,jk,jn) + rbcp * ztd ! used for both Asselin and Brown & Campana filters 230 ! 231 ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + atfp * ztd ! ptb <-- filtered ptn 232 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 233 pta(ji,jj,jk,jn) = ztm ! pta <-- Brown & Campana average 226 234 END DO 227 235 END DO … … 235 243 DO jj = 1, jpj 236 244 DO ji = 1, jpi 237 ztf = atfp * ( pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) ) ! Asselin filter on t 238 ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + ztf ! ptb <-- filtered ptn 239 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 245 ztd = pta(ji,jj,jk,jn) - 2.* ptn(ji,jj,jk,jn) + ptb(ji,jj,jk,jn) ! time laplacian on tracers 246 ! 247 ptb(ji,jj,jk,jn) = ptn(ji,jj,jk,jn) + atfp * ztd ! ptb <-- filtered ptn 248 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 240 249 END DO 241 250 END DO 242 251 END DO 243 252 END DO 244 !245 253 ENDIF 246 254 ! … … 250 258 251 259 252 SUBROUTINE tra_nxt_vvl( kt, ptb, ptn, pta, kjpt )260 SUBROUTINE tra_nxt_vvl( kt, cdtype, ptb, ptn, pta, kjpt ) 253 261 !!---------------------------------------------------------------------- 254 262 !! *** ROUTINE tra_nxt_vvl *** … … 263 271 !! - swap tracer fields to prepare the next time_step. 264 272 !! This can be summurized for tempearture as: 265 !! ztm = ( e3t_a*ta+2*e3t_n*tn+e3t_b*tb) ln_dynhpg_imp = T266 !! /( e3t_a +2*e3t_n +e3t_b )267 !! ztm = 0 otherwise273 !! ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) ln_dynhpg_imp = T 274 !! /( e3t_n + rbcp*[ e3t_b - 2 e3t_n + e3t_a ] ) 275 !! ztm = 0 otherwise 268 276 !! tb = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) 269 277 !! /( e3t_n + atfp*[ e3t_b - 2 e3t_n + e3t_a ] ) … … 274 282 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 275 283 !!---------------------------------------------------------------------- 276 INTEGER , INTENT(in ) :: kt ! ocean time-step index 277 INTEGER , INTENT(in ) :: kjpt ! number of tracers 278 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields 279 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields 280 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend 284 INTEGER , INTENT(in ) :: kt ! ocean time-step index 285 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 286 INTEGER , INTENT(in ) :: kjpt ! number of tracers 287 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields 288 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields 289 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend 281 290 !! 282 INTEGER :: ji, jj, jk, jn ! dummy loop indices 283 REAL(wp) :: ztm , ztc_f , ztf , ztca, ztcn, ztcb ! temporary scalar 284 REAL(wp) :: ze3mr, ze3fr ! - - 285 REAL(wp) :: ze3t_b, ze3t_n, ze3t_a, ze3t_f ! - - 291 INTEGER :: ji, jj, jk, jn ! dummy loop indices 292 REAL(wp) :: ztc_a , ztc_n , ztc_b ! temporary scalar 293 REAL(wp) :: ztc_f , ztc_d , ztc_m ! - - 294 REAL(wp) :: ze3t_b, ze3t_n, ze3t_a ! - - 295 REAL(wp) :: ze3t_f, ze3t_d, ze3t_m ! - - 296 REAL :: zfact1, zfact2 ! - - 286 297 !!---------------------------------------------------------------------- 287 298 … … 293 304 ! 294 305 ! 295 IF( neuler == 0 .AND. kt == nit000 ) THEN 296 ! 306 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 307 ! ! (only swap) 297 308 DO jn = 1, kjpt 298 DO jk = 1, jpkm1 299 ptn(:,:,jk,jn) = pta(:,:,jk,jn) 309 DO jk = 1, jpkm1 310 ptn(:,:,jk,jn) = pta(:,:,jk,jn) ! ptb <-- ptn 300 311 END DO 301 312 END DO 302 313 ! 303 ELSE ! general case (Leapfrog + Asselin filter314 ELSE ! general case (Leapfrog + Asselin filter) 304 315 ! 305 ! 306 IF( ln_dynhpg_imp ) THEN 307 ! 308 DO jn = 1, kjpt 316 ! ! ----------------------- ! 317 IF( ln_dynhpg_imp ) THEN ! semi-implicite hpg case ! 318 ! ! ----------------------- ! 319 DO jn = 1, kjpt 309 320 DO jk = 1, jpkm1 310 321 DO jj = 1, jpj … … 314 325 ze3t_a = fse3t_a(ji,jj,jk) 315 326 ! ! tracer content at Before, now and after 316 ztcb = ptb(ji,jj,jk,jn) * ze3t_b 317 ztcn = ptn(ji,jj,jk,jn) * ze3t_n 318 ztca = pta(ji,jj,jk,jn) * ze3t_a 319 ! 320 ! ! Asselin filter on thickness and tracer content 321 ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b ) 322 ztc_f = atfp * ( ztca - 2.* ztcn + ztcb ) 323 ! 324 ! ! filtered tracer including the correction 325 ze3fr = 1.e0 / ( ze3t_n + ze3t_f ) 326 ztf = ze3fr * ( ztcn + ztc_f ) 327 ! ! mean thickness and tracer 328 ze3mr = 1.e0 / ( ze3t_a + 2.* ze3t_n + ze3t_b ) 329 ztm = ze3mr * ( ztca + 2.* ztcn + ztcb ) 330 !!gm mean e3t have to be saved and used in dynhpg or it can be recomputed in dynhpg !! 331 !!gm e3t_m(ji,jj,jk) = 0.25 / ze3mr 327 ztc_b = ptb(ji,jj,jk,jn) * ze3t_b 328 ztc_n = ptn(ji,jj,jk,jn) * ze3t_n 329 ztc_a = pta(ji,jj,jk,jn) * ze3t_a 330 ! ! Time laplacian on tracer contents 331 ! ! used for both Asselin and Brown & Campana filters 332 ze3t_d = ze3t_b - 2.* ze3t_n + ze3t_a 333 ztc_d = ztc_b - 2.* ztc_n + ztc_a 334 ! ! Asselin Filter on thicknesses and tracer contents 335 ztc_f = ztc_n + atfp * ztc_d 336 ztc_m = ztc_n + rbcp * ztc_d 337 ! 338 ze3t_f = 1.0 / ( ze3t_n + atfp * ze3t_d ) 339 ze3t_m = 1.0 / ( ze3t_n + rbcp * ze3t_d ) 332 340 ! ! swap of arrays 333 ptb(ji,jj,jk,jn) = zt f ! ptb <-- ptn + filter334 pt n(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta335 pt a(ji,jj,jk,jn) = ztm ! pta <-- mean t341 ptb(ji,jj,jk,jn) = ztc_f * ze3t_f ! ptb <-- ptn filtered 342 pta(ji,jj,jk,jn) = ztc_m * ze3t_m ! pta <-- Brown & Campana average 343 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 336 344 END DO 337 345 END DO 338 346 END DO 339 347 END DO 340 ! ! ----------------------- ! 341 ELSE ! explicit hpg case ! 342 ! ! ----------------------- ! 343 DO jn = 1, kjpt 344 DO jk = 1, jpkm1 345 DO jj = 1, jpj 346 DO ji = 1, jpi 347 ze3t_b = fse3t_b(ji,jj,jk) 348 ze3t_n = fse3t_n(ji,jj,jk) 349 ze3t_a = fse3t_a(ji,jj,jk) 350 ! ! tracer content at Before, now and after 351 ztcb = ptb(ji,jj,jk,jn) * ze3t_b 352 ztcn = ptn(ji,jj,jk,jn) * ze3t_n 353 ztca = pta(ji,jj,jk,jn) * ze3t_a 354 ! 355 ! ! Asselin filter on thickness and tracer content 356 ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b ) 357 ztc_f = atfp * ( ztca - 2.* ztcn + ztcb ) 358 ! 359 ! ! filtered tracer including the correction 360 ze3fr = 1.e0 / ( ze3t_n + ze3t_f ) 361 ztf = ( ztcn + ztc_f ) * ze3fr 362 ! ! swap of arrays 363 ptb(ji,jj,jk,jn) = ztf ! tb <-- tn filtered 364 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! tn <-- ta 348 ! ! ----------------------- ! 349 ELSE ! explicit hpg case ! 350 ! ! ----------------------- ! 351 IF( cdtype == 'TRA' ) THEN 352 ! 353 DO jn = 1, kjpt 354 DO jk = 1, jpkm1 355 zfact1 = atfp * rdttra(jk) 356 zfact2 = zfact1 / rau0 357 DO jj = 1, jpj 358 DO ji = 1, jpi 359 ze3t_b = fse3t_b(ji,jj,jk) 360 ze3t_n = fse3t_n(ji,jj,jk) 361 ze3t_a = fse3t_a(ji,jj,jk) 362 ! ! tracer content at Before, now and after 363 ztc_b = ptb(ji,jj,jk,jn) * ze3t_b 364 ztc_n = ptn(ji,jj,jk,jn) * ze3t_n 365 ztc_a = pta(ji,jj,jk,jn) * ze3t_a 366 ! 367 ze3t_f = ze3t_n + atfp * ( ze3t_b - 2. * ze3t_n + ze3t_a ) 368 ztc_f = ztc_n + atfp * ( ztc_a - 2. * ztc_n + ztc_b ) 369 370 IF( jk == 1 ) THEN 371 ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) ) 372 IF( jn == jp_tem ) ztc_f = ztc_f - zfact1 * ( sbc_hc_n(ji,jj) - sbc_hc_b(ji,jj) ) 373 IF( jn == jp_sal ) ztc_f = ztc_f - zfact1 * ( sbc_sc_n(ji,jj) - sbc_sc_b(ji,jj) ) 374 ENDIF 375 IF( jn == jp_tem .AND. ln_traqsr .AND. jk <= nksr ) & 376 & ztc_f = ztc_f - zfact1 * ( qsr_hc_n(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 377 378 ze3t_f = 1.e0 / ze3t_f 379 ptb(ji,jj,jk,jn) = ztc_f * ze3t_f ! tb <-- tn filtered 380 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! tn <-- ta 381 END DO 365 382 END DO 366 383 END DO 367 384 END DO 368 END DO 385 ! 386 ELSE IF( cdtype == 'TRC' ) THEN 387 ! 388 DO jn = 1, kjpt 389 DO jk = 1, jpkm1 390 DO jj = 1, jpj 391 DO ji = 1, jpi 392 ze3t_b = fse3t_b(ji,jj,jk) 393 ze3t_n = fse3t_n(ji,jj,jk) 394 ze3t_a = fse3t_a(ji,jj,jk) 395 ! ! tracer content at Before, now and after 396 ztc_b = ptb(ji,jj,jk,jn) * ze3t_b 397 ztc_n = ptn(ji,jj,jk,jn) * ze3t_n 398 ztc_a = pta(ji,jj,jk,jn) * ze3t_a 399 ! 400 ze3t_f = ze3t_n + atfp * ( ze3t_b - 2. * ze3t_n + ze3t_a ) 401 ztc_f = ztc_n + atfp * ( ztc_a - 2. * ztc_n + ztc_b ) 402 403 ze3t_f = 1.e0 / ze3t_f 404 ptb(ji,jj,jk,jn) = ztc_f * ze3t_f ! tb <-- tn filtered 405 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! tn <-- ta 406 END DO 407 END DO 408 END DO 409 END DO 410 ! 411 END IF 369 412 ! 370 413 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
