- Timestamp:
- 2016-07-19T10:38:35+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
r5467 r6808 28 28 USE sbc_oce ! surface boundary condition: ocean 29 29 USE sbcrnf ! river runoffs 30 USE sbcisf ! ice shelf melting 30 31 USE zdf_oce ! ocean vertical mixing 31 32 USE domvvl ! variable volume 32 USE dynspg_oce ! surface pressure gradient variables33 USE dynhpg ! hydrostatic pressure gradient34 33 USE trd_oce ! trends: ocean variables 35 34 USE trdtra ! trends manager: tracers 36 35 USE traqsr ! penetrative solar radiation (needed for nksr) 37 36 USE phycst ! physical constant 38 USE ldftra_oce ! lateral physics on tracers 37 USE ldftra ! lateral physics on tracers 38 USE ldfslp 39 39 USE bdy_oce ! BDY open boundary condition variables 40 40 USE bdytra ! open boundary condition (bdy_tra routine) … … 46 46 USE timing ! Timing 47 47 #if defined key_agrif 48 USE agrif_opa_update49 48 USE agrif_opa_interp 50 49 #endif … … 57 56 PUBLIC tra_nxt_vvl ! to be used in trcnxt 58 57 59 REAL(wp) :: rbcp ! Brown & Campana parameters for semi-implicit hpg60 61 58 !! * Substitutions 62 # include " domzgr_substitute.h90"59 # include "vectopt_loop_substitute.h90" 63 60 !!---------------------------------------------------------------------- 64 61 !! NEMO/OPA 3.3 , NEMO-Consortium (2010) … … 88 85 !! domains (lk_agrif=T) 89 86 !! 90 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step 91 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 87 !! ** Action : - tsb & tsn ready for the next time step 92 88 !!---------------------------------------------------------------------- 93 89 INTEGER, INTENT(in) :: kt ! ocean time-step index 94 90 !! 95 INTEGER :: j k, jn! dummy loop indices96 REAL(wp) :: zfact ! local scalars91 INTEGER :: ji, jj, jk, jn ! dummy loop indices 92 REAL(wp) :: zfact ! local scalars 97 93 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt, ztrds 98 94 !!---------------------------------------------------------------------- … … 104 100 IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap' 105 101 IF(lwp) WRITE(numout,*) '~~~~~~~' 106 !107 rbcp = 0.25_wp * (1._wp + atfp) * (1._wp + atfp) * ( 1._wp - atfp) ! Brown & Campana parameter for semi-implicit hpg108 102 ENDIF 109 103 110 104 ! Update after tracer on domain lateral boundaries 111 105 ! 106 #if defined key_agrif 107 CALL Agrif_tra ! AGRIF zoom boundaries 108 #endif 109 ! 112 110 CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1._wp ) ! local domain boundaries (T-point, unchanged sign) 113 111 CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1._wp ) … … 116 114 IF( lk_bdy ) CALL bdy_tra( kt ) ! BDY open boundaries 117 115 #endif 118 #if defined key_agrif119 CALL Agrif_tra ! AGRIF zoom boundaries120 #endif121 116 122 117 ! set time step size (Euler/Leapfrog) 123 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt ra(:) = rdttra(:)! at nit000 (Euler)124 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt ra(:) = 2._wp* rdttra(:)! at nit000 or nit000+1 (Leapfrog)118 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! at nit000 (Euler) 119 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2._wp* rdt ! at nit000 or nit000+1 (Leapfrog) 125 120 ENDIF 126 121 … … 142 137 END DO 143 138 END DO 139 ! 144 140 ELSE ! Leap-Frog + Asselin filter time stepping 145 141 ! 146 IF( l k_vvl ) THEN ; CALL tra_nxt_vvl( kt, nit000, rdttra, 'TRA', tsb, tsn, tsa, &147 & sbc_tsc, sbc_tsc_b, jpts ) ! variable volume level (vvl)148 ELSE ; CALL tra_nxt_fix( kt, nit000, 'TRA', tsb, tsn, tsa, jpts ) ! fixed volume level142 IF( ln_linssh ) THEN ; CALL tra_nxt_fix( kt, nit000, 'TRA', tsb, tsn, tsa, jpts ) ! linear free surface 143 ELSE ; CALL tra_nxt_vvl( kt, nit000, rdt, 'TRA', tsb, tsn, tsa, & 144 & sbc_tsc, sbc_tsc_b, jpts ) ! non-linear free surface 149 145 ENDIF 150 ENDIF151 !152 #if defined key_agrif 153 ! Update tracer at AGRIF zoom boundaries154 IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tra( kt ) ! children only155 #endif 156 !157 ! trends computation146 ! 147 DO jn = 1, jpts 148 CALL lbc_lnk( tsb(:,:,:,jn), 'T', 1._wp ) 149 CALL lbc_lnk( tsn(:,:,:,jn), 'T', 1._wp ) 150 CALL lbc_lnk( tsa(:,:,:,jn), 'T', 1._wp ) 151 END DO 152 ENDIF 153 ! 158 154 IF( l_trdtra ) THEN ! trend of the Asselin filter (tb filtered - tb)/dt 159 155 DO jk = 1, jpkm1 160 zfact = 1._wp / r2dt ra(jk)156 zfact = 1._wp / r2dt 161 157 ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact 162 158 ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact … … 184 180 !! 185 181 !! ** Method : - Apply a Asselin time filter on now fields. 186 !! - save in (ta,sa) an average over the three time levels187 !! which will be used to compute rdn and thus the semi-implicit188 !! hydrostatic pressure gradient (ln_dynhpg_imp = T)189 182 !! - swap tracer fields to prepare the next time_step. 190 !! This can be summurized for tempearture as: 191 !! ztm = tn + rbcp * [ta -2 tn + tb ] ln_dynhpg_imp = T 192 !! ztm = 0 otherwise 193 !! with rbcp=1/4 * (1-atfp^4) / (1-atfp) 194 !! tb = tn + atfp*[ tb - 2 tn + ta ] 195 !! tn = ta 196 !! ta = ztm (NB: reset to 0 after eos_bn2 call) 197 !! 198 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step 199 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 200 !!---------------------------------------------------------------------- 201 INTEGER , INTENT(in ) :: kt ! ocean time-step index 202 INTEGER , INTENT(in ) :: kit000 ! first time step index 203 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 204 INTEGER , INTENT(in ) :: kjpt ! number of tracers 205 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields 206 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields 207 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend 183 !! 184 !! ** Action : - tsb & tsn ready for the next time step 185 !!---------------------------------------------------------------------- 186 INTEGER , INTENT(in ) :: kt ! ocean time-step index 187 INTEGER , INTENT(in ) :: kit000 ! first time step index 188 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 189 INTEGER , INTENT(in ) :: kjpt ! number of tracers 190 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: ptb ! before tracer fields 191 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: ptn ! now tracer fields 192 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 208 193 ! 209 194 INTEGER :: ji, jj, jk, jn ! dummy loop indices 210 LOGICAL :: ll_tra_hpg ! local logical211 195 REAL(wp) :: ztn, ztd ! local scalars 212 196 !!---------------------------------------------------------------------- 213 197 ! 214 198 IF( kt == kit000 ) THEN 215 199 IF(lwp) WRITE(numout,*) … … 218 202 ENDIF 219 203 ! 220 IF( cdtype == 'TRA' ) THEN ; ll_tra_hpg = ln_dynhpg_imp ! active tracers case and semi-implicit hpg221 ELSE ; ll_tra_hpg = .FALSE. ! passive tracers case or NO semi-implicit hpg222 ENDIF223 !224 204 DO jn = 1, kjpt 225 205 ! 226 206 DO jk = 1, jpkm1 227 DO jj = 1, jpj228 DO ji = 1, jpi207 DO jj = 2, jpjm1 208 DO ji = fs_2, fs_jpim1 229 209 ztn = ptn(ji,jj,jk,jn) 230 ztd = pta(ji,jj,jk,jn) - 2. * ztn + ptb(ji,jj,jk,jn) ! time laplacian on tracers 231 ! 232 ptb(ji,jj,jk,jn) = ztn + atfp * ztd ! ptb <-- filtered ptn 233 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 234 ! 235 IF( ll_tra_hpg ) pta(ji,jj,jk,jn) = ztn + rbcp * ztd ! pta <-- Brown & Campana average 210 ztd = pta(ji,jj,jk,jn) - 2._wp * ztn + ptb(ji,jj,jk,jn) ! time laplacian on tracers 211 ! 212 ptb(ji,jj,jk,jn) = ztn + atfp * ztd ! ptb <-- filtered ptn 213 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 236 214 END DO 237 215 END DO … … 251 229 !! 252 230 !! ** Method : - Apply a thickness weighted Asselin time filter on now fields. 253 !! - save in (ta,sa) a thickness weighted average over the three254 !! time levels which will be used to compute rdn and thus the semi-255 !! implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)256 231 !! - swap tracer fields to prepare the next time_step. 257 !! This can be summurized for tempearture as:258 !! ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) ln_dynhpg_imp = T259 !! /( e3t_n + rbcp*[ e3t_b - 2 e3t_n + e3t_a ] )260 !! ztm = 0 otherwise261 232 !! tb = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) 262 233 !! /( e3t_n + atfp*[ e3t_b - 2 e3t_n + e3t_a ] ) 263 234 !! tn = ta 264 !! ta = zt (NB: reset to 0 after eos_bn2 call) 265 !! 266 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step 267 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 268 !!---------------------------------------------------------------------- 269 INTEGER , INTENT(in ) :: kt ! ocean time-step index 270 INTEGER , INTENT(in ) :: kit000 ! first time step index 271 REAL(wp) , INTENT(in ), DIMENSION(jpk) :: p2dt ! time-step 272 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 273 INTEGER , INTENT(in ) :: kjpt ! number of tracers 274 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields 275 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields 276 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend 277 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,kjpt) :: psbc_tc ! surface tracer content 278 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,kjpt) :: psbc_tc_b ! before surface tracer content 279 280 !! 281 LOGICAL :: ll_tra_hpg, ll_traqsr, ll_rnf ! local logical 235 !! 236 !! ** Action : - tsb & tsn ready for the next time step 237 !!---------------------------------------------------------------------- 238 INTEGER , INTENT(in ) :: kt ! ocean time-step index 239 INTEGER , INTENT(in ) :: kit000 ! first time step index 240 REAL(wp) , INTENT(in ) :: p2dt ! time-step 241 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 242 INTEGER , INTENT(in ) :: kjpt ! number of tracers 243 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: ptb ! before tracer fields 244 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: ptn ! now tracer fields 245 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 246 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: psbc_tc ! surface tracer content 247 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: psbc_tc_b ! before surface tracer content 248 ! 249 LOGICAL :: ll_traqsr, ll_rnf, ll_isf ! local logical 282 250 INTEGER :: ji, jj, jk, jn ! dummy loop indices 283 251 REAL(wp) :: zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d ! local scalar … … 292 260 ! 293 261 IF( cdtype == 'TRA' ) THEN 294 ll_tra_hpg = ln_dynhpg_imp ! active tracers case and semi-implicit hpg295 262 ll_traqsr = ln_traqsr ! active tracers case and solar penetration 296 263 ll_rnf = ln_rnf ! active tracers case and river runoffs 297 ELSE 298 ll_tra_hpg = .FALSE. ! passive tracers case or NO semi-implicit hpg 299 ll_traqsr = .FALSE. ! active tracers case and NO solar penetration 300 ll_rnf = .FALSE. ! passive tracers or NO river runoffs 264 ll_isf = ln_isf ! active tracers case and ice shelf melting 265 ELSE ! passive tracers case 266 ll_traqsr = .FALSE. ! NO solar penetration 267 ll_rnf = .FALSE. ! NO river runoffs ???? !!gm BUG ? 268 ll_isf = .FALSE. ! NO ice shelf melting/freezing !!gm BUG ?? 301 269 ENDIF 302 270 ! 303 271 DO jn = 1, kjpt 304 272 DO jk = 1, jpkm1 305 zfact1 = atfp * p2dt (jk)306 zfact2 = zfact1 /rau0307 DO jj = 1, jpj308 DO ji = 1, jpi309 ze3t_b = fse3t_b(ji,jj,jk)310 ze3t_n = fse3t_n(ji,jj,jk)311 ze3t_a = fse3t_a(ji,jj,jk)273 zfact1 = atfp * p2dt 274 zfact2 = zfact1 * r1_rau0 275 DO jj = 2, jpjm1 276 DO ji = fs_2, fs_jpim1 277 ze3t_b = e3t_b(ji,jj,jk) 278 ze3t_n = e3t_n(ji,jj,jk) 279 ze3t_a = e3t_a(ji,jj,jk) 312 280 ! ! tracer content at Before, now and after 313 281 ztc_b = ptb(ji,jj,jk,jn) * ze3t_b … … 321 289 ztc_f = ztc_n + atfp * ztc_d 322 290 ! 323 IF( jk == 1 ) THEN ! first level 324 ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) + rnf(ji,jj) - rnf_b(ji,jj) ) 291 IF( jk == mikt(ji,jj) ) THEN ! first level 292 ze3t_f = ze3t_f - zfact2 * ( (emp_b(ji,jj) - emp(ji,jj) ) & 293 & - (rnf_b(ji,jj) - rnf(ji,jj) ) & 294 & + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) 325 295 ztc_f = ztc_f - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) ) 326 296 ENDIF 327 328 IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr ) & ! solar penetration (temperature only) 297 ! 298 ! solar penetration (temperature only) 299 IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr ) & 329 300 & ztc_f = ztc_f - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 330 331 IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) ) & ! river runoffs 301 ! 302 ! river runoff 303 IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) ) & 332 304 & ztc_f = ztc_f - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 333 & * fse3t_n(ji,jj,jk) / h_rnf(ji,jj) 334 305 & * e3t_n(ji,jj,jk) / h_rnf(ji,jj) 306 ! 307 ! ice shelf 308 IF( ll_isf ) THEN 309 ! level fully include in the Losch_2008 ice shelf boundary layer 310 IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) ) & 311 ztc_f = ztc_f - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) ) & 312 & * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) 313 ! level partially include in Losch_2008 ice shelf boundary layer 314 IF ( jk == misfkb(ji,jj) ) & 315 ztc_f = ztc_f - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) ) & 316 & * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj) 317 END IF 318 ! 335 319 ze3t_f = 1.e0 / ze3t_f 336 320 ptb(ji,jj,jk,jn) = ztc_f * ze3t_f ! ptb <-- ptn filtered 337 321 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 338 322 ! 339 IF( ll_tra_hpg ) THEN ! semi-implicit hpg (T & S only)340 ze3t_d = 1.e0 / ( ze3t_n + rbcp * ze3t_d )341 pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n + rbcp * ztc_d ) ! ta <-- Brown & Campana average342 ENDIF343 323 END DO 344 324 END DO
Note: See TracChangeset
for help on using the changeset viewer.