Changeset 1110
- Timestamp:
- 2008-06-13T16:10:35+02:00 (16 years ago)
- Location:
- trunk/NEMO/OPA_SRC/TRA
- Files:
-
- 3 edited
-
tranxt.F90 (modified) (6 diffs)
-
trazdf.F90 (modified) (1 diff)
-
trazdf_exp.F90 (modified) (4 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 -
trunk/NEMO/OPA_SRC/TRA/trazdf.F90
r888 r1110 92 92 CASE ( 0 ) ! explicit scheme 93 93 CALL tra_zdf_exp ( kt, r2dt ) 94 IF( l_trdtra ) THEN ! save the vertical diffusive trends for further diagnostics95 ztrdt(:,:,:) = ta(:,:,:) - ztrdt(:,:,:)96 ztrds(:,:,:) = sa(:,:,:) - ztrds(:,:,:)97 CALL trd_mod( ztrdt, ztrds, jptra_trd_zdf, 'TRA', kt )98 ENDIF99 94 100 95 CASE ( 1 ) ! implicit scheme 101 96 CALL tra_zdf_imp ( kt, r2dt ) 102 IF( l_trdtra ) THEN ! save the vertical diffusive trends for further diagnostics103 DO jk = 1, jpkm1104 ztrdt(:,:,jk) = ( ( ta(:,:,jk) - tb(:,:,jk) ) / r2dt(jk) ) - ztrdt(:,:,jk)105 ztrds(:,:,jk) = ( ( sa(:,:,jk) - sb(:,:,jk) ) / r2dt(jk) ) - ztrds(:,:,jk)106 END DO107 CALL trd_mod( ztrdt, ztrds, jptra_trd_zdf, 'TRA', kt )108 ENDIF109 97 110 98 END SELECT 111 99 112 ! ! print mean trends (used for debugging) 100 IF( l_trdtra ) THEN ! save the vertical diffusive trends for further diagnostics 101 DO jk = 1, jpkm1 102 ztrdt(:,:,jk) = ( ( ta(:,:,jk) - tb(:,:,jk) ) / r2dt(jk) ) - ztrdt(:,:,jk) 103 ztrds(:,:,jk) = ( ( sa(:,:,jk) - sb(:,:,jk) ) / r2dt(jk) ) - ztrds(:,:,jk) 104 END DO 105 CALL trd_mod( ztrdt, ztrds, jptra_trd_zdf, 'TRA', kt ) 106 ENDIF 107 108 ! ! print mean trends (used for debugging) 113 109 IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' zdf - Ta: ', mask1=tmask, & 114 110 & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) -
trunk/NEMO/OPA_SRC/TRA/trazdf_exp.F90
r789 r1110 3 3 !! *** MODULE trazdf_exp *** 4 4 !! Ocean active tracers: vertical component of the tracer mixing trend using 5 !! a n explicit time-stepping (time spllitting scheme)5 !! a split-explicit time-stepping 6 6 !!============================================================================== 7 !! History : 8 !! 6.0 ! 90-10 (B. Blanke) Original code9 !! 7.0 ! 91-11 (G. Madec)10 !! ! 92-06 (M. Imbard) correction on tracer trend loops11 !! ! 96-01 (G. Madec) statement function for e312 !! ! 97-05 (G. Madec) vertical component of isopycnal13 !! ! 97-07 (G. Madec) geopotential diffusion in s-coord14 !! ! 00-08 (G. Madec) double diffusive mixing15 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module16 !! 9.0 ! 04-08 (C. Talandier) New trendsorganisation17 !! ! 05-11 (G. Madec) New organisation7 !! History : OPA ! 1990-10 (B. Blanke) Original code 8 !! 7.0 ! 1991-11 (G. Madec) 9 !! ! 1992-06 (M. Imbard) correction on tracer trend loops 10 !! ! 1996-01 (G. Madec) statement function for e3 11 !! ! 1997-05 (G. Madec) vertical component of isopycnal 12 !! ! 1997-07 (G. Madec) geopotential diffusion in s-coord 13 !! ! 2000-08 (G. Madec) double diffusive mixing 14 !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module 15 !! - ! 2004-08 (C. Talandier) New trends organisation 16 !! - ! 2005-11 (G. Madec) New organisation 17 !! 3.0 ! 2008-04 (G. Madec) leap-frog time stepping done in trazdf 18 18 !!---------------------------------------------------------------------- 19 !! tra_zdf_exp : update the tracer trend with the vertical diffusion 20 !! using an explicit time stepping 19 21 20 !!---------------------------------------------------------------------- 22 !! * Modules used 21 !! tra_zdf_exp : compute the tracer the vertical diffusion trend using a 22 !! split-explicit time stepping and provide the after tracer 23 !!---------------------------------------------------------------------- 23 24 USE oce ! ocean dynamics and active tracers 24 25 USE dom_oce ! ocean space and time domain 26 USE domvvl ! variablevolume levels 25 27 USE trdmod ! ocean active tracers trends 26 28 USE trdmod_oce ! ocean variables trends … … 33 35 PRIVATE 34 36 35 !! * Routine accessibility 36 PUBLIC tra_zdf_exp ! routine called by step.F90 37 PUBLIC tra_zdf_exp ! routine called by step.F90 37 38 38 39 !! * Substitutions 39 40 # include "domzgr_substitute.h90" 40 41 # include "zdfddm_substitute.h90" 42 # include "vectopt_loop_substitute.h90" 41 43 !!---------------------------------------------------------------------- 42 !! OPA 9.0 , LOCEAN-IPSL (2005)43 !! $ Header$44 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt44 !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 45 !! $Id:$ 46 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 45 47 !!---------------------------------------------------------------------- 46 48 … … 51 53 !! *** ROUTINE tra_zdf_exp *** 52 54 !! 53 !! ** Purpose : Compute the trend due to the vertical tracer mixing 54 !! using an explicit time stepping and add it to the general trend 55 !! of the tracer equations. 55 !! ** Purpose : Compute the after tracer fields due to the vertical 56 !! tracer mixing alone, and then due to the whole tracer trend. 56 57 !! 57 !! ** Method : The vertical diffusion of tracers (t & s) is given by: 58 !! difft = dz( avt dz(tb) ) = 1/e3t dk+1( avt/e3w dk(tb) ) 59 !! It is evaluated with an Euler scheme, using a time splitting 60 !! technique. 61 !! Surface and bottom boundary conditions: no diffusive flux on 62 !! both tracers (bottom, applied through the masked field avt). 63 !! Add this trend to the general trend ta,sa : 64 !! ta = ta + dz( avt dz(t) ) 65 !! (sa = sa + dz( avs dz(t) ) if lk_zdfddm= T) 58 !! ** Method : - The after tracer fields due to the vertical diffusion 59 !! of tracers alone is given by: 60 !! zwx = tb + p2dt difft 61 !! where difft = dz( avt dz(tb) ) = 1/e3t dk+1( avt/e3w dk(tb) ) 62 !! (if lk_zdfddm=T use avs on salinity instead of avt) 63 !! difft is evaluated with an Euler split-explit scheme using a 64 !! no flux boundary condition at both surface and bottomi boundaries. 65 !! (N.B. bottom condition is applied through the masked field avt). 66 !! - the after tracer fields due to the whole trend is 67 !! obtained in leap-frog environment by : 68 !! ta = zwx + p2dt ta 69 !! - in case of variable level thickness (lk_vvl=T) the 70 !! the leap-frog is applied on thickness weighted tracer. That is: 71 !! ta = [ tb*e3tb + e3tn*( zwx - tb + p2dt ta ) ] / e3tn 66 72 !! 67 !! ** Action : - Update (ta,sa) with the before vertical diffusion trend 73 !! ** Action : - after tracer fields (ta,sa) 74 !!--------------------------------------------------------------------- 75 INTEGER , INTENT(in) :: kt ! ocean time-step index 76 REAL(wp), INTENT(in), DIMENSION(jpk) :: p2dt ! vertical profile of tracer time-step 68 77 !! 69 !!--------------------------------------------------------------------- 70 !! * Arguments 71 INTEGER, INTENT( in ) :: kt ! ocean time-step index 72 REAL(wp), DIMENSION(jpk), INTENT( in ) :: & 73 p2dt ! vertical profile of tracer time-step 74 75 !! * Local declarations 76 INTEGER :: ji, jj, jk, jl ! dummy loop indices 77 REAL(wp) :: & 78 zlavmr, & ! temporary scalars 79 zave3r, ze3tr, & ! " " 80 zta, zsa ! " " 81 REAL(wp), DIMENSION(jpi,jpk) :: & 82 zwx, zwy, zwz, zww 78 INTEGER :: ji, jj, jk, jl ! dummy loop indices 79 REAL(wp) :: zlavmr, zave3r, ze3tr ! temporary scalars 80 REAL(wp) :: zta, zsa, ze3tb, zcoef ! temporary scalars 81 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwx, zwy, zwz, zww ! 3D workspace 83 82 !!--------------------------------------------------------------------- 84 83 … … 89 88 ENDIF 90 89 90 ! Initializations 91 ! --------------- 92 zlavmr = 1. / float( n_zdfexp ) ! Local constant 93 ! 94 zwy(:,:, 1 ) = 0.e0 ; zww(:,:, 1 ) = 0.e0 ! surface boundary conditions: no flux 95 zwy(:,:,jpk) = 0.e0 ; zww(:,:,jpk) = 0.e0 ! bottom boundary conditions: no flux 96 ! 97 zwx(:,:,:) = tb(:,:,:) ; zwz(:,:,:) = sb(:,:,:) ! zwx and zwz arrays set to before tracer values 91 98 92 ! 0. Local constant initialization 93 ! -------------------------------- 94 zlavmr = 1. / float( n_zdfexp ) 95 96 ! ! =============== 97 DO jj = 2, jpjm1 ! Vertical slab 98 ! ! =============== 99 ! 1. Initializations 100 ! ------------------ 101 102 ! Surface & bottom boundary conditions: no flux 103 DO ji = 2, jpim1 104 zwy(ji, 1 ) = 0.e0 105 zwy(ji,jpk) = 0.e0 106 zww(ji, 1 ) = 0.e0 107 zww(ji,jpk) = 0.e0 108 END DO 109 110 ! zwx and zwz arrays set to before tracer values 111 DO jk = 1, jpk 112 DO ji = 2, jpim1 113 zwx(ji,jk) = tb(ji,jj,jk) 114 zwz(ji,jk) = sb(ji,jj,jk) 115 END DO 116 END DO 117 118 119 ! 2. Time splitting loop 120 ! ---------------------- 121 122 DO jl = 1, n_zdfexp 123 124 ! first vertical derivative 125 IF( lk_zdfddm ) THEN ! double diffusion: avs /= avt 126 DO jk = 2, jpk 127 DO ji = 2, jpim1 128 zave3r = 1.e0 / fse3w(ji,jj,jk) 129 zwy(ji,jk) = avt(ji,jj,jk) * ( zwx(ji,jk-1) - zwx(ji,jk) ) * zave3r 130 zww(ji,jk) = fsavs(ji,jj,jk) * ( zwz(ji,jk-1) - zwz(ji,jk) ) * zave3r 131 END DO 132 END DO 133 ELSE ! default : avs = avt 134 DO jk = 2, jpk 135 DO ji = 2, jpim1 136 zave3r = avt(ji,jj,jk) / fse3w(ji,jj,jk) 137 zwy(ji,jk) = zave3r *(zwx(ji,jk-1) - zwx(ji,jk) ) 138 zww(ji,jk) = zave3r *(zwz(ji,jk-1) - zwz(ji,jk) ) 139 END DO 140 END DO 141 ENDIF 142 143 ! trend estimation at kt+l*2*rdt/n_zdfexp 144 DO jk = 1, jpkm1 145 DO ji = 2, jpim1 146 ze3tr = zlavmr / fse3t(ji,jj,jk) 147 ! 2nd vertical derivative 148 zta = ( zwy(ji,jk) - zwy(ji,jk+1) ) * ze3tr 149 zsa = ( zww(ji,jk) - zww(ji,jk+1) ) * ze3tr 150 ! update the tracer trends 151 ta(ji,jj,jk) = ta(ji,jj,jk) + zta 152 sa(ji,jj,jk) = sa(ji,jj,jk) + zsa 153 ! update tracer fields at kt+l*2*rdt/n_zdfexp 154 zwx(ji,jk) = zwx(ji,jk) + p2dt(jk) * zta * tmask(ji,jj,jk) 155 zwz(ji,jk) = zwz(ji,jk) + p2dt(jk) * zsa * tmask(ji,jj,jk) 99 ! Split-explicit loop (after tracer due to the vertical diffusion alone) 100 ! ------------------- 101 ! 102 DO jl = 1, n_zdfexp 103 ! ! first vertical derivative 104 DO jk = 2, jpk 105 DO jj = 2, jpjm1 106 DO ji = fs_2, fs_jpim1 ! vector opt. 107 zave3r = 1.e0 / fse3w(ji,jj,jk) 108 zwy(ji,jj,jk) = avt(ji,jj,jk) * ( zwx(ji,jj,jk-1) - zwx(ji,jj,jk) ) * zave3r 109 zww(ji,jj,jk) = fsavs(ji,jj,jk) * ( zwz(ji,jj,jk-1) - zwz(ji,jj,jk) ) * zave3r 156 110 END DO 157 111 END DO 158 112 END DO 159 ! ! =============== 160 END DO ! End of slab 161 ! ! =============== 113 ! 114 DO jk = 1, jpkm1 ! second vertical derivative ==> tracer at kt+l*2*rdt/n_zdfexp 115 DO jj = 2, jpjm1 116 DO ji = fs_2, fs_jpim1 ! vector opt. 117 ze3tr = zlavmr / fse3t(ji,jj,jk) 118 zwx(ji,jj,jk) = zwx(ji,jj,jk) + p2dt(jk) * ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) ) * ze3tr 119 zwz(ji,jj,jk) = zwz(ji,jj,jk) + p2dt(jk) * ( zww(ji,jj,jk) - zww(ji,jj,jk+1) ) * ze3tr 120 END DO 121 END DO 122 END DO 123 ! 124 END DO 162 125 126 ! After tracer due to all trends 127 ! ------------------------------ 128 IF( lk_vvl ) THEN ! variable level thickness : leap-frog on tracer*e3t 129 DO jk = 1, jpkm1 130 DO jj = 2, jpjm1 131 DO ji = fs_2, fs_jpim1 ! vector opt. 132 ze3tb = fsve3t(ji,jj,jk) * ( 1 + sshb(ji,jj) * mut(ji,jj,jk) ) ! before e3t 133 zta = zwx(ji,jj,jk) - tb(ji,jj,jk) + p2dt(jk) * ta(ji,jj,jk) ! total trends * 2*rdt 134 zsa = zwz(ji,jj,jk) - sb(ji,jj,jk) + p2dt(jk) * sa(ji,jj,jk) 135 zcoef = 1.e0 / fse3t(ji,jj,jk) * tmask(ji,jj,jk) 136 ta(ji,jj,jk) = ( ze3tb * tb(ji,jj,jk) + fse3t(ji,jj,jk) * zta ) * zcoef 137 sa(ji,jj,jk) = ( ze3tb * sb(ji,jj,jk) + fse3t(ji,jj,jk) * zsa ) * zcoef 138 END DO 139 END DO 140 END DO 141 ELSE ! fixed level thickness : leap-frog on tracers 142 DO jk = 1, jpkm1 143 DO jj = 2, jpjm1 144 DO ji = fs_2, fs_jpim1 ! vector opt. 145 ta(ji,jj,jk) = ( zwx(ji,jj,jk) + p2dt(jk) * ta(ji,jj,jk) ) *tmask(ji,jj,jk) 146 sa(ji,jj,jk) = ( zwz(ji,jj,jk) + p2dt(jk) * sa(ji,jj,jk) ) *tmask(ji,jj,jk) 147 END DO 148 END DO 149 END DO 150 ENDIF 151 ! 163 152 END SUBROUTINE tra_zdf_exp 164 153
Note: See TracChangeset
for help on using the changeset viewer.
