Changeset 2678 for branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO
- Timestamp:
- 2011-03-09T16:02:46+01:00 (13 years ago)
- Location:
- branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO
- Files:
-
- 3 edited
-
LIM_SRC_3/limthd_dif.F90 (modified) (2 diffs)
-
OPA_SRC/LDF/ldfslp.F90 (modified) (1 diff)
-
OPA_SRC/TRA/trazdf_imp.F90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r2612 r2678 79 79 !! profile of the ice/snow layers : z_i, z_s 80 80 !! total ice/snow thickness : ht_i_b, ht_s_b 81 !! 82 !! ** External : 83 !! 84 !! ** References : 85 !! 86 !! ** History : 87 !! (02-2003) Martin Vancoppenolle, Louvain-la-Neuve, Belgium 88 !! (06-2005) Martin Vancoppenolle, 3d version 89 !! (11-2006) Vectorized by Xavier Fettweis (UCL-ASTR) 90 !! (04-2007) Energy conservation tested by M. Vancoppenolle 81 91 !!------------------------------------------------------------------ 82 92 INTEGER , INTENT (in) :: & … … 333 343 END DO 334 344 END DO 345 ENDIF 346 347 IF ( thcon_i_swi .EQ. 1 ) THEN ! Pringle (0.011/2=0.0055) 348 DO layer = 1, nlay_i-1 349 DO ji = kideb , kiut 350 ztcond_i(ji,layer) = rcdic + 0.09*( s_i_b(ji,layer) & 351 + s_i_b(ji,layer+1) ) / MIN(-2.0*zeps, & 352 t_i_b(ji,layer)+t_i_b(ji,layer+1)-2.0*rtt) - & 353 0.0055* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt ) 354 ztcond_i(ji,layer) = MAX(ztcond_i(ji,layer),zkimin) 355 END DO 356 END DO 357 ENDIF 358 359 IF ( thcon_i_swi .EQ. 0 ) THEN ! Untersteiner 335 360 DO ji = kideb , kiut 336 361 ztcond_i(ji,nlay_i) = rcdic + zbeta*s_i_b(ji,nlay_i) / & -
branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r2633 r2678 815 815 WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1) 816 816 END SUBROUTINE ldf_slp 817 SUBROUTINE ldf_slp_grif( kt ) ! Dummy routine 818 INTEGER, INTENT(in) :: kt 819 WRITE(*,*) 'ldf_slp_grif: You should not have seen this print! error?', kt 820 END SUBROUTINE ldf_slp_grif 817 821 SUBROUTINE ldf_slp_init ! Dummy routine 818 822 END SUBROUTINE ldf_slp_init -
branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90
r2636 r2678 15 15 !! 3.2 ! 2009-03 (G. Madec) heat and salt content trends 16 16 !! 3.3 ! 2010-06 (C. Ethe, G. Madec) Merge TRA-TRC 17 !! - ! 2011-02 (A. Coward, C. Ethe, G. Madec) improvment of surface boundary condition 17 18 !!---------------------------------------------------------------------- 18 19 … … 40 41 PUBLIC tra_zdf_imp ! routine called by step.F90 41 42 43 REAL(wp) :: r_vvl ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise 44 42 45 !! * Substitutions 43 46 # include "domzgr_substitute.h90" … … 56 59 !! *** ROUTINE tra_zdf_imp *** 57 60 !! 58 !! ** Purpose : Compute the trend due to the vertical tracer diffusion 59 !! including the vertical component of lateral mixing (only for 2nd 60 !! order operator, for fourth order it is already computed and add 61 !! to the general trend in traldf.F) and add it to the general trend 62 !! of the tracer equations. 63 !! 64 !! ** Method : The vertical component of the lateral diffusive trends 65 !! is provided by a 2nd order operator rotated along neutral or geo- 66 !! potential surfaces to which an eddy induced advection can be 67 !! added. It is computed using before fields (forward in time) and 68 !! isopycnal or geopotential slopes computed in routine ldfslp. 69 !! 70 !! Second part: vertical trend associated with the vertical physics 71 !! =========== (including the vertical flux proportional to dk[t] 72 !! associated with the lateral mixing, through the 73 !! update of avt) 74 !! The vertical diffusion of the tracer t is given by: 75 !! difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) ) 61 !! ** Purpose : Compute the after tracer through a implicit computation 62 !! of the vertical tracer diffusion (including the vertical component 63 !! of lateral mixing (only for 2nd order operator, for fourth order 64 !! it is already computed and add to the general trend in traldf) 65 !! 66 !! ** Method : The vertical diffusion of the tracer t is given by: 67 !! difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) ) 76 68 !! It is computed using a backward time scheme (t=ta). 69 !! If lk_zdfddm=T, use avs for salinity or for passive tracers 77 70 !! Surface and bottom boundary conditions: no diffusive flux on 78 71 !! both tracers (bottom, applied through the masked field avt). 79 !! Add this trend to the general trend ta,sa : 80 !! ta = ta + dz( avt dz(t) ) 81 !! if lk_zdfddm=T, use avs for salinity or for passive tracers 82 !! (sa = sa + dz( avs dz(t) ) 83 !! 84 !! Third part: recover avt resulting from the vertical physics 85 !! ========== alone, for further diagnostics (for example to 86 !! compute the turbocline depth in zdfmxl.F90). 87 !! if lk_zdfddm=T, use avt = zavt 88 !! (avs = zavs if lk_zdfddm=T ) 89 !! 90 !! ** Action : - Update (ta) with before vertical diffusion trend 72 !! If iso-neutral mixing, add to avt the contribution due to lateral mixing. 73 !! 74 !! ** Action : - pta becomes the after tracer 91 75 !!--------------------------------------------------------------------- 92 76 USE oce , ONLY : zwd => ua ! ua used as workspace … … 103 87 !! 104 88 INTEGER :: ji, jj, jk, jn ! dummy loop indices 105 REAL(wp) :: z avi, zrhs, znvvl! local scalars89 REAL(wp) :: zrhs ! local scalars 106 90 REAL(wp) :: ze3tb, ze3tn, ze3ta ! variable vertical scale factors 107 91 !!--------------------------------------------------------------------- … … 116 100 IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype 117 101 IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ ' 118 zavi = 0._wp ! avoid warning at compilation phase when lk_ldfslp=F 102 ! 103 IF( lk_vvl ) THEN ; r_vvl = 1._wp ! Variable volume indicator 104 ELSE ; r_vvl = 0._wp 105 ENDIF 119 106 ENDIF 120 107 ! 121 ! I. Local initialization 122 ! ----------------------- 123 zwd(1,:, : ) = 0._wp ; zwd(jpi,:,:) = 0._wp 124 zws(1,:, : ) = 0._wp ; zws(jpi,:,:) = 0._wp 125 zwi(1,:, : ) = 0._wp ; zwi(jpi,:,:) = 0._wp 126 zwt(1,:, : ) = 0._wp ; zwt(jpi,:,:) = 0._wp 127 zwt(:,:,jpk) = 0._wp ; zwt( : ,:,1) = 0._wp 128 129 ! I.1 Variable volume : to take into account vertical variable vertical scale factors 130 ! ------------------- 131 IF( lk_vvl ) THEN ; znvvl = 1._wp 132 ELSE ; znvvl = 0._wp 133 ENDIF 134 135 ! II. Vertical trend associated with the vertical physics 136 ! ======================================================= 137 ! (including the vertical flux proportional to dk[t] associated 138 ! with the lateral mixing, through the avt update) 139 ! dk[ avt dk[ (t,s) ] ] diffusive trends 140 141 ! 142 ! II.0 Matrix construction 143 ! ------------------------ 144 DO jn = 1, kjpt 108 ! ! ============= ! 109 DO jn = 1, kjpt ! tracer loop ! 110 ! ! ============= ! 145 111 ! 146 112 ! Matrix construction 147 ! ------------------------ 148 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN 113 ! -------------------- 114 ! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer 115 ! 116 IF( ( cdtype == 'TRA' .AND. ( ( jn == jp_tem ) .OR. ( jn == jp_sal .AND. lk_zdfddm ) ) ) .OR. & 117 & ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN 118 ! 119 ! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers 120 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN ; zwt(:,:,2:jpk) = avt (:,:,2:jpk) 121 ELSE ; zwt(:,:,2:jpk) = fsavs(:,:,2:jpk) 122 ENDIF 123 zwt(:,:,1) = 0._wp 124 ! 149 125 #if defined key_ldfslp 150 IF( ln_traldf_grif ) THEN 126 ! isoneutral diffusion: add the contribution 127 IF( ln_traldf_grif ) THEN ! Griffies isoneutral diff 151 128 DO jk = 2, jpkm1 152 129 DO jj = 2, jpjm1 153 130 DO ji = fs_2, fs_jpim1 ! vector opt. 154 ! zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk) ! vertical mixing coef. due to lateral mixing 155 zavi = ah_wslp2(ji,jj,jk) ! vertical mixing coef. due to lateral mixing 156 zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on temperature) 131 zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk) 157 132 END DO 158 133 END DO 159 134 END DO 160 ! update and save of avt (and avs if double diffusive mixing) 161 ELSE IF( l_traldf_rot ) THEN 135 ELSE IF( l_traldf_rot ) THEN ! standard isoneutral diff 162 136 DO jk = 2, jpkm1 163 137 DO jj = 2, jpjm1 164 138 DO ji = fs_2, fs_jpim1 ! vector opt. 165 zavi = fsahtw(ji,jj,jk) & ! vertical mixing coef. due to lateral mixing 166 & * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 167 & + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) 168 zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on temperature) 139 zwt(ji,jj,jk) = zwt(ji,jj,jk) + fsahtw(ji,jj,jk) & 140 & * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 141 & + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) 169 142 END DO 170 143 END DO 171 144 END DO 172 ELSE ! no rotation but key_ldfslp defined173 zwt(:,:,:) = avt(:,:,:)174 145 ENDIF 175 #else176 ! No isopycnal diffusion177 zwt(:,:,:) = avt(:,:,:)178 146 #endif 179 ! Diagonal, inferior, superior (including the bottom boundary condition via avtmasked)147 ! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked) 180 148 DO jk = 1, jpkm1 181 149 DO jj = 2, jpjm1 182 150 DO ji = fs_2, fs_jpim1 ! vector opt. 183 ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,jk) ! after scale factor at T-point184 ze3tn = znvvl + ( 1. - znvvl ) * fse3t_n(ji,jj,jk) ! now scale factor at T-point151 ze3ta = ( 1. - r_vvl ) + r_vvl * fse3t_a(ji,jj,jk) ! after scale factor at T-point 152 ze3tn = r_vvl + ( 1. - r_vvl ) * fse3t_n(ji,jj,jk) ! now scale factor at T-point 185 153 zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) ) 186 154 zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) ) … … 189 157 END DO 190 158 END DO 159 ! 160 !! Matrix inversion from the first level 161 !!---------------------------------------------------------------------- 162 ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk ) 163 ! 164 ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 ) 165 ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 ) 166 ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 ) 167 ! ( ... )( ... ) ( ... ) 168 ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) 169 ! 170 ! m is decomposed in the product of an upper and lower triangular matrix. 171 ! The 3 diagonal terms are in 3d arrays: zwd, zws, zwi. 172 ! Suffices i,s and d indicate "inferior" (below diagonal), diagonal 173 ! and "superior" (above diagonal) components of the tridiagonal system. 174 ! The solution will be in the 4d array pta. 175 ! The 3d array zwt is used as a work space array. 176 ! En route to the solution pta is used a to evaluate the rhs and then 177 ! used as a work space array: its value is modified. 178 ! 191 179 ! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) 180 ! done once for all passive tracers (so included in the IF instruction) 192 181 DO jj = 2, jpjm1 193 182 DO ji = fs_2, fs_jpim1 … … 203 192 END DO 204 193 ! 205 ELSE IF( ( cdtype == 'TRA' .AND. jn == jp_sal ) .OR. ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN206 #if defined key_ldfslp207 IF( ln_traldf_grif ) THEN208 DO jk = 2, jpkm1209 DO jj = 2, jpjm1210 DO ji = fs_2, fs_jpim1 ! vector opt.211 zavi = ah_wslp2(ji,jj,jk) ! vertical mixing coef. due to lateral mixing212 ! zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk) ! vertical mixing coef. due to lateral mixing213 zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on temperature)214 END DO215 END DO216 END DO217 ELSE IF( l_traldf_rot ) THEN218 DO jk = 2, jpkm1219 DO jj = 2, jpjm1220 DO ji = fs_2, fs_jpim1 ! vector opt.221 zavi = fsahtw(ji,jj,jk) & ! vertical mixing coef. due to lateral mixing222 & * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) &223 & + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) )224 zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on salinity)225 END DO226 END DO227 END DO228 ELSE ! no rotation but key_ldfslp defined229 zwt(:,:,:) = fsavs(:,:,:)230 ENDIF231 #else232 ! No isopycnal diffusion233 zwt(:,:,:) = fsavs(:,:,:)234 #endif235 ! Diagonal, inferior, superior (including the bottom boundary condition via avt masked)236 DO jk = 1, jpkm1237 DO jj = 2, jpjm1238 DO ji = fs_2, fs_jpim1 ! vector opt.239 ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,jk) ! after scale factor at T-point240 ze3tn = znvvl + ( 1. - znvvl ) * fse3t_n(ji,jj,jk) ! now scale factor at T-point241 zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) )242 zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )243 zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)244 END DO245 END DO246 END DO247 ! Surface boudary conditions248 DO jj = 2, jpjm1249 DO ji = fs_2, fs_jpim1 ! vector opt.250 ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,1) ! after scale factor at T-point251 zwi(ji,jj,1) = 0._wp252 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)253 END DO254 END DO255 !256 ! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k)257 DO jj = 2, jpjm1258 DO ji = fs_2, fs_jpim1259 zwt(ji,jj,1) = zwd(ji,jj,1)260 END DO261 END DO262 DO jk = 2, jpkm1263 DO jj = 2, jpjm1264 DO ji = fs_2, fs_jpim1265 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)266 END DO267 END DO268 END DO269 !270 194 END IF 271 272 ! II.1. Vertical diffusion on tracer 273 ! --------------------------- 274 195 ! 275 196 ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1 276 197 DO jj = 2, jpjm1 277 198 DO ji = fs_2, fs_jpim1 278 ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1)279 ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1)199 ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,1) 200 ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,1) 280 201 pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn) 281 202 END DO … … 284 205 DO jj = 2, jpjm1 285 206 DO ji = fs_2, fs_jpim1 286 ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk)287 ze3tn = ( 1. - znvvl ) + znvvl * fse3t (ji,jj,jk)207 ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,jk) 208 ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t (ji,jj,jk) 288 209 zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn) ! zrhs=right hand side 289 210 pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn) … … 292 213 END DO 293 214 294 ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 215 ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer) 295 216 DO jj = 2, jpjm1 296 217 DO ji = fs_2, fs_jpim1 … … 306 227 END DO 307 228 END DO 308 ! 309 END DO 229 ! ! ================= ! 230 END DO ! end tracer loop ! 231 ! ! ================= ! 310 232 ! 311 233 IF(wrk_not_released(3, 1,2))THEN
Note: See TracChangeset
for help on using the changeset viewer.
