Changeset 1110 for trunk/NEMO/OPA_SRC/TRA/tranxt.F90
- Timestamp:
- 2008-06-13T16:10:35+02:00 (16 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/TRA/tranxt.F90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/TRA/tranxt.F90
r911 r1110 4 4 !! Ocean active tracers: time stepping on temperature and salinity 5 5 !!====================================================================== 6 !! History : 7.0 ! 91-11 (G. Madec) Original code 7 !! ! 93-03 (M. Guyon) symetrical conditions 8 !! ! 96-02 (G. Madec & M. Imbard) opa release 8.0 9 !! 8.0 ! 96-04 (A. Weaver) Euler forward step 10 !! 8.2 ! 99-02 (G. Madec, N. Grima) semi-implicit pressure grad. 11 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 12 !! ! 02-11 (C. Talandier, A-M Treguier) Open boundaries 13 !! ! 05-04 (C. Deltel) Add Asselin trend in the ML budget 14 !! 9.0 ! 06-02 (L. Debreu, C. Mazauric) Agrif implementation 6 !! History : OPA ! 1991-11 (G. Madec) Original code 7 !! 7.0 ! 1993-03 (M. Guyon) symetrical conditions 8 !! 8.0 ! 1996-02 (G. Madec & M. Imbard) opa release 8.0 9 !! - ! 1996-04 (A. Weaver) Euler forward step 10 !! 8.2 ! 1999-02 (G. Madec, N. Grima) semi-implicit pressure grad. 11 !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module 12 !! - ! 2002-11 (C. Talandier, A-M Treguier) Open boundaries 13 !! - ! 2005-04 (C. Deltel) Add Asselin trend in the ML budget 14 !! 2.0 ! 2006-02 (L. Debreu, C. Mazauric) Agrif implementation 15 !! 3.0 ! 2008-06 (G. Madec) time stepping always done in trazdf 15 16 !!---------------------------------------------------------------------- 16 17 17 18 !!---------------------------------------------------------------------- 18 !! tra_nxt : time stepping on temperature and salinity19 !! tra_nxt : time stepping on temperature and salinity 19 20 !!---------------------------------------------------------------------- 20 21 USE oce ! ocean dynamics and tracers variables … … 25 26 USE trdmod ! ocean active tracers trends 26 27 USE phycst 27 USE domvvl ! variable volume28 28 USE obctra ! open boundary condition (obc_tra routine) 29 29 USE bdytra ! Unstructured open boundary condition (bdy_tra routine) … … 34 34 USE agrif_opa_interp 35 35 36 37 36 IMPLICIT NONE 38 37 PRIVATE 39 38 40 !! * Routine accessibility 41 PUBLIC tra_nxt ! routine called by step.F90 42 43 REAL(wp) :: vemp ! total amount of volume added or removed by E-P forcing 39 PUBLIC tra_nxt ! routine called by step.F90 44 40 45 41 !! * Substitutions 46 42 # include "domzgr_substitute.h90" 47 43 !!---------------------------------------------------------------------- 48 !! OPA 9.0 , LOCEAN-IPSL (2006)49 !! $Id $44 !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 45 !! $Id:$ 50 46 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 51 47 !!---------------------------------------------------------------------- … … 57 53 !! *** ROUTINE tranxt *** 58 54 !! 59 !! ** Purpose : Compute the temperature and salinity fields at the 60 !! next time-step from their temporal trends and swap the fields. 55 !! ** Purpose : Apply the boundary condition on the after temperature 56 !! and salinity fields, achieved the time stepping by adding 57 !! the Asselin filter on now fields and swapping the fields. 61 58 !! 62 !! ** Method : Apply lateral boundary conditions on (ua,va) through 63 !! call to lbc_lnk routine 64 !! After t and s are compute using a leap-frog scheme environment: 65 !! ta = tb + 2 rdttra(k) * ta 66 !! sa = sb + 2 rdttra(k) * sa 67 !! Compute and save in (ta,sa) an average over three time levels 68 !! (before,now and after) of temperature and salinity which is 69 !! used to compute rhd in eos routine and thus the hydrostatic 70 !! pressure gradient (ln_dynhpg_imp = T) 71 !! Apply an Asselin time filter on now tracers (tn,sn) to avoid 72 !! the divergence of two consecutive time-steps and swap tracer 73 !! arrays to prepare the next time_step: 74 !! (zt,zs) = (ta+2tn+tb,sa+2sn+sb)/4 (ln_dynhpg_imp = T) 75 !! (zt,zs) = (0,0) (default option) 76 !! (tb,sb) = (tn,vn) + atfp [ (tb,sb) + (ta,sa) - 2 (tn,sn) ] 77 !! (tn,sn) = (ta,sa) 78 !! (ta,sa) = (zt,zs) (NB: reset to 0 after use in eos.F) 59 !! ** Method : At this stage of the computation, ta and sa are the 60 !! after temperature and salinity as the time stepping has 61 !! been performed in trazdf_imp or trazdf_exp module. 79 62 !! 80 !! ** Action : - update (tb,sb) and (tn,sn) 81 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 63 !! - Apply lateral boundary conditions on (ta,sa) 64 !! at the local domain boundaries through lbc_lnk call, 65 !! at the radiative open boundaries (lk_obc=T), 66 !! at the relaxed open boundaries (lk_bdy=T), and 67 !! at the AGRIF zoom boundaries (lk_agrif=T) 68 !! 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) 80 !! 81 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step 82 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 82 83 !!---------------------------------------------------------------------- 83 84 USE oce, ONLY : ztrdt => ua ! use ua as 3D workspace … … 87 88 !! 88 89 INTEGER :: ji, jj, jk ! dummy loop indices 89 REAL(wp) :: zt, zs, zssh1 ! temporary scalars 90 REAL(wp) :: zfact ! temporary scalar 91 !! Variable volume 92 REAL(wp), DIMENSION(jpi,jpj) :: zssh ! temporary scalars 93 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfse3tb, zfse3tn, zfse3ta ! 3D workspace 94 90 REAL(wp) :: zt, zs, zfact ! temporary scalars 95 91 !!---------------------------------------------------------------------- 96 92 97 !! Explicit physics with thickness weighted updates 98 IF( lk_vvl .AND. ln_zdfexp ) THEN 99 100 ! Scale factors at before and after time step 101 ! ------------------------------------------- 102 CALL dom_vvl_sf( sshb, 'T', zfse3tb ) ; CALL dom_vvl_sf( ssha, 'T', zfse3ta ) 103 104 ! Asselin filtered scale factor at now time step 105 ! ---------------------------------------------- 106 IF( (neuler == 0 .AND. kt == nit000) .OR. lk_dynspg_ts ) THEN 107 CALL dom_vvl_sf_ini( 'T', zfse3tn ) 108 ELSE 109 zssh(:,:) = atfp * ( sshb(:,:) + ssha(:,:) ) + atfp1 * sshn(:,:) 110 CALL dom_vvl_sf( zssh, 'T', zfse3tn ) 111 ENDIF 112 113 ! Thickness weighting 114 ! ------------------- 115 DO jk = 1, jpkm1 116 DO jj = 1, jpj 117 DO ji = 1, jpi 118 ta(ji,jj,jk) = ta(ji,jj,jk) * fse3t(ji,jj,jk) 119 sa(ji,jj,jk) = sa(ji,jj,jk) * fse3t(ji,jj,jk) 120 121 tn(ji,jj,jk) = tn(ji,jj,jk) * fse3t(ji,jj,jk) 122 sn(ji,jj,jk) = sn(ji,jj,jk) * fse3t(ji,jj,jk) 123 124 tb(ji,jj,jk) = tb(ji,jj,jk) * zfse3tb(ji,jj,jk) 125 sb(ji,jj,jk) = sb(ji,jj,jk) * zfse3tb(ji,jj,jk) 126 END DO 127 END DO 128 END DO 129 93 IF( kt == nit000 ) THEN 94 IF(lwp) WRITE(numout,*) 95 IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap' 96 IF(lwp) WRITE(numout,*) '~~~~~~~' 130 97 ENDIF 131 98 132 IF( l_trdtra ) THEN 133 ztrdt(:,:,jpk) = 0.e0 134 ztrds(:,:,jpk) = 0.e0 99 ! Update after tracer on domain lateral boundaries 100 ! 101 CALL lbc_lnk( ta, 'T', 1. ) ! local domain boundaries (T-point, unchanged sign) 102 CALL lbc_lnk( sa, 'T', 1. ) 103 ! 104 #if defined key_obc 105 CALL obc_tra( kt ) ! OBC open boundaries 106 #endif 107 #if defined key_bdy 108 CALL bdy_tra( kt ) ! BDY open boundaries 109 #endif 110 #if defined key_agrif 111 CALL Agrif_tra ! AGRIF zoom boundaries 112 #endif 113 114 ! trends computation initialisation 115 IF( l_trdtra ) THEN ! store now fields before applying the Asselin filter 116 ztrdt(:,:,:) = tn(:,:,:) 117 ztrds(:,:,:) = sn(:,:,:) 135 118 ENDIF 136 119 137 ! 0. Lateral boundary conditions on ( ta, sa ) (T-point, unchanged sign) 138 ! ---------------------------------============ 139 CALL lbc_lnk( ta, 'T', 1. ) 140 CALL lbc_lnk( sa, 'T', 1. ) 141 142 ! ! =============== 143 DO jk = 1, jpkm1 ! Horizontal slab 144 ! ! =============== 145 146 ! 1. Leap-frog scheme (only in explicit case, otherwise the 147 ! ------------------- time stepping is already done in trazdf) 148 IF( ln_zdfexp ) THEN 149 zfact = 2. * rdttra(jk) 150 IF( neuler == 0 .AND. kt == nit000 ) zfact = rdttra(jk) 151 ta(:,:,jk) = ( tb(:,:,jk) + zfact * ta(:,:,jk) ) * tmask(:,:,jk) 152 sa(:,:,jk) = ( sb(:,:,jk) + zfact * sa(:,:,jk) ) * tmask(:,:,jk) 153 IF(l_trdtra) CALL ctl_stop( 'tranxt: Asselin ML trend not yet accounted for.' ) 154 ENDIF 155 156 #if defined key_obc 157 ! ! =============== 158 END DO ! End of slab 159 ! ! =============== 160 ! Update tracers on open boundaries. 161 CALL obc_tra( kt ) 162 163 ! ! =============== 164 DO jk = 1, jpkm1 ! Horizontal slab 165 ! ! =============== 166 #elif defined key_bdy 167 ! ! =============== 168 END DO ! End of slab 169 ! ! =============== 170 171 ! Update tracers on open boundaries. 172 CALL bdy_tra( kt ) 173 174 ! ! =============== 175 DO jk = 1, jpkm1 ! Horizontal slab 176 ! ! =============== 177 #endif 178 #if defined key_agrif 179 ! ! =============== 180 END DO ! End of slab 181 ! ! =============== 182 ! Update tracers on open boundaries. 183 CALL Agrif_tra 184 ! ! =============== 185 DO jk = 1, jpkm1 ! Horizontal slab 186 ! ! =============== 187 #endif 188 ! 2. Time filter and swap of arrays 189 ! --------------------------------- 190 191 IF( ln_dynhpg_imp ) THEN ! semi-implicite hpg 192 IF( neuler == 0 .AND. kt == nit000 ) THEN 193 DO jj = 1, jpj 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) 194 135 DO ji = 1, jpi 195 136 zt = ( ta(ji,jj,jk) + 2. * tn(ji,jj,jk) + tb(ji,jj,jk) ) * 0.25 196 137 zs = ( sa(ji,jj,jk) + 2. * sn(ji,jj,jk) + sb(ji,jj,jk) ) * 0.25 197 tb(ji,jj,jk) = tn(ji,jj,jk)198 sb(ji,jj,jk) = sn(ji,jj,jk)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) 199 140 tn(ji,jj,jk) = ta(ji,jj,jk) 200 141 sn(ji,jj,jk) = sa(ji,jj,jk) … … 203 144 END DO 204 145 END DO 205 IF( l_trdtra ) THEN 206 ztrdt(:,:,jk) = 0.e0 207 ztrds(:,:,jk) = 0.e0 208 END IF 209 ELSE 146 END DO 147 ELSE ! explicit hpg case 148 DO jk = 1, jpkm1 210 149 DO jj = 1, jpj 211 150 DO ji = 1, jpi 212 zt = ( ta(ji,jj,jk) + 2. * tn(ji,jj,jk) + tb(ji,jj,jk) ) * 0.25 213 zs = ( sa(ji,jj,jk) + 2. * sn(ji,jj,jk) + sb(ji,jj,jk) ) * 0.25 214 tb(ji,jj,jk) = atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk) 215 sb(ji,jj,jk) = atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk) 216 IF( l_trdtra ) THEN ! ChD ceci est a optimiser, mais ca marche 217 ztrdt(ji,jj,jk) = tb(ji,jj,jk) - tn(ji,jj,jk) 218 ztrds(ji,jj,jk) = sb(ji,jj,jk) - sn(ji,jj,jk) 219 END IF 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) 220 153 tn(ji,jj,jk) = ta(ji,jj,jk) 221 154 sn(ji,jj,jk) = sa(ji,jj,jk) 222 ta(ji,jj,jk) = zt223 sa(ji,jj,jk) = zs224 155 END DO 225 156 END DO 226 ENDIF 227 ELSE ! Default case 228 IF( neuler == 0 .AND. kt == nit000 ) THEN 229 IF( (lk_vvl .AND. ln_zdfexp) ) THEN ! Varying levels 230 DO jj = 1, jpj 231 DO ji = 1, jpi 232 zssh1 = tmask(ji,jj,jk) / fse3t(ji,jj,jk) 233 tb(ji,jj,jk) = tn(ji,jj,jk) * zssh1 * tmask(ji,jj,jk) 234 sb(ji,jj,jk) = sn(ji,jj,jk) * zssh1 * tmask(ji,jj,jk) 235 zssh1 = tmask(ji,jj,jk) / zfse3ta(ji,jj,jk) 236 tn(ji,jj,jk) = ta(ji,jj,jk) * zssh1 * tmask(ji,jj,jk) 237 sn(ji,jj,jk) = sa(ji,jj,jk) * zssh1 * tmask(ji,jj,jk) 238 END DO 239 END DO 240 ELSE ! Fixed levels 241 DO jj = 1, jpj 242 DO ji = 1, jpi 243 tb(ji,jj,jk) = tn(ji,jj,jk) 244 sb(ji,jj,jk) = sn(ji,jj,jk) 245 tn(ji,jj,jk) = ta(ji,jj,jk) 246 sn(ji,jj,jk) = sa(ji,jj,jk) 247 END DO 248 END DO 249 ENDIF 250 IF( l_trdtra ) THEN 251 ztrdt(:,:,jk) = 0.e0 252 ztrds(:,:,jk) = 0.e0 253 END IF 254 ELSE 255 IF( l_trdtra ) THEN 256 DO jj = 1, jpj 257 DO ji = 1, jpi 258 ztrdt(ji,jj,jk) = atfp * ( tb(ji,jj,jk) - 2*tn(ji,jj,jk) + ta(ji,jj,jk) ) 259 ztrds(ji,jj,jk) = atfp * ( sb(ji,jj,jk) - 2*sn(ji,jj,jk) + sa(ji,jj,jk) ) 260 END DO 261 END DO 262 END IF 263 IF( (lk_vvl .AND. ln_zdfexp) ) THEN ! Varying levels 264 DO jj = 1, jpj 265 DO ji = 1, jpi 266 zssh1 = tmask(ji,jj,jk) / zfse3tn(ji,jj,jk) 267 tb(ji,jj,jk) = ( atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) & 268 & + atfp1 * tn(ji,jj,jk) ) * zssh1 269 sb(ji,jj,jk) = ( atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) & 270 & + atfp1 * sn(ji,jj,jk) ) * zssh1 271 zssh1 = tmask(ji,jj,1) / zfse3ta(ji,jj,jk) 272 tn(ji,jj,jk) = ta(ji,jj,jk) * zssh1 273 sn(ji,jj,jk) = sa(ji,jj,jk) * zssh1 274 END DO 275 END DO 276 ELSE ! Fixed levels or first varying level 277 DO jj = 1, jpj 278 DO ji = 1, jpi 279 tb(ji,jj,jk) = atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk) 280 sb(ji,jj,jk) = atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk) 281 tn(ji,jj,jk) = ta(ji,jj,jk) 282 sn(ji,jj,jk) = sa(ji,jj,jk) 283 END DO 284 END DO 285 ENDIF 286 ENDIF 157 END DO 287 158 ENDIF 288 ! ! =============== 289 END DO ! End of slab 290 ! ! =============== 159 ! 160 ENDIF 291 161 292 IF( l_trdtra ) THEN ! Take the Asselin trend into account 293 ztrdt(:,:,:) = ztrdt(:,:,:) / ( 2.*rdt ) 294 ztrds(:,:,:) = ztrds(:,:,:) / ( 2.*rdt ) 162 #if defined key_agrif 163 ! Update tracer at AGRIF zoom boundaries 164 IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tra( kt ) ! children only 165 #endif 166 167 ! trends diagnostics : Asselin filter trend : (tb filtered - tb)/2dt 168 IF( l_trdtra ) THEN 169 DO jk = 1, jpkm1 170 zfact = 1.e0 / ( 2.*rdttra(jk) ) ! NB: euler case, (tb filtered - tb)=0 so 2dt always OK 171 ztrdt(:,:,jk) = ( tb(:,:,jk) - ztrdt(:,:,jk) ) * zfact 172 ztrds(:,:,jk) = ( sb(:,:,jk) - ztrds(:,:,jk) ) * zfact 173 END DO 295 174 CALL trd_mod( ztrdt, ztrds, jptra_trd_atf, 'TRA', kt ) 296 175 END IF 297 176 ! ! control print 298 177 IF(ln_ctl) CALL prt_ctl( tab3d_1=tn, clinfo1=' nxt - Tn: ', mask1=tmask, & 299 178 & tab3d_2=sn, clinfo2= ' Sn: ', mask2=tmask ) 300 #if defined key_agrif301 IF (.NOT.Agrif_Root()) CALL Agrif_Update_Tra( kt )302 #endif303 179 ! 304 180 END SUBROUTINE tra_nxt
Note: See TracChangeset
for help on using the changeset viewer.
