Changeset 1847
- Timestamp:
- 2010-04-21T15:38:47+02:00 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90
r1601 r1847 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-04 (M. Leclair, G. Madec) semi-implicit hpg with asselin filter 17 18 !!---------------------------------------------------------------------- 18 19 … … 44 45 PUBLIC tra_nxt ! routine called by step.F90 45 46 46 REAL(wp), DIMENSION(jpk) :: r2dt_t ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler) 47 REAL(wp) :: rbcp, r1_2bcp ! Brown & Campana parameters for semi-implicit hpg 48 REAL(wp), DIMENSION(jpk) :: r2dt_t ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler) 47 49 48 50 !! * Substitutions 49 51 # include "domzgr_substitute.h90" 50 52 !!---------------------------------------------------------------------- 51 !! NEMO/OPA 3. 0 , LOCEAN-IPSL (2008)53 !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 52 54 !! $Id$ 53 55 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 85 87 INTEGER, INTENT(in) :: kt ! ocean time-step index 86 88 !! 87 INTEGER :: jk ! dummy loop indices88 REAL(wp) :: zfact ! temporary scalars89 INTEGER :: jk ! dummy loop indices 90 REAL(wp) :: zfact ! temporary scalars 89 91 !!---------------------------------------------------------------------- 90 92 … … 93 95 IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap' 94 96 IF(lwp) WRITE(numout,*) '~~~~~~~' 97 ! 98 rbcp = 0.25 * (1 + atfp) * (1 + atfp * atfp) ! Brown & Campana parameter for semi-implicit hpg 99 r1_2bcp = 1.e0 - 2.e0 * rbcp 95 100 ENDIF 96 101 … … 160 165 !! hydrostatic pressure gradient (ln_dynhpg_imp = T) 161 166 !! - swap tracer fields to prepare the next time_step. 162 !! This can be summurized for tempearture as: 163 !! ztm = (ta+2tn+tb)/4 ln_dynhpg_imp = T 164 !! ztm = 0 otherwise 167 !! This can be summurized for temperature as: 168 !! ztm = (rbcp*ta+(2-rbcp)*tn+rbcp*tb) ln_dynhpg_imp = T 169 !! ztm = 0 otherwise 170 !! with rbcp=(1+atfp)*(1+atfp*atfp)/4 165 171 !! tb = tn + atfp*[ tb - 2 tn + ta ] 166 172 !! tn = ta … … 189 195 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 190 196 DO jk = 1, jpkm1 ! (only swap) 197 tn(:,:,jk) = ta(:,:,jk) ! tb <-- tn 198 sn(:,:,jk) = sa(:,:,jk) 199 END DO 200 ELSE ! general case (Leapfrog + Asselin filter) 201 DO jk = 1, jpkm1 191 202 DO jj = 1, jpj 192 203 DO ji = 1, jpi 193 tn(ji,jj,jk) = ta(ji,jj,jk) ! tb <-- tn 204 ztm = rbcp * ( ta(ji,jj,jk) + tb(ji,jj,jk) ) + r1_2bcp * tn(ji,jj,jk) ! mean t 205 zsm = rbcp * ( sa(ji,jj,jk) + sb(ji,jj,jk) ) + r1_2bcp * sn(ji,jj,jk) 206 ztf = atfp * ( ta(ji,jj,jk) + tb(ji,jj,jk) - 2. * tn(ji,jj,jk) ) ! Asselin filter on t 207 zsf = atfp * ( sa(ji,jj,jk) + sb(ji,jj,jk) - 2. * sn(ji,jj,jk) ) 208 tb(ji,jj,jk) = tn(ji,jj,jk) + ztf ! tb <-- filtered tn 209 sb(ji,jj,jk) = sn(ji,jj,jk) + zsf 210 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 194 211 sn(ji,jj,jk) = sa(ji,jj,jk) 195 END DO 196 END DO 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 = 0.25 * ( ta(ji,jj,jk) + 2.* tn(ji,jj,jk) + tb(ji,jj,jk) ) ! mean t 203 zsm = 0.25 * ( sa(ji,jj,jk) + 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 204 ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) ) ! Asselin filter on t 205 zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(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 212 ta(ji,jj,jk) = ztm ! ta <-- mean t 211 213 sa(ji,jj,jk) = zsm 212 214 END DO … … 220 222 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 221 223 DO jk = 1, jpkm1 224 tn(:,:,jk) = ta(:,:,jk) ! tn <-- ta 225 sn(:,:,jk) = sa(:,:,jk) 226 END DO 227 ELSE ! general case (Leapfrog + Asselin filter 228 DO jk = 1, jpkm1 222 229 DO jj = 1, jpj 223 230 DO ji = 1, jpi 224 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 225 sn(ji,jj,jk) = sa(ji,jj,jk) 226 END DO 227 END DO 228 END DO 229 ELSE ! general case (Leapfrog + Asselin filter 230 DO jk = 1, jpkm1 231 DO jj = 1, jpj 232 DO ji = 1, jpi 233 ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) ) ! Asselin filter on t 231 ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) ) ! Asselin filter on t 234 232 zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 235 tb(ji,jj,jk) = tn(ji,jj,jk) + ztf 233 tb(ji,jj,jk) = tn(ji,jj,jk) + ztf ! tb <-- filtered tn 236 234 sb(ji,jj,jk) = sn(ji,jj,jk) + zsf 237 tn(ji,jj,jk) = ta(ji,jj,jk) 235 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 238 236 sn(ji,jj,jk) = sa(ji,jj,jk) 239 237 END DO … … 258 256 !! implicit hydrostatic pressure gradient (ln_dynhpg_imp = T) 259 257 !! - swap tracer fields to prepare the next time_step. 260 !! This can be summurized for tempe arture as:261 !! ztm = ( e3t_a*ta+2*e3t_n*tn+e3t_b*tb) ln_dynhpg_imp = T262 !! /( e3t_a +2*e3t_n +e3t_b )263 !! ztm = 0 otherwise258 !! This can be summurized for temperature as: 259 !! ztm = ( rbcp*e3t_a*ta + (2-rbtp)*e3t_n*tn + rbcp*e3t_b*tb ) ln_dynhpg_imp = T 260 !! /( rbcp*e3t_a + (2-rbcp)*e3t_n + rbcp*e3t_b ) 261 !! ztm = 0 otherwise 264 262 !! tb = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) 265 263 !! /( e3t_n + atfp*[ e3t_b - 2 e3t_n + e3t_a ] ) … … 286 284 287 285 ! ! ----------------------- ! 288 IF( ln_dynhpg_imp ) THEN ! semi-implicite hpg case ! 286 IF( ln_dynhpg_imp ) THEN ! semi-implicite hpg case ! (mean tracer computation) 289 287 ! ! ----------------------- ! 290 288 ! 291 289 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 292 DO jk = 1, jpkm1 ! (only swap) 293 DO jj = 1, jpj 294 DO ji = 1, jpi 295 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 296 sn(ji,jj,jk) = sa(ji,jj,jk) 297 END DO 298 END DO 299 END DO 300 ELSE 290 DO jk = 1, jpkm1 ! (only swap) 291 tn(:,:,jk) = ta(:,:,jk) ! tn <-- ta 292 sn(:,:,jk) = sa(:,:,jk) 293 END DO 294 ! ! general case (Leapfrog + Asselin filter) 295 ELSE ! apply filter on thickness weighted tracer and swap 301 296 DO jk = 1, jpkm1 302 297 DO jj = 1, jpj … … 311 306 ! 312 307 ! ! Asselin filter on thickness and tracer content 313 ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b)314 ztc_f = atfp * ( ztca - 2.* ztcn + ztcb)315 zsc_f = atfp * ( zsca - 2.* zscn + zscb)308 ze3t_f = atfp * ( ze3t_a + ze3t_b - 2.* ze3t_n ) 309 ztc_f = atfp * ( ztca + ztcb - 2.* ztcn ) 310 zsc_f = atfp * ( zsca + zscb - 2.* zscn ) 316 311 ! 317 312 ! ! filtered tracer including the correction … … 320 315 zsf = ze3fr * ( zscn + zsc_f ) 321 316 ! ! mean thickness and tracer 322 ze3mr = 1.e0 / ( ze3t_a + 2.* ze3t_n + ze3t_b)323 ztm = ze3mr * ( ztca + 2.* ztcn + ztcb)324 zsm = ze3mr * ( zsca + 2.* zscn + zscb)317 ze3mr = 1.e0 / ( rbcp * ( ze3t_a + ze3t_b ) + r1_2bcp * ze3t_n ) 318 ztm = ze3mr * ( rbcp * ( ztca + ztcb ) + r1_2bcp * ztcn ) 319 zsm = ze3mr * ( rbcp * ( zsca + zscb ) + r1_2bcp * zscn ) 325 320 !!gm mean e3t have to be saved and used in dynhpg or it can be recomputed in dynhpg !! 326 !!gm e3t_m(ji,jj,jk) = 0.25/ ze3mr321 !!gm e3t_m(ji,jj,jk) = 1.e0 / ze3mr 327 322 ! ! swap of arrays 328 323 tb(ji,jj,jk) = ztf ! tb <-- tn + filter … … 337 332 ENDIF 338 333 ! ! ----------------------- ! 339 ELSE ! explicit hpg case ! 334 ELSE ! explicit hpg case ! (NO mean tracer computation) 340 335 ! ! ----------------------- ! 341 336 ! 342 337 IF( neuler == 0 .AND. kt == nit000 ) THEN ! case of Euler time-stepping at first time-step 343 DO jk = 1, jpkm1 ! No filter nor thickness weighting computation required 344 DO jj = 1, jpj ! ONLY swap 345 DO ji = 1, jpi 346 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 347 sn(ji,jj,jk) = sa(ji,jj,jk) 348 END DO 349 END DO 338 DO jk = 1, jpkm1 ! (only swap) 339 tn(:,:,jk) = ta(:,:,jk) ! tn <-- ta 340 sn(:,:,jk) = sa(:,:,jk) 350 341 END DO 351 342 ! ! general case (Leapfrog + Asselin filter)
Note: See TracChangeset
for help on using the changeset viewer.
