Changeset 1361
- Timestamp:
- 2009-03-31T18:26:41+02:00 (15 years ago)
- Location:
- branches/dev_004_VVL/NEMO/OPA_SRC
- Files:
-
- 2 edited
-
DYN/dynnxt.F90 (modified) (7 diffs)
-
TRA/tranxt.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynnxt.F90
r1152 r1361 4 4 !! Ocean dynamics: time stepping 5 5 !!====================================================================== 6 6 !!====================================================================== 7 !! History : OPA ! 1987-02 (P. Andrich, D. L Hostis) Original code 8 !! ! 1990-10 (C. Levy, G. Madec) 9 !! 7.0 ! 1993-03 (M. Guyon) symetrical conditions 10 !! 8.0 ! 1997-02 (G. Madec & M. Imbard) opa, release 8.0 11 !! 8.2 ! 1997-04 (A. Weaver) Euler forward step 12 !! - ! 1997-06 (G. Madec) lateral boudary cond., lbc routine 13 !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module 14 !! - ! 2002-10 (C. Talandier, A-M. Treguier) Open boundary cond. 15 !! 2.0 ! 2005-11 (V. Garnier) Surface pressure gradient organization 16 !! 2.3 ! 2007-07 (D. Storkey) Calls to BDY routines. 17 !! 3.1 ! 2009-02 (G. Madec) re-introduce the vvl option 18 !!---------------------------------------------------------------------- 19 7 20 !!---------------------------------------------------------------------- 8 21 !! dyn_nxt : update the horizontal velocity from the momentum trend … … 34 47 # include "domzgr_substitute.h90" 35 48 !!---------------------------------------------------------------------- 49 !! OPA 9.0 , LOCEAN-IPSL (2005) 50 !! $Id$ 51 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 52 !!---------------------------------------------------------------------- 36 53 37 54 CONTAINS … … 59 76 !! - Update un,vn arrays, the now horizontal velocity 60 77 !! 61 !! History :62 !! ! 87-02 (P. Andrich, D. L Hostis) Original code63 !! ! 90-10 (C. Levy, G. Madec)64 !! ! 93-03 (M. Guyon) symetrical conditions65 !! ! 97-02 (G. Madec & M. Imbard) opa, release 8.066 !! ! 97-04 (A. Weaver) Euler forward step67 !! ! 97-06 (G. Madec) lateral boudary cond., lbc routine68 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module69 !! ! 02-10 (C. Talandier, A-M. Treguier) Open boundary cond.70 !! 9.0 ! 05-11 (V. Garnier) Surface pressure gradient organization71 !! " ! 07-07 (D. Storkey) Calls to BDY routines.72 78 !!---------------------------------------------------------------------- 73 !! * Arguments74 79 INTEGER, INTENT( in ) :: kt ! ocean time-step index 75 76 !! * Local declarations 77 INTEGER :: ji, jj, jk ! dummy loop indices 80 !! 81 INTEGER :: jk ! dummy loop indices 78 82 REAL(wp) :: z2dt ! temporary scalar 79 REAL(wp) :: zsshun1, zsshvn1 ! temporary scalar80 !! Variable volume81 REAL(wp), DIMENSION(jpi,jpj) :: & ! 2D workspace82 zsshub, zsshun, zsshua, &83 zsshvb, zsshvn, zsshva84 REAL(wp), DIMENSION(jpi,jpj,jpk) :: &85 zfse3ub, zfse3un, zfse3ua, & ! 3D workspace86 zfse3vb, zfse3vn, zfse3va87 !!----------------------------------------------------------------------88 !! OPA 9.0 , LOCEAN-IPSL (2005)89 !! $Id$90 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt91 83 !!---------------------------------------------------------------------- 92 84 … … 102 94 103 95 !! Explicit physics with thickness weighted updates 104 IF( lk_vvl .AND. .NOT. lk_dynspg_flt ) THEN105 106 ! Sea surface elevation time stepping107 ! -----------------------------------108 !109 DO jj = 1, jpjm1110 DO ji = 1,jpim1111 112 ! Sea Surface Height at u-point before113 zsshub(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &114 & * ( e1t(ji ,jj ) * e2t(ji ,jj ) * sshbb(ji ,jj ) &115 & + e1t(ji+1,jj ) * e2t(ji+1,jj ) * sshbb(ji+1,jj ) )116 117 ! Sea Surface Height at v-point before118 zsshvb(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &119 & * ( e1t(ji ,jj ) * e2t(ji ,jj ) * sshbb(ji ,jj ) &120 & + e1t(ji ,jj+1) * e2t(ji ,jj+1) * sshbb(ji ,jj+1) )121 122 ! Sea Surface Height at u-point after123 zsshua(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &124 & * ( e1t(ji ,jj ) * e2t(ji ,jj ) * ssha(ji ,jj ) &125 & + e1t(ji+1,jj ) * e2t(ji+1,jj ) * ssha(ji+1,jj ) )126 127 ! Sea Surface Height at v-point after128 zsshva(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &129 & * ( e1t(ji ,jj ) * e2t(ji ,jj ) * ssha(ji ,jj ) &130 & + e1t(ji ,jj+1) * e2t(ji ,jj+1) * ssha(ji ,jj+1) )131 132 END DO133 END DO134 135 ! Boundaries conditions136 CALL lbc_lnk( zsshua, 'U', 1. ) ; CALL lbc_lnk( zsshva, 'V', 1. )137 CALL lbc_lnk( zsshub, 'U', 1. ) ; CALL lbc_lnk( zsshvb, 'V', 1. )138 139 ! Scale factors at before and after time step140 ! -------------------------------------------141 CALL dom_vvl_sf( zsshub, 'U', zfse3ub ) ; CALL dom_vvl_sf( zsshua, 'U', zfse3ua )142 CALL dom_vvl_sf( zsshvb, 'V', zfse3vb ) ; CALL dom_vvl_sf( zsshva, 'V', zfse3va )143 144 ! Asselin filtered scale factor at now time step145 ! ----------------------------------------------146 IF( (neuler == 0 .AND. kt == nit000) .OR. lk_dynspg_ts ) THEN147 CALL dom_vvl_sf_ini( 'U', zfse3un ) ; CALL dom_vvl_sf_ini( 'V', zfse3vn )148 ELSE149 zsshun(:,:) = atfp * ( zsshub(:,:) + zsshua(:,:) ) + atfp1 * sshu(:,:)150 zsshvn(:,:) = atfp * ( zsshvb(:,:) + zsshva(:,:) ) + atfp1 * sshv(:,:)151 CALL dom_vvl_sf( zsshun, 'U', zfse3un ) ; CALL dom_vvl_sf( zsshvn, 'V', zfse3vn )152 ENDIF153 154 ! Thickness weighting155 ! -------------------156 DO jk = 1, jpkm1157 DO jj = 1, jpj158 DO ji = 1, jpi159 ua(ji,jj,jk) = ua(ji,jj,jk) * fse3u(ji,jj,jk)160 va(ji,jj,jk) = va(ji,jj,jk) * fse3v(ji,jj,jk)161 162 un(ji,jj,jk) = un(ji,jj,jk) * fse3u(ji,jj,jk)163 vn(ji,jj,jk) = vn(ji,jj,jk) * fse3v(ji,jj,jk)164 165 ub(ji,jj,jk) = ub(ji,jj,jk) * zfse3ub(ji,jj,jk)166 vb(ji,jj,jk) = vb(ji,jj,jk) * zfse3vb(ji,jj,jk)167 END DO168 END DO169 END DO170 171 ENDIF172 96 173 97 ! Lateral boundary conditions on ( ua, va ) … … 175 99 CALL lbc_lnk( va, 'V', -1. ) 176 100 177 ! ! =============== 178 DO jk = 1, jpkm1 ! Horizontal slab 179 ! ! =============== 180 ! Next velocity 181 ! ------------- 101 ! Next velocity 102 ! ------------- 182 103 #if defined key_dynspg_flt 183 ! Leap-frog time stepping already done in dynspg.F routine104 ! Leap-frog time stepping already done in dynspg_flt.F routine 184 105 #else 185 DO jj = 1, jpj ! caution: don't use (:,:) for this loop 186 DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking 187 ! Leap-frog time stepping 188 ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk) 189 va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk) 190 END DO 191 END DO 192 193 IF( lk_vvl ) THEN 194 ! Unweight velocities prior to updating open boundaries. 195 196 DO jj = 1, jpj ! caution: don't use (:,:) for this loop 197 DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking 198 ua(ji,jj,jk) = ua(ji,jj,jk) / fse3u(ji,jj,jk) 199 va(ji,jj,jk) = va(ji,jj,jk) / fse3v(ji,jj,jk) 200 201 un(ji,jj,jk) = un(ji,jj,jk) / fse3u(ji,jj,jk) 202 vn(ji,jj,jk) = vn(ji,jj,jk) / fse3v(ji,jj,jk) 203 204 ub(ji,jj,jk) = ub(ji,jj,jk) / zfse3ub(ji,jj,jk) 205 vb(ji,jj,jk) = vb(ji,jj,jk) / zfse3vb(ji,jj,jk) 206 END DO 207 END DO 208 209 ENDIF 106 IF( lk_vvl ) THEN ! Varying levels 107 !RB_vvl scale factors are wrong at this point 108 DO jk = 1, jpkm1 109 ua(ji,jj,jk) = ( ub(:,:,jk) * fse3u(:,:,jk) & 110 & + z2dt * ua(:,:,jk) * fse3u(:,:,jk) ) & 111 & / fse3u(:,:,jk) * umask(:,:,jk) 112 va(ji,jj,jk) = ( vb(:,:,jk) * fse3v(:,:,jk) & 113 & + z2dt * va(:,:,jk) * fse3v(:,:,jk) ) & 114 & / fse3v(:,:,jk) * vmask(:,:,jk) 115 END DO 116 ELSE 117 DO jk = 1, jpkm1 118 ! Leap-frog time stepping 119 ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk) 120 va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk) 121 END DO 122 ENDIF 210 123 211 124 # if defined key_obc 212 ! ! ===============213 END DO ! End of slab214 ! ! ===============215 125 ! Update (ua,va) along open boundaries (only in the rigid-lid case) 216 126 CALL obc_dyn( kt ) … … 235 145 ENDIF 236 146 237 ! ! ===============238 DO jk = 1, jpkm1 ! Horizontal slab239 ! ! ===============240 147 # elif defined key_bdy 241 ! ! ===============242 END DO ! End of slab243 ! ! ===============244 148 ! Update (ua,va) along open boundaries (for exp or ts options). 245 149 IF ( lk_dynspg_exp .or. lk_dynspg_ts ) THEN … … 279 183 ENDIF 280 184 281 ! ! ===============282 DO jk = 1, jpkm1 ! Horizontal slab283 ! ! ===============284 185 # endif 285 186 # if defined key_agrif 286 ! ! ===============287 END DO ! End of slab288 ! ! ===============289 ! Update (ua,va) along open boundaries (only in the rigid-lid case)290 187 CALL Agrif_dyn( kt ) 291 ! ! ===============292 DO jk = 1, jpkm1 ! Horizontal slab293 ! ! ===============294 188 # endif 295 189 #endif 296 190 297 ! Time filter and swap of dynamics arrays 298 ! ------------------------------------------ 299 IF( neuler == 0 .AND. kt == nit000 ) THEN 300 IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN ! Varying levels 301 ! caution: don't use (:,:) for this loop 302 ! it causes optimization problems on NEC in auto-tasking 303 DO jj = 1, jpj 304 DO ji = 1, jpi 305 zsshun1 = umask(ji,jj,jk) / fse3u(ji,jj,jk) 306 zsshvn1 = vmask(ji,jj,jk) / fse3v(ji,jj,jk) 307 ub(ji,jj,jk) = un(ji,jj,jk) * zsshun1 * umask(ji,jj,jk) 308 vb(ji,jj,jk) = vn(ji,jj,jk) * zsshvn1 * vmask(ji,jj,jk) 309 zsshun1 = umask(ji,jj,jk) / zfse3ua(ji,jj,jk) 310 zsshvn1 = vmask(ji,jj,jk) / zfse3va(ji,jj,jk) 311 un(ji,jj,jk) = ua(ji,jj,jk) * zsshun1 * umask(ji,jj,jk) 312 vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn1 * vmask(ji,jj,jk) 313 END DO 314 END DO 315 ELSE ! Fixed levels 316 DO jj = 1, jpj 317 DO ji = 1, jpi 318 ! Euler (forward) time stepping 319 ub(ji,jj,jk) = un(ji,jj,jk) 320 vb(ji,jj,jk) = vn(ji,jj,jk) 321 un(ji,jj,jk) = ua(ji,jj,jk) 322 vn(ji,jj,jk) = va(ji,jj,jk) 323 END DO 324 END DO 325 ENDIF 326 ELSE 327 IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN ! Varying levels 328 ! caution: don't use (:,:) for this loop 329 ! it causes optimization problems on NEC in auto-tasking 330 DO jj = 1, jpj 331 DO ji = 1, jpi 332 zsshun1 = umask(ji,jj,jk) / zfse3un(ji,jj,jk) 333 zsshvn1 = vmask(ji,jj,jk) / zfse3vn(ji,jj,jk) 334 ub(ji,jj,jk) = ( atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) & 335 & + atfp1 * un(ji,jj,jk) ) * zsshun1 336 vb(ji,jj,jk) = ( atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) & 337 & + atfp1 * vn(ji,jj,jk) ) * zsshvn1 338 zsshun1 = umask(ji,jj,jk) / zfse3ua(ji,jj,jk) 339 zsshvn1 = vmask(ji,jj,jk) / zfse3va(ji,jj,jk) 340 un(ji,jj,jk) = ua(ji,jj,jk) * zsshun1 341 vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn1 342 END DO 343 END DO 344 ELSE ! Fixed levels 345 DO jj = 1, jpj 346 DO ji = 1, jpi 347 ! Leap-frog time stepping 348 ub(ji,jj,jk) = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk) 349 vb(ji,jj,jk) = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk) 350 un(ji,jj,jk) = ua(ji,jj,jk) 351 vn(ji,jj,jk) = va(ji,jj,jk) 352 END DO 353 END DO 354 ENDIF 191 ! Time filter and swap of dynamics arrays 192 ! ------------------------------------------ 193 IF( neuler == 0 .AND. kt == nit000 ) THEN 194 DO jk = 1, jpkm1 195 ub(:,:,jk) = un(:,:,jk) 196 vb(:,:,jk) = vn(:,:,jk) 197 un(:,:,jk) = ua(:,:,jk) 198 vn(:,:,jk) = va(:,:,jk) 199 END DO 200 ELSE 201 IF( lk_vvl ) THEN ! Varying levels 202 ! Not done 203 ELSE ! Fixed levels 204 !RB_vvl : should be done as in tranxt ? 205 DO jk = 1, jpkm1 ! filter applied on velocities 206 ub(:,:,jk) = atfp * ( ub(:,:,jk) + ua(:,:,jk) ) + atfp1 * un(:,:,jk) 207 vb(:,:,jk) = atfp * ( vb(:,:,jk) + va(:,:,jk) ) + atfp1 * vn(:,:,jk) 208 un(:,:,jk) = ua(:,:,jk) 209 vn(:,:,jk) = va(:,:,jk) 210 END DO 355 211 ENDIF 356 ! ! =============== 357 END DO ! End of slab 358 ! ! =============== 212 ENDIF 359 213 360 214 IF(ln_ctl) THEN -
branches/dev_004_VVL/NEMO/OPA_SRC/TRA/tranxt.F90
r1146 r1361 39 39 PUBLIC tra_nxt ! routine called by step.F90 40 40 41 REAL(wp), DIMENSION(jpk) :: r2dt_t ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler) 42 41 43 !! * Substitutions 42 44 # include "domzgr_substitute.h90" … … 67 69 !! at the AGRIF zoom boundaries (lk_agrif=T) 68 70 !! 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) 71 !! - Update lateral boundary conditions on AGRIF children 72 !! domains (lk_agrif=T) 80 73 !! 81 74 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step … … 87 80 INTEGER, INTENT(in) :: kt ! ocean time-step index 88 81 !! 89 INTEGER :: j i, jj, jk ! dummy loop indices90 REAL(wp) :: z t, zs, zfact ! temporary scalars82 INTEGER :: jk ! dummy loop indices 83 REAL(wp) :: zfact ! temporary scalars 91 84 !!---------------------------------------------------------------------- 92 85 … … 111 104 CALL Agrif_tra ! AGRIF zoom boundaries 112 105 #endif 106 107 ! set time step size (Euler/Leapfrog) 108 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt_t(:) = rdttra(:) ! at nit000 (Euler) 109 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt_t(:) = 2.* rdttra(:) ! at nit000 or nit000+1 (Leapfrog) 110 ENDIF 113 111 114 112 ! trends computation initialisation … … 118 116 ENDIF 119 117 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 ! 118 ! Leap-Frog + Asselin filter time stepping 119 IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt ) ! variable volume level (vvl) 120 ELSE ; CALL tra_nxt_fix( kt ) ! fixed volume level 160 121 ENDIF 161 122 … … 168 129 IF( l_trdtra ) THEN 169 130 DO jk = 1, jpkm1 170 zfact = 1.e0 / ( 2.*rdttra(jk) ) ! NB: euler case, (tb filtered - tb)=0 so 2dt always OK131 zfact = 1.e0 / r2dt_t(jk) 171 132 ztrdt(:,:,jk) = ( tb(:,:,jk) - ztrdt(:,:,jk) ) * zfact 172 133 ztrds(:,:,jk) = ( sb(:,:,jk) - ztrds(:,:,jk) ) * zfact … … 180 141 END SUBROUTINE tra_nxt 181 142 143 144 SUBROUTINE tra_nxt_fix( kt ) 145 !!---------------------------------------------------------------------- 146 !! *** ROUTINE tra_nxt_fix *** 147 !! 148 !! ** Purpose : fixed volume: apply the Asselin time filter and 149 !! swap the tracer fields. 150 !! 151 !! ** Method : - Apply a Asselin time filter on now fields. 152 !! - save in (ta,sa) an average over the three time levels 153 !! which will be used to compute rdn and thus the semi-implicit 154 !! hydrostatic pressure gradient (ln_dynhpg_imp = T) 155 !! - swap tracer fields to prepare the next time_step. 156 !! This can be summurized for tempearture as: 157 !! ztm = (ta+2tn+tb)/4 ln_dynhpg_imp = T 158 !! ztm = 0 otherwise 159 !! tb = tn + atfp*[ tb - 2 tn + ta ] 160 !! tn = ta 161 !! ta = ztm (NB: reset to 0 after eos_bn2 call) 162 !! 163 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step 164 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 165 !!---------------------------------------------------------------------- 166 INTEGER, INTENT(in) :: kt ! ocean time-step index 167 !! 168 INTEGER :: ji, jj, jk ! dummy loop indices 169 REAL(wp) :: ztm, ztf ! temporary scalars 170 REAL(wp) :: zsm, zsf ! - - 171 !!---------------------------------------------------------------------- 172 173 IF( kt == nit000 ) THEN 174 IF(lwp) WRITE(numout,*) 175 IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping' 176 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 177 ENDIF 178 ! 179 ! ! ----------------------- ! 180 IF( ln_dynhpg_imp ) THEN ! semi-implicite hpg case ! 181 ! ! ----------------------- ! 182 ! 183 IF( neuler == 0 .AND. kt == nit000 ) THEN ! case of Euler time-stepping at first time-step 184 DO jk = 1, jpkm1 185 DO jj = 1, jpj 186 DO ji = 1, jpi 187 ztm = 0.25 * ( ta(ji,jj,jk) + 2. * tn(ji,jj,jk) + tb(ji,jj,jk) ) ! mean t 188 zsm = 0.25 * ( sa(ji,jj,jk) + 2. * sn(ji,jj,jk) + sb(ji,jj,jk) ) 189 tb(ji,jj,jk) = tn(ji,jj,jk) ! tb <-- tn 190 sb(ji,jj,jk) = sn(ji,jj,jk) 191 tn(ji,jj,jk) = ta(ji,jj,jk) ! tb <-- tn 192 sn(ji,jj,jk) = sa(ji,jj,jk) 193 ta(ji,jj,jk) = ztm ! ta <-- mean t 194 sa(ji,jj,jk) = zsm 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 211 sa(ji,jj,jk) = zsm 212 END DO 213 END DO 214 END DO 215 ENDIF 216 ! ! ----------------------- ! 217 ELSE ! explicit hpg case ! 218 ! ! ----------------------- ! 219 ! 220 IF( neuler == 0 .AND. kt == nit000 ) THEN ! case of Euler time-stepping at first time-step 221 DO jk = 1, jpkm1 222 DO jj = 1, jpj 223 DO ji = 1, jpi 224 tb(ji,jj,jk) = tn(ji,jj,jk) ! tb <-- tn 225 sb(ji,jj,jk) = sn(ji,jj,jk) 226 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 227 sn(ji,jj,jk) = sa(ji,jj,jk) 228 END DO 229 END DO 230 END DO 231 ELSE ! general case (Leapfrog + Asselin filter 232 DO jk = 1, jpkm1 233 DO jj = 1, jpj 234 DO ji = 1, jpi 235 !RBvvl for reproducibility 236 ! ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) ) ! Asselin filter on t 237 ! zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 238 ! tb(ji,jj,jk) = tn(ji,jj,jk) + ztf ! tb <-- filtered tn 239 ! sb(ji,jj,jk) = sn(ji,jj,jk) + zsf 240 tb(ji,jj,jk) = atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk) 241 sb(ji,jj,jk) = atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk) 242 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 243 sn(ji,jj,jk) = sa(ji,jj,jk) 244 END DO 245 END DO 246 END DO 247 ENDIF 248 ENDIF 249 ! 250 END SUBROUTINE tra_nxt_fix 251 252 253 SUBROUTINE tra_nxt_vvl( kt ) 254 !!---------------------------------------------------------------------- 255 !! *** ROUTINE tra_nxt_vvl *** 256 !! 257 !! ** Purpose : Time varying volume: apply the Asselin time filter 258 !! and swap the tracer fields. 259 !! 260 !! ** Method : - Apply a thickness weighted Asselin time filter on now fields. 261 !! - save in (ta,sa) a thickness weighted average over the three 262 !! time levels which will be used to compute rdn and thus the semi- 263 !! implicit hydrostatic pressure gradient (ln_dynhpg_imp = T) 264 !! - swap tracer fields to prepare the next time_step. 265 !! This can be summurized for tempearture as: 266 !! ztm = (e3t_a*ta+2*e3t_n*tn+e3t_b*tb) ln_dynhpg_imp = T 267 !! /(e3t_a +2*e3t_n +e3t_b ) 268 !! ztm = 0 otherwise 269 !! tb = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) 270 !! /( e3t_n + atfp*[ e3t_b - 2 e3t_n + e3t_a ] ) 271 !! tn = ta 272 !! ta = zt (NB: reset to 0 after eos_bn2 call) 273 !! 274 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step 275 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 276 !!---------------------------------------------------------------------- 277 INTEGER, INTENT(in) :: kt ! ocean time-step index 278 !! 279 280 ! Not yet ready 281 WRITE(*,*) 'tra_next_vvl : you should not be there' 282 STOP 283 ! 284 END SUBROUTINE tra_nxt_vvl 285 182 286 !!====================================================================== 183 287 END MODULE tranxt
Note: See TracChangeset
for help on using the changeset viewer.
