Changeset 1975 for branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90
- Timestamp:
- 2010-06-28T19:22:14+02:00 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90
r1870 r1975 25 25 USE oce ! ocean dynamics and tracers variables 26 26 USE dom_oce ! ocean space and time domain variables 27 USE sbc_oce ! surface boundary condition: ocean 27 28 USE zdf_oce ! ??? 28 29 USE domvvl ! variable volume … … 37 38 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 38 39 USE prtctl ! Print control 40 USE traqsr ! penetrative solar radiation (needed for nksr) 39 41 USE agrif_opa_update 40 42 USE agrif_opa_interp … … 45 47 PUBLIC tra_nxt ! routine called by step.F90 46 48 47 REAL(wp) :: rbcp , r1_2bcp! Brown & Campana parameters for semi-implicit hpg49 REAL(wp) :: rbcp ! Brown & Campana parameters for semi-implicit hpg 48 50 REAL(wp), DIMENSION(jpk) :: r2dt_t ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler) 49 51 … … 98 100 IF(lwp) WRITE(numout,*) '~~~~~~~' 99 101 ! 100 rbcp = 0.25 * (1 + atfp) * (1 + atfp * atfp) ! Brown & Campana parameter for semi-implicit hpg 101 r1_2bcp = 1.e0 - 2.e0 * rbcp 102 rbcp = 0.25 * (1. + atfp) * (1. + atfp) * ( 1. - atfp) ! Brown & Campana parameter for semi-implicit hpg 102 103 ENDIF 103 104 ! ! set time step size (Euler/Leapfrog) … … 164 165 !! - swap tracer fields to prepare the next time_step. 165 166 !! This can be summurized for temperature as: 166 !! ztm = (rbcp*ta+(2-rbcp)*tn+rbcp*tb)ln_dynhpg_imp = T167 !! ztm = tn + rbcp * [ta -2 tn + tb ] ln_dynhpg_imp = T 167 168 !! ztm = 0 otherwise 168 !! with rbcp= (1+atfp)*(1+atfp*atfp)/4169 !! with rbcp=1/4 * (1-atfp^4) / (1-atfp) 169 170 !! tb = tn + atfp*[ tb - 2 tn + ta ] 170 171 !! tn = ta … … 177 178 !! 178 179 INTEGER :: ji, jj, jk ! dummy loop indices 179 REAL(wp) :: zt m, ztf, zac! temporary scalars180 REAL(wp) :: z sm, zsf! - -180 REAL(wp) :: zt_m, zs_m ! temporary scalars 181 REAL(wp) :: ztn, zsn ! - - 181 182 !!---------------------------------------------------------------------- 182 183 … … 187 188 ENDIF 188 189 ! 189 ! ! ----------------------- ! 190 IF( ln_dynhpg_imp ) THEN ! semi-implicite hpg case ! 191 ! ! ----------------------- ! 192 ! 193 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 194 DO jk = 1, jpkm1 ! (only swap) 195 tn(:,:,jk) = ta(:,:,jk) ! tb <-- tn 196 sn(:,:,jk) = sa(:,:,jk) 197 END DO 198 ELSE ! general case (Leapfrog + Asselin filter) 199 DO jk = 1, jpkm1 200 DO jj = 1, jpj 201 DO ji = 1, jpi 202 ztm = rbcp * ( ta(ji,jj,jk) + tb(ji,jj,jk) ) + r1_2bcp * tn(ji,jj,jk) ! mean t 203 zsm = rbcp * ( sa(ji,jj,jk) + sb(ji,jj,jk) ) + r1_2bcp * sn(ji,jj,jk) 204 ztf = atfp * ( ta(ji,jj,jk) + tb(ji,jj,jk) - 2. * tn(ji,jj,jk) ) ! Asselin filter on t 205 zsf = atfp * ( sa(ji,jj,jk) + sb(ji,jj,jk) - 2. * sn(ji,jj,jk) ) 206 tb(ji,jj,jk) = tn(ji,jj,jk) + ztf ! tb <-- filtered tn 207 sb(ji,jj,jk) = sn(ji,jj,jk) + zsf 208 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 209 sn(ji,jj,jk) = sa(ji,jj,jk) 210 ta(ji,jj,jk) = ztm ! ta <-- mean t 211 sa(ji,jj,jk) = zsm 212 END DO 190 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 191 DO jk = 1, jpkm1 ! (only swap) 192 tn(:,:,jk) = ta(:,:,jk) ! tn <-- ta 193 sn(:,:,jk) = sa(:,:,jk) ! sn <-- sa 194 END DO 195 ELSE ! General case (Leapfrog + Asselin filter) 196 DO jk = 1, jpkm1 197 DO jj = 1, jpj 198 DO ji = 1, jpi 199 IF( ln_dynhpg_imp ) THEN ! implicit hpg: keep tn, sn in memory 200 ztn = tn(ji,jj,jk) 201 zsn = sn(ji,jj,jk) 202 ENDIF 203 ! ! time laplacian on tracers 204 ! ! used for both Asselin and Brown & Campana filters 205 zt_m = ta(ji,jj,jk) - 2. * tn(ji,jj,jk) + tb(ji,jj,jk) 206 zs_m = sa(ji,jj,jk) - 2. * sn(ji,jj,jk) + sb(ji,jj,jk) 207 ! 208 ! ! swap of arrays 209 tb(ji,jj,jk) = tn(ji,jj,jk) + atfp * zt_m ! tb <-- tn filtered 210 sb(ji,jj,jk) = sn(ji,jj,jk) + atfp * zs_m ! sb <-- sn filtered 211 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 212 sn(ji,jj,jk) = sa(ji,jj,jk) ! sn <-- sa 213 ! ! semi imlicit hpg computation (Brown & Campana) 214 IF( ln_dynhpg_imp ) THEN 215 ta(ji,jj,jk) = ztn + rbcp * zt_m ! ta <-- Brown & Campana average 216 sa(ji,jj,jk) = zsn + rbcp * zs_m ! sa <-- Brown & Campana average 217 ENDIF 213 218 END DO 214 219 END DO 215 ENDIF 216 ! ! ----------------------- ! 217 ELSE ! explicit hpg case ! 218 ! ! ----------------------- ! 219 ! 220 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 221 DO jk = 1, jpkm1 222 tn(:,:,jk) = ta(:,:,jk) ! tn <-- ta 223 sn(:,:,jk) = sa(:,:,jk) 224 END DO 225 ELSE ! general case (Leapfrog + Asselin filter) 226 DO jk = 1, jpkm1 227 DO jj = 1, jpj 228 DO ji = 1, jpi 229 ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) ) ! Asselin filter on t 230 zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 231 tb(ji,jj,jk) = tn(ji,jj,jk) + ztf ! tb <-- filtered tn 232 sb(ji,jj,jk) = sn(ji,jj,jk) + zsf 233 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 234 sn(ji,jj,jk) = sa(ji,jj,jk) 235 END DO 236 END DO 237 END DO 238 ENDIF 239 ENDIF 240 ! 241 !!gm from Matthieu : unchecked 242 IF( neuler /= 0 .OR. kt /= nit000 ) THEN ! remove the forcing from the Asselin filter 243 zac = atfp * rdttra(1) 244 tb(:,:,1) = tb(:,:,1) - zac * ( qns (:,:) - qns_b (:,:) ) ! non solar surface flux 245 sb(:,:,1) = sn(:,:,1) - zac * ( emps(:,:) - emps_b(:,:) ) ! surface salt flux 246 ! 247 IF( ln_traqsr ) ! solar penetrating flux 248 DO jk = 1, nksr 249 DO jj = 1, jpj 250 DO ji = 1, jpi 251 IF( ln_sco ) THEN 252 z_cor_qsr = etot3(ji,jj,jk) * ( qsr(ji,jj) - qsrb(ji,jj) ) 253 ELSEIF( ln_zps .OR. ln_zco ) THEN 254 z_cor_qsr = ( qsr(ji,jj) - qsrb(ji,jj) ) * & 255 & ( gdsr(jk) * tmask(ji,jj,jk) - gdsr(jk+1) * tmask(ji,jj,jk+1) ) 256 ENDIF 257 zt = zt - zfact1 * z_cor_qsr 258 END DO 259 END DO 220 END DO 260 221 ENDIF 261 222 ! … … 276 237 !! - swap tracer fields to prepare the next time_step. 277 238 !! This can be summurized for temperature as: 278 !! ztm = ( rbcp*e3t_a*ta + (2-rbtp)*e3t_n*tn + rbcp*e3t_b*tb) ln_dynhpg_imp = T279 !! /( rbcp*e3t_a + (2-rbcp)*e3t_n + rbcp*e3t_b)239 !! ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) ln_dynhpg_imp = T 240 !! /( e3t_n + rbcp*[ e3t_b - 2 e3t_n + e3t_a ] ) 280 241 !! ztm = 0 otherwise 281 242 !! tb = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) … … 287 248 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 288 249 !!---------------------------------------------------------------------- 289 INTEGER, INTENT(in) :: kt ! ocean time-step index250 INTEGER, INTENT(in) :: kt ! ocean time-step index 290 251 !! 291 INTEGER :: ji, jj, jk ! dummy loop indices 292 REAL(wp) :: ztm , ztc_f , ztf , ztca, ztcn, ztcb ! temporary scalar 293 REAL(wp) :: zsm , zsc_f , zsf , zsca, zscn, zscb ! - - 294 REAL(wp) :: ze3mr, ze3fr ! - - 295 REAL(wp) :: ze3t_b, ze3t_n, ze3t_a, ze3t_f ! - - 296 !!---------------------------------------------------------------------- 297 252 INTEGER :: ji, jj, jk ! dummy loop indices 253 REAL :: ze3t_a, ze3t_n, ze3t_b ! temporary scalars 254 REAL :: ztc_a, ztc_n, ztc_b ! - - 255 REAL :: zsc_a, zsc_n, zsc_b ! - - 256 REAL :: ztc_f, zsc_f, ztc_m, zsc_m ! - - 257 REAL :: ze3t_f, ze3t_m ! - - 258 REAL :: zfact1, zfact2 ! - - 259 !!---------------------------------------------------------------------- 260 298 261 IF( kt == nit000 ) THEN 299 262 IF(lwp) WRITE(numout,*) … … 302 265 ENDIF 303 266 304 ! ! ----------------------- ! 305 IF( ln_dynhpg_imp ) THEN ! semi-implicite hpg case ! (mean tracer computation) 306 ! ! ----------------------- ! 307 ! 308 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 309 DO jk = 1, jpkm1 ! (only swap) 310 tn(:,:,jk) = ta(:,:,jk) ! tn <-- ta 311 sn(:,:,jk) = sa(:,:,jk) 312 END DO 313 ! ! general case (Leapfrog + Asselin filter) 314 ELSE ! apply filter on thickness weighted tracer and swap 315 DO jk = 1, jpkm1 316 DO jj = 1, jpj 317 DO ji = 1, jpi 318 ze3t_b = fse3t_b(ji,jj,jk) 319 ze3t_n = fse3t_n(ji,jj,jk) 320 ze3t_a = fse3t_a(ji,jj,jk) 321 ! ! tracer content at Before, now and after 322 ztcb = tb(ji,jj,jk) * ze3t_b ; zscb = sb(ji,jj,jk) * ze3t_b 323 ztcn = tn(ji,jj,jk) * ze3t_n ; zscn = sn(ji,jj,jk) * ze3t_n 324 ztca = ta(ji,jj,jk) * ze3t_a ; zsca = sa(ji,jj,jk) * ze3t_a 325 ! 326 ! ! Asselin filter on thickness and tracer content 327 ze3t_f = atfp * ( ze3t_a + ze3t_b - 2.* ze3t_n ) 328 ztc_f = atfp * ( ztca + ztcb - 2.* ztcn ) 329 zsc_f = atfp * ( zsca + zscb - 2.* zscn ) 330 ! 331 ! ! filtered tracer including the correction 332 ze3fr = 1.e0 / ( ze3t_n + ze3t_f ) 333 ztf = ze3fr * ( ztcn + ztc_f ) 334 zsf = ze3fr * ( zscn + zsc_f ) 335 ! ! mean thickness and tracer 336 ze3mr = 1.e0 / ( rbcp * ( ze3t_a + ze3t_b ) + r1_2bcp * ze3t_n ) 337 ztm = ze3mr * ( rbcp * ( ztca + ztcb ) + r1_2bcp * ztcn ) 338 zsm = ze3mr * ( rbcp * ( zsca + zscb ) + r1_2bcp * zscn ) 339 !!gm mean e3t have to be saved and used in dynhpg or it can be recomputed in dynhpg !! 340 !!gm e3t_m(ji,jj,jk) = 1.e0 / ze3mr 341 ! ! swap of arrays 342 tb(ji,jj,jk) = ztf ! tb <-- tn + filter 343 sb(ji,jj,jk) = zsf 344 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 345 sn(ji,jj,jk) = sa(ji,jj,jk) 346 ta(ji,jj,jk) = ztm ! ta <-- mean t 347 sa(ji,jj,jk) = zsm 348 END DO 267 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 268 ! ! (only swap) 269 DO jk = 1, jpkm1 ! tn <-- ta 270 tn(:,:,jk) = ta(:,:,jk) ! sn <-- sa 271 sn(:,:,jk) = sa(:,:,jk) 272 END DO 273 ! ! General case (Leapfrog + Asselin filter) 274 ELSE ! apply filter on thickness weighted tracer and swap 275 DO jk = 1, jpkm1 276 zfact1 = atfp * r2dt_t(jk) 277 zfact2 = zfact1 / rau0 278 DO jj = 1, jpj 279 DO ji = 1, jpi 280 ! ! scale factors at Before, now and after 281 ze3t_b = fse3t_b(ji,jj,jk) 282 ze3t_n = fse3t_n(ji,jj,jk) 283 ze3t_a = fse3t_a(ji,jj,jk) 284 ! ! tracer content at Before, now and after 285 ztc_b = tb(ji,jj,jk) * ze3t_b ; zsc_b = sb(ji,jj,jk) * ze3t_b 286 ztc_n = tn(ji,jj,jk) * ze3t_n ; zsc_n = sn(ji,jj,jk) * ze3t_n 287 ztc_a = ta(ji,jj,jk) * ze3t_a ; zsc_a = sa(ji,jj,jk) * ze3t_a 288 ! 289 ! ! Time laplacian on thickness and tracer content 290 ! ! used for both Asselin and Brown & Campana filters 291 ze3t_m = ze3t_a - 2. * ze3t_n + ze3t_b 292 ztc_m = ztc_a - 2. * ztc_n + ztc_b 293 zsc_m = zsc_a - 2. * zsc_n + zsc_b 294 ! ! Asselin Filter + correction 295 ze3t_f = ze3t_n + atfp * ze3t_m 296 ztc_f = ztc_n + atfp * ztc_m 297 zsc_f = zsc_n + atfp * zsc_m 298 299 IF( jk == 1 ) THEN 300 ze3t_f = ze3t_f - zfact2 * ( emp_b (ji,jj) - emp (ji,jj) ) 301 ztc_f = ztc_f - zfact1 * ( sbc_trd_hc_n(ji,jj) - sbc_trd_hc_b(ji,jj) ) 302 ENDIF 303 IF( ln_traqsr .AND. ( jk .LE. nksr ) ) THEN 304 ztc_f = ztc_f - zfact1 * ( qsr_trd_hc_n(ji,jj,jk) - qsr_trd_hc_b(ji,jj,jk) ) 305 ENDIF 306 ! ! swap of arrays 307 ze3t_f = 1.e0 / ze3t_f 308 tb(ji,jj,jk) = ztc_f * ze3t_f ! tb <-- tn filtered 309 sb(ji,jj,jk) = zsc_f * ze3t_f ! sb <-- sn filtered 310 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 311 sn(ji,jj,jk) = sa(ji,jj,jk) ! sn <-- sa 312 ! ! semi imlicit hpg computation (Brown & Campana) 313 IF( ln_dynhpg_imp ) THEN 314 ze3t_m = 1.e0 / ( ze3t_n + rbcp * ze3t_m ) 315 ta(ji,jj,jk) = ( ztc_n + rbcp * ztc_m ) * ze3t_m ! ta <-- Brown & Campana average 316 sa(ji,jj,jk) = ( zsc_n + rbcp * zsc_m ) * ze3t_m ! sa <-- Brown & Campana average 317 ENDIF 349 318 END DO 350 319 END DO 351 ENDIF 352 ! ! ----------------------- ! 353 ELSE ! explicit hpg case ! (NO mean tracer computation) 354 ! ! ----------------------- ! 355 ! 356 IF( neuler == 0 .AND. kt == nit000 ) THEN ! case of Euler time-stepping at first time-step 357 DO jk = 1, jpkm1 ! (only swap) 358 tn(:,:,jk) = ta(:,:,jk) ! tn <-- ta 359 sn(:,:,jk) = sa(:,:,jk) 360 END DO 361 ! ! general case (Leapfrog + Asselin filter) 362 ELSE ! apply filter on thickness weighted tracer and swap 363 DO jk = 1, jpkm1 364 DO jj = 1, jpj 365 DO ji = 1, jpi 366 ze3t_b = fse3t_b(ji,jj,jk) 367 ze3t_n = fse3t_n(ji,jj,jk) 368 ze3t_a = fse3t_a(ji,jj,jk) 369 ! ! tracer content at Before, now and after 370 ztcb = tb(ji,jj,jk) * ze3t_b ; zscb = sb(ji,jj,jk) * ze3t_b 371 ztcn = tn(ji,jj,jk) * ze3t_n ; zscn = sn(ji,jj,jk) * ze3t_n 372 ztca = ta(ji,jj,jk) * ze3t_a ; zsca = sa(ji,jj,jk) * ze3t_a 373 ! 374 ! ! Asselin filter on thickness and tracer content 375 ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b ) 376 ztc_f = atfp * ( ztca - 2.* ztcn + ztcb ) 377 zsc_f = atfp * ( zsca - 2.* zscn + zscb ) 378 ! 379 ! ! filtered tracer including the correction 380 ze3fr = 1.e0 / ( ze3t_n + ze3t_f ) 381 ztf = ( ztcn + ztc_f ) * ze3fr 382 zsf = ( zscn + zsc_f ) * ze3fr 383 ! ! swap of arrays 384 tb(ji,jj,jk) = ztf ! tb <-- tn filtered 385 sb(ji,jj,jk) = zsf 386 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 387 sn(ji,jj,jk) = sa(ji,jj,jk) 388 END DO 389 END DO 390 END DO 391 ENDIF 392 ENDIF 393 !!gm from Matthieu : unchecked 394 ! 395 IF( neuler /= 0 .OR. kt /= nit000 ) THEN ! remove the forcing from the Asselin filter 396 IF( lk_vvl ) THEN ! Varying levels 397 DO jj = 1, jpj 398 DO ji = 1, jpi 399 ! - ML - modification for varaiance diagnostics 400 zssh1 = tmask(ji,jj,jk) / ( fse3t(ji,jj,jk) + atfp * ( zfse3ta(ji,jj,jk) - 2*fse3t(ji,jj,jk) & 401 & + fse3tb(ji,jj,jk) ) ) 402 zt = zssh1 * ( tn(ji,jj,jk) + atfp * ( ta(ji,jj,jk) - 2 * tn(ji,jj,jk) + tb(ji,jj,jk) ) ) 403 zs = zssh1 * ( sn(ji,jj,jk) + atfp * ( sa(ji,jj,jk) - 2 * sn(ji,jj,jk) + sb(ji,jj,jk) ) ) 404 ! filter correction for global conservation 405 !------------------------------------------ 406 zfact1 = atfp * rdttra(1) * zssh1 407 IF (jk == 1) THEN ! remove sbc trend from time filter 408 zt = zt - zfact1 * ( qn(ji,jj) - qb(ji,jj) ) 409 !!??? zs = zs - zfact1 * ( - emp( ji,jj) + empb(ji,jj) ) * zsrau * 35.5e0 410 ENDIF 411 ! remove qsr trend from time filter 412 IF (jk <= nksr) THEN 413 IF( ln_sco ) THEN 414 z_cor_qsr = etot3(ji,jj,jk) * ( qsr(ji,jj) - qsrb(ji,jj) ) 415 ELSEIF( ln_zps .OR. ln_zco ) THEN 416 z_cor_qsr = ( qsr(ji,jj) - qsrb(ji,jj) ) * & 417 & ( gdsr(jk) * tmask(ji,jj,jk) - gdsr(jk+1) * tmask(ji,jj,jk+1) ) 418 ENDIF 419 zt = zt - zfact1 * z_cor_qsr 420 IF (jk == nksr) qsrb(ji,jj) = qsr(ji,jj) 421 ENDIF 422 ! - ML - end of modification 423 zssh1 = tmask(ji,jj,1) / zfse3ta(ji,jj,jk) 424 tn(ji,jj,jk) = ta(ji,jj,jk) * zssh1 425 sn(ji,jj,jk) = sa(ji,jj,jk) * zssh1 426 END DO 427 END DO 428 ENDIF 429 !!gm end 320 END DO 321 ENDIF 430 322 ! 431 323 END SUBROUTINE tra_nxt_vvl
Note: See TracChangeset
for help on using the changeset viewer.
