- Timestamp:
- 2020-05-14T21:46:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/TRA/tramle.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Property svn:externals
-
old new 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRA/tramle.F90
r12384 r12928 45 45 46 46 REAL(wp) :: r5_21 = 5.e0 / 21.e0 ! factor used in mle streamfunction computation 47 REAL(wp) :: rb_c ! ML buoyancy criteria = g rho_c /r au0 where rho_c is defined in zdfmld47 REAL(wp) :: rb_c ! ML buoyancy criteria = g rho_c /rho0 where rho_c is defined in zdfmld 48 48 REAL(wp) :: rc_f ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_mle=1 case 49 49 … … 52 52 53 53 !! * Substitutions 54 # include " vectopt_loop_substitute.h90"54 # include "do_loop_substitute.h90" 55 55 !!---------------------------------------------------------------------- 56 56 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 60 60 CONTAINS 61 61 62 SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype ) 63 !!---------------------------------------------------------------------- 64 !! *** ROUTINE tra_mle_trp *** 65 !! 66 !! ** Purpose : Add to the transport the Mixed Layer Eddy induced transport 67 !! 68 !! ** Method : The 3 components of the Mixed Layer Eddy (MLE) induced 69 !! transport are computed as follows : 70 !! zu_mle = dk[ zpsi_uw ] 71 !! zv_mle = dk[ zpsi_vw ] 72 !! zw_mle = - di[ zpsi_uw ] - dj[ zpsi_vw ] 73 !! where zpsi is the MLE streamfunction at uw and vw points (see the doc) 74 !! and added to the input velocity : 75 !! p.n = p.n + z._mle 76 !! 77 !! ** Action : - (pun,pvn,pwn) increased by the mle transport 78 !! CAUTION, the transport is not updated at the last line/raw 79 !! this may be a problem for some advection schemes 80 !! 81 !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 82 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 83 !!---------------------------------------------------------------------- 84 INTEGER , INTENT(in ) :: kt ! ocean time-step index 85 INTEGER , INTENT(in ) :: kit000 ! first time step index 86 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 87 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu ! in : 3 ocean transport components 88 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pv ! out: same 3 transport components 89 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pw ! increased by the MLE induced transport 90 ! 91 INTEGER :: ji, jj, jk ! dummy loop indices 92 INTEGER :: ii, ij, ik, ikmax ! local integers 93 REAL(wp) :: zcuw, zmuw, zc ! local scalar 94 REAL(wp) :: zcvw, zmvw ! - - 95 INTEGER , DIMENSION(jpi,jpj) :: inml_mle 96 REAL(wp), DIMENSION(jpi,jpj) :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH 97 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw 98 !!---------------------------------------------------------------------- 99 ! 62 SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype, Kmm ) 63 !!---------------------------------------------------------------------- 64 !! *** ROUTINE tra_mle_trp *** 65 !! 66 !! ** Purpose : Add to the transport the Mixed Layer Eddy induced transport 67 !! 68 !! ** Method : The 3 components of the Mixed Layer Eddy (MLE) induced 69 !! transport are computed as follows : 70 !! zu_mle = dk[ zpsi_uw ] 71 !! zv_mle = dk[ zpsi_vw ] 72 !! zw_mle = - di[ zpsi_uw ] - dj[ zpsi_vw ] 73 !! where zpsi is the MLE streamfunction at uw and vw points (see the doc) 74 !! and added to the input velocity : 75 !! p.n = p.n + z._mle 76 !! 77 !! ** Action : - (pu,pv,pw) increased by the mle transport 78 !! CAUTION, the transport is not updated at the last line/raw 79 !! this may be a problem for some advection schemes 80 !! 81 !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 82 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 83 !!---------------------------------------------------------------------- 84 INTEGER , INTENT(in ) :: kt ! ocean time-step index 85 INTEGER , INTENT(in ) :: kit000 ! first time step index 86 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 87 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 88 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu ! in : 3 ocean transport components 89 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pv ! out: same 3 transport components 90 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pw ! increased by the MLE induced transport 91 ! 92 INTEGER :: ji, jj, jk ! dummy loop indices 93 INTEGER :: ii, ij, ik, ikmax ! local integers 94 REAL(wp) :: zcuw, zmuw, zc ! local scalar 95 REAL(wp) :: zcvw, zmvw ! - - 96 INTEGER , DIMENSION(jpi,jpj) :: inml_mle 97 REAL(wp), DIMENSION(jpi,jpj) :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH 98 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw 99 !!---------------------------------------------------------------------- 100 ! 100 101 ! 101 102 IF(ln_osm_mle.and.ln_zdfosm) THEN … … 105 106 SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts 106 107 CASE ( 0 ) != min of the 2 neighbour MLDs 107 DO jj = 1, jpjm1 108 DO ji = 1, fs_jpim1 ! vector opt. 109 zhu(ji,jj) = MIN( hmle(ji+1,jj), hmle(ji,jj) ) 110 zhv(ji,jj) = MIN( hmle(ji,jj+1), hmle(ji,jj) ) 111 END DO 112 END DO 108 DO_2D_10_10 109 zhu(ji,jj) = MIN( hmle(ji+1,jj), hmle(ji,jj) ) 110 zhv(ji,jj) = MIN( hmle(ji,jj+1), hmle(ji,jj) ) 111 END_2D 113 112 CASE ( 1 ) != average of the 2 neighbour MLDs 114 DO jj = 1, jpjm1 115 DO ji = 1, fs_jpim1 ! vector opt. 116 zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 117 zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) 118 END DO 119 END DO 113 DO_2D_10_10 114 zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 115 zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) 116 END_2D 120 117 CASE ( 2 ) != max of the 2 neighbour MLDs 121 DO jj = 1, jpjm1 122 DO ji = 1, fs_jpim1 ! vector opt. 123 zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 124 zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) 125 END DO 126 END DO 118 DO_2D_10_10 119 zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 120 zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) 121 END_2D 127 122 END SELECT 128 123 IF( nn_mle == 0 ) THEN ! Fox-Kemper et al. 2010 formulation 129 DO jj = 1, jpjm1 130 DO ji = 1, fs_jpim1 ! vector opt. 131 zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj) * e2u(ji,jj) & 132 & * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) & 133 & / ( MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) ) ) 134 ! 135 zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj) * e1v(ji,jj) & 136 & * dbdy_mle(ji,jj) * MIN( 111.e3_wp , e2v(ji,jj) ) & 137 & / ( MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) ) ) 138 END DO 139 END DO 124 DO_2D_10_10 125 zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj) * e2u(ji,jj) & 126 & * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) & 127 & / ( MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) ) ) 128 ! 129 zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj) * e1v(ji,jj) & 130 & * dbdy_mle(ji,jj) * MIN( 111.e3_wp , e2v(ji,jj) ) & 131 & / ( MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) ) ) 132 END_2D 140 133 ! 141 134 ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 142 DO jj = 1, jpjm1 143 DO ji = 1, fs_jpim1 ! vector opt. 144 zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2u(ji,jj) & 145 & * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) 146 ! 147 zpsim_v(ji,jj) = rc_f * zhv(ji,jj) * zhv(ji,jj) * e1v(ji,jj) & 148 & * dbdy_mle(ji,jj) * MIN( 111.e3_wp , e2v(ji,jj) ) 149 END DO 150 END DO 135 DO_2D_10_10 136 zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2u(ji,jj) & 137 & * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) 138 ! 139 zpsim_v(ji,jj) = rc_f * zhv(ji,jj) * zhv(ji,jj) * e1v(ji,jj) & 140 & * dbdy_mle(ji,jj) * MIN( 111.e3_wp , e2v(ji,jj) ) 141 END_2D 151 142 ENDIF 152 153 143 ELSE !do not use osn_mle 154 ! !== MLD used for MLE ==! 155 ! ! compute from the 10m density to deal with the diurnal cycle 156 inml_mle(:,:) = mbkt(:,:) + 1 ! init. to number of ocean w-level (T-level + 1) 157 IF ( nla10 > 0 ) THEN ! avoid case where first level is thicker than 10m 158 DO jk = jpkm1, nlb10, -1 ! from the bottom to nlb10 (10m) 159 DO jj = 1, jpj 160 DO ji = 1, jpi ! index of the w-level at the ML based 161 IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle ) inml_mle(ji,jj) = jk ! Mixed layer 162 END DO 163 END DO 164 END DO 165 ENDIF 166 ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 ) ! max level of the computation 167 168 ! 169 ! 170 zmld(:,:) = 0._wp !== Horizontal shape of the MLE ==! 171 zbm (:,:) = 0._wp 172 zn2 (:,:) = 0._wp 173 DO jk = 1, ikmax ! MLD and mean buoyancy and N2 over the mixed layer 174 DO jj = 1, jpj 175 DO ji = 1, jpi 176 zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 177 zmld(ji,jj) = zmld(ji,jj) + zc 178 zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0 179 zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 180 END DO 181 END DO 182 END DO 183 184 SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts 185 CASE ( 0 ) != min of the 2 neighbour MLDs 186 DO jj = 1, jpjm1 187 DO ji = 1, fs_jpim1 ! vector opt. 188 zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) ) 189 zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) ) 190 END DO 191 END DO 192 CASE ( 1 ) != average of the 2 neighbour MLDs 193 DO jj = 1, jpjm1 194 DO ji = 1, fs_jpim1 ! vector opt. 195 zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp 196 zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp 197 END DO 198 END DO 199 CASE ( 2 ) != max of the 2 neighbour MLDs 200 DO jj = 1, jpjm1 201 DO ji = 1, fs_jpim1 ! vector opt. 202 zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) ) 203 zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) ) 204 END DO 205 END DO 206 END SELECT 207 ! ! convert density into buoyancy 208 zbm(:,:) = + grav * zbm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 209 ! 210 ! 211 ! !== Magnitude of the MLE stream function ==! 212 ! 213 ! di[bm] Ds 214 ! Psi = Ce H^2 ---------------- e2u mu(z) where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 ) 215 ! e1u Lf fu and the e2u for the "transport" 216 ! (not *e3u as divided by e3u at the end) 217 ! 218 IF( nn_mle == 0 ) THEN ! Fox-Kemper et al. 2010 formulation 219 DO jj = 1, jpjm1 220 DO ji = 1, fs_jpim1 ! vector opt. 221 zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 222 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) & 223 & / ( MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) ) ) 224 ! 225 zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) & 226 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) & 227 & / ( MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) ) ) 228 END DO 229 END DO 230 ! 231 ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 232 DO jj = 1, jpjm1 233 DO ji = 1, fs_jpim1 ! vector opt. 234 zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 235 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) 236 ! 237 zpsim_v(ji,jj) = rc_f * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) & 238 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) 239 END DO 240 END DO 241 ENDIF 242 ! 243 IF( nn_conv == 1 ) THEN ! No MLE in case of convection 244 DO jj = 1, jpjm1 245 DO ji = 1, fs_jpim1 ! vector opt. 246 IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp ) zpsim_u(ji,jj) = 0._wp 247 IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp ) zpsim_v(ji,jj) = 0._wp 248 END DO 249 END DO 250 ENDIF 144 ! !== MLD used for MLE ==! 145 ! ! compute from the 10m density to deal with the diurnal cycle 146 inml_mle(:,:) = mbkt(:,:) + 1 ! init. to number of ocean w-level (T-level + 1) 147 IF ( nla10 > 0 ) THEN ! avoid case where first level is thicker than 10m 148 DO_3DS_11_11( jpkm1, nlb10, -1 ) 149 IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle ) inml_mle(ji,jj) = jk ! Mixed layer 150 END_3D 151 ENDIF 152 ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 ) ! max level of the computation 153 154 ! 155 ! 156 zmld(:,:) = 0._wp !== Horizontal shape of the MLE ==! 157 zbm (:,:) = 0._wp 158 zn2 (:,:) = 0._wp 159 DO_3D_11_11( 1, ikmax ) 160 zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 161 zmld(ji,jj) = zmld(ji,jj) + zc 162 zbm (ji,jj) = zbm (ji,jj) + zc * (rho0 - rhop(ji,jj,jk) ) * r1_rho0 163 zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 164 END_3D 165 166 SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts 167 CASE ( 0 ) != min of the 2 neighbour MLDs 168 DO_2D_10_10 169 zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) ) 170 zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) ) 171 END_2D 172 CASE ( 1 ) != average of the 2 neighbour MLDs 173 DO_2D_10_10 174 zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp 175 zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp 176 END_2D 177 CASE ( 2 ) != max of the 2 neighbour MLDs 178 DO_2D_10_10 179 zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) ) 180 zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) ) 181 END_2D 182 END SELECT 183 ! ! convert density into buoyancy 184 zbm(:,:) = + grav * zbm(:,:) / MAX( e3t(:,:,1,Kmm), zmld(:,:) ) 185 ! 186 ! 187 ! !== Magnitude of the MLE stream function ==! 188 ! 189 ! di[bm] Ds 190 ! Psi = Ce H^2 ---------------- e2u mu(z) where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 ) 191 ! e1u Lf fu and the e2u for the "transport" 192 ! (not *e3u as divided by e3u at the end) 193 ! 194 IF( nn_mle == 0 ) THEN ! Fox-Kemper et al. 2010 formulation 195 DO_2D_10_10 196 zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 197 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) & 198 & / ( MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) ) ) 199 ! 200 zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) & 201 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) & 202 & / ( MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) ) ) 203 END_2D 204 ! 205 ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 206 DO_2D_10_10 207 zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 208 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) 209 ! 210 zpsim_v(ji,jj) = rc_f * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) & 211 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) 212 END_2D 213 ENDIF 214 ! 215 IF( nn_conv == 1 ) THEN ! No MLE in case of convection 216 DO_2D_10_10 217 IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp ) zpsim_u(ji,jj) = 0._wp 218 IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp ) zpsim_v(ji,jj) = 0._wp 219 END_2D 220 ENDIF 251 221 ! 252 222 ENDIF ! end of lm_osm_mle loop 253 ! !== structure function value at uw- and vw-points ==! 254 DO jj = 1, jpjm1 255 DO ji = 1, fs_jpim1 ! vector opt. 256 zhu(ji,jj) = 1._wp / MAX(zhu(ji,jj), rsmall) ! hu --> 1/hu 257 zhv(ji,jj) = 1._wp / MAX(zhv(ji,jj), rsmall) 258 END DO 259 END DO 260 ! 261 zpsi_uw(:,:,:) = 0._wp 262 zpsi_vw(:,:,:) = 0._wp 263 ! 264 DO jk = 2, ikmax ! start from 2 : surface value = 0 265 DO jj = 1, jpjm1 266 DO ji = 1, fs_jpim1 ! vector opt. 267 zcuw = 1._wp - ( gdepw_n(ji+1,jj,jk) + gdepw_n(ji,jj,jk) ) * zhu(ji,jj) 268 zcvw = 1._wp - ( gdepw_n(ji,jj+1,jk) + gdepw_n(ji,jj,jk) ) * zhv(ji,jj) 269 zcuw = zcuw * zcuw 270 zcvw = zcvw * zcvw 271 zmuw = MAX( 0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw ) ) 272 zmvw = MAX( 0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw ) ) 273 ! 274 zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk) 275 zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk) 276 END DO 277 END DO 278 END DO 279 ! 280 ! !== transport increased by the MLE induced transport ==! 281 DO jk = 1, ikmax 282 DO jj = 1, jpjm1 ! CAUTION pu,pv must be defined at row/column i=1 / j=1 283 DO ji = 1, fs_jpim1 ! vector opt. 284 pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 285 pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 286 END DO 287 END DO 288 DO jj = 2, jpjm1 289 DO ji = fs_2, fs_jpim1 ! vector opt. 290 pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk) & 291 & + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) ) 292 END DO 293 END DO 294 END DO 223 ! 224 ! !== structure function value at uw- and vw-points ==! 225 DO_2D_10_10 226 zhu(ji,jj) = 1._wp / MAX(zhu(ji,jj), rsmall) ! hu --> 1/hu 227 zhv(ji,jj) = 1._wp / MAX(zhv(ji,jj), rsmall) 228 END_2D 229 ! 230 zpsi_uw(:,:,:) = 0._wp 231 zpsi_vw(:,:,:) = 0._wp 232 ! 233 DO_3D_10_10( 2, ikmax ) 234 zcuw = 1._wp - ( gdepw(ji+1,jj,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhu(ji,jj) 235 zcvw = 1._wp - ( gdepw(ji,jj+1,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhv(ji,jj) 236 zcuw = zcuw * zcuw 237 zcvw = zcvw * zcvw 238 zmuw = MAX( 0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw ) ) 239 zmvw = MAX( 0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw ) ) 240 ! 241 zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk) 242 zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk) 243 END_3D 244 ! 245 ! !== transport increased by the MLE induced transport ==! 246 DO jk = 1, ikmax 247 DO_2D_10_10 248 pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 249 pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 250 END_2D 251 DO_2D_00_00 252 pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk) & 253 & + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) ) 254 END_2D 255 END DO 295 256 296 257 IF( cdtype == 'TRA') THEN !== outputs ==! … … 330 291 !!---------------------------------------------------------------------- 331 292 332 REWIND( numnam_ref ) ! Namelist namtra_mle in reference namelist : Tracer advection scheme333 293 READ ( numnam_ref, namtra_mle, IOSTAT = ios, ERR = 901) 334 294 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_mle in reference namelist' ) 335 295 336 REWIND( numnam_cfg ) ! Namelist namtra_mle in configuration namelist : Tracer advection scheme337 296 READ ( numnam_cfg, namtra_mle, IOSTAT = ios, ERR = 902 ) 338 297 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_mle in configuration namelist' ) … … 368 327 IF( ln_mle ) THEN ! MLE initialisation 369 328 ! 370 rb_c = grav * rn_rho_c_mle /r au0 ! Mixed Layer buoyancy criteria329 rb_c = grav * rn_rho_c_mle /rho0 ! Mixed Layer buoyancy criteria 371 330 IF(lwp) WRITE(numout,*) 372 331 IF(lwp) WRITE(numout,*) ' ML buoyancy criteria = ', rb_c, ' m/s2 ' … … 377 336 IF( ierr /= 0 ) CALL ctl_stop( 'tra_adv_mle_init: failed to allocate arrays' ) 378 337 z1_t2 = 1._wp / ( rn_time * rn_time ) 379 DO jj = 2, jpj ! "coriolis+ time^-1" at u- & v-points 380 DO ji = fs_2, jpi ! vector opt. 381 zfu = ( ff_f(ji,jj) + ff_f(ji,jj-1) ) * 0.5_wp 382 zfv = ( ff_f(ji,jj) + ff_f(ji-1,jj) ) * 0.5_wp 383 rfu(ji,jj) = SQRT( zfu * zfu + z1_t2 ) 384 rfv(ji,jj) = SQRT( zfv * zfv + z1_t2 ) 385 END DO 386 END DO 338 DO_2D_01_01 339 zfu = ( ff_f(ji,jj) + ff_f(ji,jj-1) ) * 0.5_wp 340 zfv = ( ff_f(ji,jj) + ff_f(ji-1,jj) ) * 0.5_wp 341 rfu(ji,jj) = SQRT( zfu * zfu + z1_t2 ) 342 rfv(ji,jj) = SQRT( zfv * zfv + z1_t2 ) 343 END_2D 387 344 CALL lbc_lnk_multi( 'tramle', rfu, 'U', 1. , rfv, 'V', 1. ) 388 345 !
Note: See TracChangeset
for help on using the changeset viewer.
