Changeset 1438 for trunk/NEMO/OPA_SRC/TRA/tranxt.F90
- Timestamp:
- 2009-05-11T16:34:47+02:00 (15 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/TRA/tranxt.F90 (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/TRA/tranxt.F90
r1146 r1438 14 14 !! 2.0 ! 2006-02 (L. Debreu, C. Mazauric) Agrif implementation 15 15 !! 3.0 ! 2008-06 (G. Madec) time stepping always done in trazdf 16 !! 3.1 ! 2009-02 (G. Madec, R. Benshila) re-introduce the vvl option 16 17 !!---------------------------------------------------------------------- 17 18 18 19 !!---------------------------------------------------------------------- 19 20 !! tra_nxt : time stepping on temperature and salinity 21 !! tra_nxt_fix : time stepping on temperature and salinity : fixed volume case 22 !! tra_nxt_vvl : time stepping on temperature and salinity : variable volume case 20 23 !!---------------------------------------------------------------------- 21 24 USE oce ! ocean dynamics and tracers variables 22 25 USE dom_oce ! ocean space and time domain variables 23 26 USE zdf_oce ! ??? 27 USE domvvl ! variable volume 24 28 USE dynspg_oce ! surface pressure gradient variables 25 29 USE trdmod_oce ! ocean variables trends … … 39 43 PUBLIC tra_nxt ! routine called by step.F90 40 44 45 REAL(wp), DIMENSION(jpk) :: r2dt_t ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler) 46 41 47 !! * Substitutions 42 48 # include "domzgr_substitute.h90" … … 67 73 !! at the AGRIF zoom boundaries (lk_agrif=T) 68 74 !! 69 !! - Apply the Asselin time filter on now fields, 70 !! save in (ta,sa) an average over the three time levels 71 !! which will be used to compute rdn and thus the semi-implicit 72 !! hydrostatic pressure gradient (ln_dynhpg_imp = T), and 73 !! swap tracer fields to prepare the next time_step. 74 !! This can be summurized for tempearture as: 75 !! zt = (ta+2tn+tb)/4 ln_dynhpg_imp = T 76 !! zt = 0 otherwise 77 !! tb = tn + atfp*[ tb - 2 tn + ta ] 78 !! tn = ta 79 !! ta = zt (NB: reset to 0 after eos_bn2 call) 75 !! - Update lateral boundary conditions on AGRIF children 76 !! domains (lk_agrif=T) 80 77 !! 81 78 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step … … 87 84 INTEGER, INTENT(in) :: kt ! ocean time-step index 88 85 !! 89 INTEGER :: j i, jj, jk ! dummy loop indices90 REAL(wp) :: z t, zs, zfact ! temporary scalars86 INTEGER :: jk ! dummy loop indices 87 REAL(wp) :: zfact ! temporary scalars 91 88 !!---------------------------------------------------------------------- 92 89 … … 111 108 CALL Agrif_tra ! AGRIF zoom boundaries 112 109 #endif 110 111 ! set time step size (Euler/Leapfrog) 112 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt_t(:) = rdttra(:) ! at nit000 (Euler) 113 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt_t(:) = 2.* rdttra(:) ! at nit000 or nit000+1 (Leapfrog) 114 ENDIF 113 115 114 116 ! trends computation initialisation … … 118 120 ENDIF 119 121 120 ! Asselin time filter and swap of arrays 121 ! 122 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler 1st time step : swap only 123 DO jk = 1, jpkm1 124 tb(:,:,jk) = tn(:,:,jk) ! ta, sa remain at their values which 125 sb(:,:,jk) = sn(:,:,jk) ! correspond to tn, sn after the sawp 126 tn(:,:,jk) = ta(:,:,jk) 127 sn(:,:,jk) = sa(:,:,jk) 128 END DO 129 ! 130 ELSE ! Leap-Frog : filter + swap 131 ! 132 IF( ln_dynhpg_imp ) THEN ! semi-implicite hpg case 133 DO jk = 1, jpkm1 ! (save the averaged of the 3 time steps 134 DO jj = 1, jpj ! in the after fields) 135 DO ji = 1, jpi 136 zt = ( ta(ji,jj,jk) + 2. * tn(ji,jj,jk) + tb(ji,jj,jk) ) * 0.25 137 zs = ( sa(ji,jj,jk) + 2. * sn(ji,jj,jk) + sb(ji,jj,jk) ) * 0.25 138 tb(ji,jj,jk) = atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk) 139 sb(ji,jj,jk) = atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk) 140 tn(ji,jj,jk) = ta(ji,jj,jk) 141 sn(ji,jj,jk) = sa(ji,jj,jk) 142 ta(ji,jj,jk) = zt 143 sa(ji,jj,jk) = zs 144 END DO 145 END DO 146 END DO 147 ELSE ! explicit hpg case 148 DO jk = 1, jpkm1 149 DO jj = 1, jpj 150 DO ji = 1, jpi 151 tb(ji,jj,jk) = atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk) 152 sb(ji,jj,jk) = atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk) 153 tn(ji,jj,jk) = ta(ji,jj,jk) 154 sn(ji,jj,jk) = sa(ji,jj,jk) 155 END DO 156 END DO 157 END DO 158 ENDIF 159 ! 122 ! Leap-Frog + Asselin filter time stepping 123 IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt ) ! variable volume level (vvl) 124 ELSE ; CALL tra_nxt_fix( kt ) ! fixed volume level 160 125 ENDIF 161 126 … … 165 130 #endif 166 131 167 ! trends diagnostics : Asselin filter trend : (tb filtered - tb)/2dt168 IF( l_trdtra ) THEN 132 ! trends computation 133 IF( l_trdtra ) THEN ! trend of the Asselin filter (tb filtered - tb)/dt 169 134 DO jk = 1, jpkm1 170 zfact = 1.e0 / ( 2.*rdttra(jk) ) ! NB: euler case, (tb filtered - tb)=0 so 2dt always OK135 zfact = 1.e0 / r2dt_t(jk) 171 136 ztrdt(:,:,jk) = ( tb(:,:,jk) - ztrdt(:,:,jk) ) * zfact 172 137 ztrds(:,:,jk) = ( sb(:,:,jk) - ztrds(:,:,jk) ) * zfact … … 174 139 CALL trd_mod( ztrdt, ztrds, jptra_trd_atf, 'TRA', kt ) 175 140 END IF 141 176 142 ! ! control print 177 143 IF(ln_ctl) CALL prt_ctl( tab3d_1=tn, clinfo1=' nxt - Tn: ', mask1=tmask, & … … 180 146 END SUBROUTINE tra_nxt 181 147 148 149 SUBROUTINE tra_nxt_fix( kt ) 150 !!---------------------------------------------------------------------- 151 !! *** ROUTINE tra_nxt_fix *** 152 !! 153 !! ** Purpose : fixed volume: apply the Asselin time filter and 154 !! swap the tracer fields. 155 !! 156 !! ** Method : - Apply a Asselin time filter on now fields. 157 !! - save in (ta,sa) an average over the three time levels 158 !! which will be used to compute rdn and thus the semi-implicit 159 !! hydrostatic pressure gradient (ln_dynhpg_imp = T) 160 !! - swap tracer fields to prepare the next time_step. 161 !! This can be summurized for tempearture as: 162 !! ztm = (ta+2tn+tb)/4 ln_dynhpg_imp = T 163 !! ztm = 0 otherwise 164 !! tb = tn + atfp*[ tb - 2 tn + ta ] 165 !! tn = ta 166 !! ta = ztm (NB: reset to 0 after eos_bn2 call) 167 !! 168 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step 169 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 170 !!---------------------------------------------------------------------- 171 INTEGER, INTENT(in) :: kt ! ocean time-step index 172 !! 173 INTEGER :: ji, jj, jk ! dummy loop indices 174 REAL(wp) :: ztm, ztf ! temporary scalars 175 REAL(wp) :: zsm, zsf ! - - 176 !!---------------------------------------------------------------------- 177 178 IF( kt == nit000 ) THEN 179 IF(lwp) WRITE(numout,*) 180 IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping' 181 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 182 ENDIF 183 ! 184 ! ! ----------------------- ! 185 IF( ln_dynhpg_imp ) THEN ! semi-implicite hpg case ! 186 ! ! ----------------------- ! 187 ! 188 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 189 DO jk = 1, jpkm1 ! (only swap) 190 DO jj = 1, jpj 191 DO ji = 1, jpi 192 tn(ji,jj,jk) = ta(ji,jj,jk) ! tb <-- tn 193 sn(ji,jj,jk) = sa(ji,jj,jk) 194 END DO 195 END DO 196 END DO 197 ELSE ! general case (Leapfrog + Asselin filter 198 DO jk = 1, jpkm1 199 DO jj = 1, jpj 200 DO ji = 1, jpi 201 ztm = 0.25 * ( ta(ji,jj,jk) + 2.* tn(ji,jj,jk) + tb(ji,jj,jk) ) ! mean t 202 zsm = 0.25 * ( sa(ji,jj,jk) + 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 203 ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) ) ! Asselin filter on t 204 zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 205 tb(ji,jj,jk) = tn(ji,jj,jk) + ztf ! tb <-- filtered tn 206 sb(ji,jj,jk) = sn(ji,jj,jk) + zsf 207 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 208 sn(ji,jj,jk) = sa(ji,jj,jk) 209 ta(ji,jj,jk) = ztm ! ta <-- mean t 210 sa(ji,jj,jk) = zsm 211 END DO 212 END DO 213 END DO 214 ENDIF 215 ! ! ----------------------- ! 216 ELSE ! explicit hpg case ! 217 ! ! ----------------------- ! 218 ! 219 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 220 DO jk = 1, jpkm1 221 DO jj = 1, jpj 222 DO ji = 1, jpi 223 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 224 sn(ji,jj,jk) = sa(ji,jj,jk) 225 END DO 226 END DO 227 END DO 228 ELSE ! general case (Leapfrog + Asselin filter 229 DO jk = 1, jpkm1 230 DO jj = 1, jpj 231 DO ji = 1, jpi 232 ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) ) ! Asselin filter on t 233 zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 234 tb(ji,jj,jk) = tn(ji,jj,jk) + ztf ! tb <-- filtered tn 235 sb(ji,jj,jk) = sn(ji,jj,jk) + zsf 236 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 237 sn(ji,jj,jk) = sa(ji,jj,jk) 238 END DO 239 END DO 240 END DO 241 ENDIF 242 ENDIF 243 ! 244 END SUBROUTINE tra_nxt_fix 245 246 247 SUBROUTINE tra_nxt_vvl( kt ) 248 !!---------------------------------------------------------------------- 249 !! *** ROUTINE tra_nxt_vvl *** 250 !! 251 !! ** Purpose : Time varying volume: apply the Asselin time filter 252 !! and swap the tracer fields. 253 !! 254 !! ** Method : - Apply a thickness weighted Asselin time filter on now fields. 255 !! - save in (ta,sa) a thickness weighted average over the three 256 !! time levels which will be used to compute rdn and thus the semi- 257 !! implicit hydrostatic pressure gradient (ln_dynhpg_imp = T) 258 !! - swap tracer fields to prepare the next time_step. 259 !! This can be summurized for tempearture as: 260 !! ztm = (e3t_a*ta+2*e3t_n*tn+e3t_b*tb) ln_dynhpg_imp = T 261 !! /(e3t_a +2*e3t_n +e3t_b ) 262 !! ztm = 0 otherwise 263 !! tb = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) 264 !! /( e3t_n + atfp*[ e3t_b - 2 e3t_n + e3t_a ] ) 265 !! tn = ta 266 !! ta = zt (NB: reset to 0 after eos_bn2 call) 267 !! 268 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step 269 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 270 !!---------------------------------------------------------------------- 271 INTEGER, INTENT(in) :: kt ! ocean time-step index 272 !! 273 INTEGER :: ji, jj, jk ! dummy loop indices 274 REAL(wp) :: ztm , ztc_f , ztf , ztca, ztcn, ztcb ! temporary scalar 275 REAL(wp) :: zsm , zsc_f , zsf , zsca, zscn, zscb ! - - 276 REAL(wp) :: ze3mr, ze3fr ! - - 277 REAL(wp) :: ze3t_b, ze3t_n, ze3t_a, ze3t_f ! - - 278 !!---------------------------------------------------------------------- 279 280 IF( kt == nit000 ) THEN 281 IF(lwp) WRITE(numout,*) 282 IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping' 283 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 284 ENDIF 285 286 ! ! ----------------------- ! 287 IF( ln_dynhpg_imp ) THEN ! semi-implicite hpg case ! 288 ! ! ----------------------- ! 289 ! 290 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 291 DO jk = 1, jpkm1 ! (only swap) 292 DO jj = 1, jpj 293 DO ji = 1, jpi 294 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 295 sn(ji,jj,jk) = sa(ji,jj,jk) 296 END DO 297 END DO 298 END DO 299 ELSE 300 DO jk = 1, jpkm1 301 DO jj = 1, jpj 302 DO ji = 1, jpi 303 ze3t_b = fse3t_b(ji,jj,jk) 304 ze3t_n = fse3t_n(ji,jj,jk) 305 ze3t_a = fse3t_a(ji,jj,jk) 306 ! ! tracer content at Before, now and after 307 ztcb = tb(ji,jj,jk) * ze3t_b ; zscb = sb(ji,jj,jk) * ze3t_b 308 ztcn = tn(ji,jj,jk) * ze3t_n ; zscn = sn(ji,jj,jk) * ze3t_n 309 ztca = ta(ji,jj,jk) * ze3t_a ; zsca = sa(ji,jj,jk) * ze3t_a 310 ! 311 ! ! Asselin filter on thickness and tracer content 312 ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b ) 313 ztc_f = atfp * ( ztca - 2.* ztcn + ztcb ) 314 zsc_f = atfp * ( zsca - 2.* zscn + zscb ) 315 ! 316 ! ! filtered tracer including the correction 317 ze3fr = 1.e0 / ( ze3t_n + ze3t_f ) 318 ztf = ze3fr * ( ztcn + ztc_f ) 319 zsf = ze3fr * ( zscn + zsc_f ) 320 ! ! mean thickness and tracer 321 ze3mr = 1.e0 / ( ze3t_a + 2.* ze3t_n + ze3t_b ) 322 ztm = ze3mr * ( ztca + 2.* ztcn + ztcb ) 323 zsm = ze3mr * ( zsca + 2.* zscn + zscb ) 324 !!gm mean e3t have to be saved and used in dynhpg or it can be recomputed in dynhpg !! 325 !!gm e3t_m(ji,jj,jk) = 0.25 / ze3mr 326 ! ! swap of arrays 327 tb(ji,jj,jk) = ztf ! tb <-- tn + filter 328 sb(ji,jj,jk) = zsf 329 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 330 sn(ji,jj,jk) = sa(ji,jj,jk) 331 ta(ji,jj,jk) = ztm ! ta <-- mean t 332 sa(ji,jj,jk) = zsm 333 END DO 334 END DO 335 END DO 336 ENDIF 337 ! ! ----------------------- ! 338 ELSE ! explicit hpg case ! 339 ! ! ----------------------- ! 340 ! 341 IF( neuler == 0 .AND. kt == nit000 ) THEN ! case of Euler time-stepping at first time-step 342 DO jk = 1, jpkm1 ! No filter nor thickness weighting computation required 343 DO jj = 1, jpj ! ONLY swap 344 DO ji = 1, jpi 345 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 346 sn(ji,jj,jk) = sa(ji,jj,jk) 347 END DO 348 END DO 349 END DO 350 ! ! general case (Leapfrog + Asselin filter) 351 ELSE ! apply filter on thickness weighted tracer and swap 352 DO jk = 1, jpkm1 353 DO jj = 1, jpj 354 DO ji = 1, jpi 355 ze3t_b = fse3t_b(ji,jj,jk) 356 ze3t_n = fse3t_n(ji,jj,jk) 357 ze3t_a = fse3t_a(ji,jj,jk) 358 ! ! tracer content at Before, now and after 359 ztcb = tb(ji,jj,jk) * ze3t_b ; zscb = sb(ji,jj,jk) * ze3t_b 360 ztcn = tn(ji,jj,jk) * ze3t_n ; zscn = sn(ji,jj,jk) * ze3t_n 361 ztca = ta(ji,jj,jk) * ze3t_a ; zsca = sa(ji,jj,jk) * ze3t_a 362 ! 363 ! ! Asselin filter on thickness and tracer content 364 ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b ) 365 ztc_f = atfp * ( ztca - 2.* ztcn + ztcb ) 366 zsc_f = atfp * ( zsca - 2.* zscn + zscb ) 367 ! 368 ! ! filtered tracer including the correction 369 ze3fr = 1.e0 / ( ze3t_n + ze3t_f ) 370 ztf = ( ztcn + ztc_f ) * ze3fr 371 zsf = ( zscn + zsc_f ) * ze3fr 372 ! ! swap of arrays 373 tb(ji,jj,jk) = ztf ! tb <-- tn filtered 374 sb(ji,jj,jk) = zsf 375 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 376 sn(ji,jj,jk) = sa(ji,jj,jk) 377 END DO 378 END DO 379 END DO 380 ENDIF 381 ENDIF 382 ! 383 END SUBROUTINE tra_nxt_vvl 384 182 385 !!====================================================================== 183 386 END MODULE tranxt
Note: See TracChangeset
for help on using the changeset viewer.
