Changeset 2528 for trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
- Timestamp:
- 2010-12-27T18:33:53+01:00 (13 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90 (modified) (9 diffs, 1 prop)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
- Property svn:eol-style deleted
r1876 r2528 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 !!---------------------------------------------------------------------- 18 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 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 19 !!---------------------------------------------------------------------- 20 21 !!---------------------------------------------------------------------- 22 !! 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 23 25 !!---------------------------------------------------------------------- 24 26 USE oce ! ocean dynamics and tracers variables 25 27 USE dom_oce ! ocean space and time domain variables 28 USE sbc_oce ! surface boundary condition: ocean 26 29 USE zdf_oce ! ??? 27 30 USE domvvl ! variable volume 28 31 USE dynspg_oce ! surface pressure gradient variables 29 32 USE dynhpg ! hydrostatic pressure gradient 30 USE trdmod_oce ! ocean variables trends31 USE trd mod! ocean active tracers trends33 USE trdmod_oce ! ocean space and time domain variables 34 USE trdtra ! ocean active tracers trends 32 35 USE phycst 36 USE obc_oce 33 37 USE obctra ! open boundary condition (obc_tra routine) 34 USE bdytra ! Unstructured open boundary condition (bdy_tra routine) 38 USE bdy_par ! Unstructured open boundary condition (bdy_tra_frs routine) 39 USE bdytra ! Unstructured open boundary condition (bdy_tra_frs routine) 35 40 USE in_out_manager ! I/O manager 36 41 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 37 42 USE prtctl ! Print control 43 USE traqsr ! penetrative solar radiation (needed for nksr) 44 USE traswp ! swap array 45 USE obc_oce 46 #if defined key_agrif 38 47 USE agrif_opa_update 39 48 USE agrif_opa_interp 40 USE obc_oce 49 #endif 41 50 42 51 IMPLICIT NONE 43 52 PRIVATE 44 53 45 PUBLIC tra_nxt ! routine called by step.F90 46 47 REAL(wp), DIMENSION(jpk) :: r2dt_t ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler) 54 PUBLIC tra_nxt ! routine called by step.F90 55 PUBLIC tra_nxt_fix ! to be used in trcnxt 56 PUBLIC tra_nxt_vvl ! to be used in trcnxt 57 58 REAL(wp) :: rbcp ! Brown & Campana parameters for semi-implicit hpg 59 REAL(wp), DIMENSION(jpk) :: r2dt ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler) 48 60 49 61 !! * Substitutions 50 62 # include "domzgr_substitute.h90" 51 63 !!---------------------------------------------------------------------- 52 !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 53 !! $Id$ 54 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 55 !!---------------------------------------------------------------------- 56 64 !! NEMO/OPA 3.3 , NEMO-Consortium (2010) 65 !! $Id $ 66 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 67 !!---------------------------------------------------------------------- 57 68 CONTAINS 58 69 … … 81 92 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 82 93 !!---------------------------------------------------------------------- 83 USE oce, ONLY : ztrdt => ua ! use ua as 3D workspace84 USE oce, ONLY : ztrds => va ! use va as 3D workspace85 !!86 94 INTEGER, INTENT(in) :: kt ! ocean time-step index 87 95 !! 88 INTEGER :: jk ! dummy loop indices 89 REAL(wp) :: zfact ! temporary scalars 96 INTEGER :: jk, jn ! dummy loop indices 97 REAL(wp) :: zfact ! local scalars 98 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds 90 99 !!---------------------------------------------------------------------- 91 100 … … 94 103 IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap' 95 104 IF(lwp) WRITE(numout,*) '~~~~~~~' 105 ! 106 rbcp = 0.25 * (1. + atfp) * (1. + atfp) * ( 1. - atfp) ! Brown & Campana parameter for semi-implicit hpg 96 107 ENDIF 97 108 98 109 ! Update after tracer on domain lateral boundaries 99 110 ! 100 CALL lbc_lnk( ta, 'T', 1. ) ! local domain boundaries (T-point, unchanged sign) 101 CALL lbc_lnk( sa, 'T', 1. ) 102 ! 103 #if defined key_obc 111 CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) ! local domain boundaries (T-point, unchanged sign) 112 CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) 113 ! 114 #if defined key_obc || defined key_bdy || defined key_agrif 115 CALL tra_unswap 116 #endif 117 118 #if defined key_obc 104 119 IF( lk_obc ) CALL obc_tra( kt ) ! OBC open boundaries 105 120 #endif 106 #if defined key_bdy 107 CALL bdy_tra( kt )! BDY open boundaries121 #if defined key_bdy 122 IF( lk_bdy ) CALL bdy_tra_frs( kt ) ! BDY open boundaries 108 123 #endif 109 124 #if defined key_agrif 110 CALL Agrif_tra ! AGRIF zoom boundaries 125 CALL Agrif_tra ! AGRIF zoom boundaries 126 #endif 127 128 #if defined key_obc || defined key_bdy || defined key_agrif 129 CALL tra_swap 111 130 #endif 112 131 113 132 ! set time step size (Euler/Leapfrog) 114 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt _t(:) = rdttra(:) ! at nit000 (Euler)115 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt _t(:) = 2.* rdttra(:) ! at nit000 or nit000+1 (Leapfrog)133 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt(:) = rdttra(:) ! at nit000 (Euler) 134 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt(:) = 2.* rdttra(:) ! at nit000 or nit000+1 (Leapfrog) 116 135 ENDIF 117 136 118 137 ! trends computation initialisation 119 IF( l_trdtra ) THEN ! store now fields before applying the Asselin filter 120 ztrdt(:,:,:) = tn(:,:,:) 121 ztrds(:,:,:) = sn(:,:,:) 122 ENDIF 123 124 ! Leap-Frog + Asselin filter time stepping 125 IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt ) ! variable volume level (vvl) 126 ELSE ; CALL tra_nxt_fix( kt ) ! fixed volume level 127 ENDIF 138 IF( l_trdtra ) THEN ! store now fields before applying the Asselin filter 139 ALLOCATE( ztrdt(jpi,jpj,jpk) ) ; ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 140 ALLOCATE( ztrds(jpi,jpj,jpk) ) ; ztrds(:,:,:) = tsn(:,:,:,jp_sal) 141 ENDIF 142 143 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step (only swap) 144 DO jn = 1, jpts 145 DO jk = 1, jpkm1 146 tsn(:,:,jk,jn) = tsa(:,:,jk,jn) 147 END DO 148 END DO 149 ELSE ! Leap-Frog + Asselin filter time stepping 150 ! 151 IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt, 'TRA', tsb, tsn, tsa, jpts ) ! variable volume level (vvl) 152 ELSE ; CALL tra_nxt_fix( kt, 'TRA', tsb, tsn, tsa, jpts ) ! fixed volume level 153 ENDIF 154 ENDIF 128 155 129 156 #if defined key_agrif 130 157 ! Update tracer at AGRIF zoom boundaries 158 CALL tra_unswap 131 159 IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tra( kt ) ! children only 160 CALL tra_swap 132 161 #endif 133 162 134 163 ! trends computation 135 IF( l_trdtra ) THEN 164 IF( l_trdtra ) THEN ! trend of the Asselin filter (tb filtered - tb)/dt 136 165 DO jk = 1, jpkm1 137 zfact = 1.e0 / r2dt _t(jk)138 ztrdt(:,:,jk) = ( t b(:,:,jk) - ztrdt(:,:,jk) ) * zfact139 ztrds(:,:,jk) = ( sb(:,:,jk) - ztrds(:,:,jk) ) * zfact166 zfact = 1.e0 / r2dt(jk) 167 ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact 168 ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact 140 169 END DO 141 CALL trd_mod( ztrdt, ztrds, jptra_trd_atf, 'TRA', kt ) 170 CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_atf, ztrdt ) 171 CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_atf, ztrds ) 172 DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds ) 142 173 END IF 143 174 144 175 ! ! control print 145 IF(ln_ctl) CALL prt_ctl( tab3d_1=t n, clinfo1=' nxt - Tn: ', mask1=tmask, &146 & tab3d_2= sn, clinfo2= ' Sn: ', mask2=tmask )176 IF(ln_ctl) CALL prt_ctl( tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' nxt - Tn: ', mask1=tmask, & 177 & tab3d_2=tsn(:,:,:,jp_sal), clinfo2= ' Sn: ', mask2=tmask ) 147 178 ! 148 179 END SUBROUTINE tra_nxt 149 180 150 181 151 SUBROUTINE tra_nxt_fix( kt )182 SUBROUTINE tra_nxt_fix( kt, cdtype, ptb, ptn, pta, kjpt ) 152 183 !!---------------------------------------------------------------------- 153 184 !! *** ROUTINE tra_nxt_fix *** … … 162 193 !! - swap tracer fields to prepare the next time_step. 163 194 !! This can be summurized for tempearture as: 164 !! ztm = (ta+2tn+tb)/4 ln_dynhpg_imp = T 165 !! ztm = 0 otherwise 195 !! ztm = tn + rbcp * [ta -2 tn + tb ] ln_dynhpg_imp = T 196 !! ztm = 0 otherwise 197 !! with rbcp=1/4 * (1-atfp^4) / (1-atfp) 166 198 !! tb = tn + atfp*[ tb - 2 tn + ta ] 167 !! tn = ta 199 !! tn = ta 168 200 !! ta = ztm (NB: reset to 0 after eos_bn2 call) 169 201 !! … … 171 203 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 172 204 !!---------------------------------------------------------------------- 173 INTEGER, INTENT(in) :: kt ! ocean time-step index 174 !! 175 INTEGER :: ji, jj, jk ! dummy loop indices 176 REAL(wp) :: ztm, ztf ! temporary scalars 177 REAL(wp) :: zsm, zsf ! - - 178 !!---------------------------------------------------------------------- 179 180 IF( kt == nit000 ) THEN 205 INTEGER , INTENT(in ) :: kt ! ocean time-step index 206 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 207 INTEGER , INTENT(in ) :: kjpt ! number of tracers 208 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields 209 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields 210 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend 211 !! 212 INTEGER :: ji, jj, jk, jn ! dummy loop indices 213 LOGICAL :: ll_tra_hpg ! local logical 214 REAL(wp) :: ztn, ztd ! local scalars 215 !!---------------------------------------------------------------------- 216 217 IF( kt == nit000 ) THEN 181 218 IF(lwp) WRITE(numout,*) 182 219 IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping' … … 184 221 ENDIF 185 222 ! 186 ! ! ----------------------- ! 187 IF( ln_dynhpg_imp ) THEN ! semi-implicite hpg case ! 188 ! ! ----------------------- ! 223 IF( cdtype == 'TRA' ) THEN ; ll_tra_hpg = ln_dynhpg_imp ! active tracers case and semi-implicit hpg 224 ELSE ; ll_tra_hpg = .FALSE. ! passive tracers case or NO semi-implicit hpg 225 ENDIF 226 ! 227 DO jn = 1, kjpt 189 228 ! 190 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 191 DO jk = 1, jpkm1 ! (only swap) 192 DO jj = 1, jpj 193 DO ji = 1, jpi 194 tn(ji,jj,jk) = ta(ji,jj,jk) ! tb <-- tn 195 sn(ji,jj,jk) = sa(ji,jj,jk) 196 END DO 229 DO jk = 1, jpkm1 230 DO jj = 1, jpj 231 DO ji = 1, jpi 232 ztn = ptn(ji,jj,jk,jn) 233 ztd = pta(ji,jj,jk,jn) - 2. * ztn + ptb(ji,jj,jk,jn) ! time laplacian on tracers 234 ! 235 ptb(ji,jj,jk,jn) = ztn + atfp * ztd ! ptb <-- filtered ptn 236 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 237 ! 238 IF( ll_tra_hpg ) pta(ji,jj,jk,jn) = ztn + rbcp * ztd ! pta <-- Brown & Campana average 197 239 END DO 198 END DO 199 ELSE ! general case (Leapfrog + Asselin filter 200 DO jk = 1, jpkm1 201 DO jj = 1, jpj 202 DO ji = 1, jpi 203 ztm = 0.25 * ( ta(ji,jj,jk) + 2.* tn(ji,jj,jk) + tb(ji,jj,jk) ) ! mean t 204 zsm = 0.25 * ( sa(ji,jj,jk) + 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 205 ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) ) ! Asselin filter on t 206 zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 207 tb(ji,jj,jk) = tn(ji,jj,jk) + ztf ! tb <-- filtered tn 208 sb(ji,jj,jk) = sn(ji,jj,jk) + zsf 209 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 210 sn(ji,jj,jk) = sa(ji,jj,jk) 211 ta(ji,jj,jk) = ztm ! ta <-- mean t 212 sa(ji,jj,jk) = zsm 213 END DO 214 END DO 215 END DO 216 ENDIF 217 ! ! ----------------------- ! 218 ELSE ! explicit hpg case ! 219 ! ! ----------------------- ! 240 END DO 241 END DO 220 242 ! 221 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 222 DO jk = 1, jpkm1 223 DO jj = 1, jpj 224 DO ji = 1, jpi 225 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 226 sn(ji,jj,jk) = sa(ji,jj,jk) 227 END DO 228 END DO 229 END DO 230 ELSE ! general case (Leapfrog + Asselin filter 231 DO jk = 1, jpkm1 232 DO jj = 1, jpj 233 DO ji = 1, jpi 234 ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) ) ! Asselin filter on t 235 zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) ) 236 tb(ji,jj,jk) = tn(ji,jj,jk) + ztf ! tb <-- filtered tn 237 sb(ji,jj,jk) = sn(ji,jj,jk) + zsf 238 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 239 sn(ji,jj,jk) = sa(ji,jj,jk) 240 END DO 241 END DO 242 END DO 243 ENDIF 244 ENDIF 243 END DO 245 244 ! 246 245 END SUBROUTINE tra_nxt_fix 247 246 248 247 249 SUBROUTINE tra_nxt_vvl( kt )248 SUBROUTINE tra_nxt_vvl( kt, cdtype, ptb, ptn, pta, kjpt ) 250 249 !!---------------------------------------------------------------------- 251 250 !! *** ROUTINE tra_nxt_vvl *** … … 260 259 !! - swap tracer fields to prepare the next time_step. 261 260 !! This can be summurized for tempearture as: 262 !! ztm = ( e3t_a*ta+2*e3t_n*tn+e3t_b*tb) ln_dynhpg_imp = T263 !! /( e3t_a +2*e3t_n +e3t_b )264 !! ztm = 0 otherwise261 !! ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) ln_dynhpg_imp = T 262 !! /( e3t_n + rbcp*[ e3t_b - 2 e3t_n + e3t_a ] ) 263 !! ztm = 0 otherwise 265 264 !! tb = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) 266 265 !! /( e3t_n + atfp*[ e3t_b - 2 e3t_n + e3t_a ] ) … … 271 270 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 272 271 !!---------------------------------------------------------------------- 273 INTEGER, INTENT(in) :: kt ! ocean time-step index 272 INTEGER , INTENT(in ) :: kt ! ocean time-step index 273 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 274 INTEGER , INTENT(in ) :: kjpt ! number of tracers 275 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields 276 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields 277 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend 274 278 !! 275 INTEGER :: ji, jj, jk ! dummy loop indices 276 REAL(wp) :: ztm , ztc_f , ztf , ztca, ztcn, ztcb ! temporary scalar 277 REAL(wp) :: zsm , zsc_f , zsf , zsca, zscn, zscb ! - - 278 REAL(wp) :: ze3mr, ze3fr ! - - 279 REAL(wp) :: ze3t_b, ze3t_n, ze3t_a, ze3t_f ! - - 279 LOGICAL :: ll_tra, ll_tra_hpg, ll_traqsr ! local logical 280 INTEGER :: ji, jj, jk, jn ! dummy loop indices 281 REAL(wp) :: ztc_a , ztc_n , ztc_b ! local scalar 282 REAL(wp) :: ztc_f , ztc_d ! - - 283 REAL(wp) :: ze3t_b, ze3t_n, ze3t_a ! - - 284 REAL(wp) :: ze3t_f, ze3t_d ! - - 285 REAL(wp) :: zfact1, zfact2 ! - - 280 286 !!---------------------------------------------------------------------- 281 287 … … 285 291 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 286 292 ENDIF 287 288 ! ! ----------------------- ! 289 IF( ln_dynhpg_imp ) THEN ! semi-implicite hpg case ! 290 ! ! ----------------------- ! 291 ! 292 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step 293 DO jk = 1, jpkm1 ! (only swap) 294 DO jj = 1, jpj 295 DO ji = 1, jpi 296 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 297 sn(ji,jj,jk) = sa(ji,jj,jk) 298 END DO 293 ! 294 IF( cdtype == 'TRA' ) THEN 295 ll_tra = .TRUE. ! active tracers case 296 ll_tra_hpg = ln_dynhpg_imp ! active tracers case and semi-implicit hpg 297 ll_traqsr = ln_traqsr ! active tracers case and solar penetration 298 ELSE 299 ll_tra = .FALSE. ! passive tracers case 300 ll_tra_hpg = .FALSE. ! passive tracers case or NO semi-implicit hpg 301 ll_traqsr = .FALSE. ! active tracers case and NO solar penetration 302 ENDIF 303 ! 304 DO jn = 1, kjpt 305 DO jk = 1, jpkm1 306 zfact1 = atfp * rdttra(jk) 307 zfact2 = zfact1 / rau0 308 DO jj = 1, jpj 309 DO ji = 1, jpi 310 ze3t_b = fse3t_b(ji,jj,jk) 311 ze3t_n = fse3t_n(ji,jj,jk) 312 ze3t_a = fse3t_a(ji,jj,jk) 313 ! ! tracer content at Before, now and after 314 ztc_b = ptb(ji,jj,jk,jn) * ze3t_b 315 ztc_n = ptn(ji,jj,jk,jn) * ze3t_n 316 ztc_a = pta(ji,jj,jk,jn) * ze3t_a 317 ! 318 ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b 319 ztc_d = ztc_a - 2. * ztc_n + ztc_b 320 ! 321 ze3t_f = ze3t_n + atfp * ze3t_d 322 ztc_f = ztc_n + atfp * ztc_d 323 ! 324 IF( ll_tra .AND. jk == 1 ) THEN ! first level only for T & S 325 ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) ) 326 ztc_f = ztc_f - zfact1 * ( sbc_tsc(ji,jj,jn) - sbc_tsc_b(ji,jj,jn) ) 327 ENDIF 328 IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr ) & ! solar penetration (temperature only) 329 & ztc_f = ztc_f - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 330 331 ze3t_f = 1.e0 / ze3t_f 332 ptb(ji,jj,jk,jn) = ztc_f * ze3t_f ! ptb <-- ptn filtered 333 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 334 ! 335 IF( ll_tra_hpg ) THEN ! semi-implicit hpg (T & S only) 336 ze3t_d = 1.e0 / ( ze3t_n + rbcp * ze3t_d ) 337 pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n + rbcp * ztc_d ) ! ta <-- Brown & Campana average 338 ENDIF 299 339 END DO 300 340 END DO 301 ELSE 302 DO jk = 1, jpkm1 303 DO jj = 1, jpj 304 DO ji = 1, jpi 305 ze3t_b = fse3t_b(ji,jj,jk) 306 ze3t_n = fse3t_n(ji,jj,jk) 307 ze3t_a = fse3t_a(ji,jj,jk) 308 ! ! tracer content at Before, now and after 309 ztcb = tb(ji,jj,jk) * ze3t_b ; zscb = sb(ji,jj,jk) * ze3t_b 310 ztcn = tn(ji,jj,jk) * ze3t_n ; zscn = sn(ji,jj,jk) * ze3t_n 311 ztca = ta(ji,jj,jk) * ze3t_a ; zsca = sa(ji,jj,jk) * ze3t_a 312 ! 313 ! ! Asselin filter on thickness and tracer content 314 ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b ) 315 ztc_f = atfp * ( ztca - 2.* ztcn + ztcb ) 316 zsc_f = atfp * ( zsca - 2.* zscn + zscb ) 317 ! 318 ! ! filtered tracer including the correction 319 ze3fr = 1.e0 / ( ze3t_n + ze3t_f ) 320 ztf = ze3fr * ( ztcn + ztc_f ) 321 zsf = ze3fr * ( zscn + zsc_f ) 322 ! ! mean thickness and tracer 323 ze3mr = 1.e0 / ( ze3t_a + 2.* ze3t_n + ze3t_b ) 324 ztm = ze3mr * ( ztca + 2.* ztcn + ztcb ) 325 zsm = ze3mr * ( zsca + 2.* zscn + zscb ) 326 !!gm mean e3t have to be saved and used in dynhpg or it can be recomputed in dynhpg !! 327 !!gm e3t_m(ji,jj,jk) = 0.25 / ze3mr 328 ! ! swap of arrays 329 tb(ji,jj,jk) = ztf ! tb <-- tn + filter 330 sb(ji,jj,jk) = zsf 331 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 332 sn(ji,jj,jk) = sa(ji,jj,jk) 333 ta(ji,jj,jk) = ztm ! ta <-- mean t 334 sa(ji,jj,jk) = zsm 335 END DO 336 END DO 337 END DO 338 ENDIF 339 ! ! ----------------------- ! 340 ELSE ! explicit hpg case ! 341 ! ! ----------------------- ! 342 ! 343 IF( neuler == 0 .AND. kt == nit000 ) THEN ! case of Euler time-stepping at first time-step 344 DO jk = 1, jpkm1 ! No filter nor thickness weighting computation required 345 DO jj = 1, jpj ! ONLY swap 346 DO ji = 1, jpi 347 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 348 sn(ji,jj,jk) = sa(ji,jj,jk) 349 END DO 350 END DO 351 END DO 352 ! ! general case (Leapfrog + Asselin filter) 353 ELSE ! apply filter on thickness weighted tracer and swap 354 DO jk = 1, jpkm1 355 DO jj = 1, jpj 356 DO ji = 1, jpi 357 ze3t_b = fse3t_b(ji,jj,jk) 358 ze3t_n = fse3t_n(ji,jj,jk) 359 ze3t_a = fse3t_a(ji,jj,jk) 360 ! ! tracer content at Before, now and after 361 ztcb = tb(ji,jj,jk) * ze3t_b ; zscb = sb(ji,jj,jk) * ze3t_b 362 ztcn = tn(ji,jj,jk) * ze3t_n ; zscn = sn(ji,jj,jk) * ze3t_n 363 ztca = ta(ji,jj,jk) * ze3t_a ; zsca = sa(ji,jj,jk) * ze3t_a 364 ! 365 ! ! Asselin filter on thickness and tracer content 366 ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b ) 367 ztc_f = atfp * ( ztca - 2.* ztcn + ztcb ) 368 zsc_f = atfp * ( zsca - 2.* zscn + zscb ) 369 ! 370 ! ! filtered tracer including the correction 371 ze3fr = 1.e0 / ( ze3t_n + ze3t_f ) 372 ztf = ( ztcn + ztc_f ) * ze3fr 373 zsf = ( zscn + zsc_f ) * ze3fr 374 ! ! swap of arrays 375 tb(ji,jj,jk) = ztf ! tb <-- tn filtered 376 sb(ji,jj,jk) = zsf 377 tn(ji,jj,jk) = ta(ji,jj,jk) ! tn <-- ta 378 sn(ji,jj,jk) = sa(ji,jj,jk) 379 END DO 380 END DO 381 END DO 382 ENDIF 383 ENDIF 341 END DO 342 ! 343 END DO 384 344 ! 385 345 END SUBROUTINE tra_nxt_vvl
Note: See TracChangeset
for help on using the changeset viewer.
