Changeset 12217
- Timestamp:
- 2019-12-12T17:01:26+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRA/tramle.F90
r12178 r12217 21 21 USE lbclnk ! lateral boundary condition / mpp link 22 22 23 USE zdfosm, ONLY : ln_osm_mle, hmle, dbdx_mle, dbdy_mle, mld_prof 24 23 25 IMPLICIT NONE 24 26 PRIVATE … … 56 58 CONTAINS 57 59 58 SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype ) 59 !!---------------------------------------------------------------------- 60 !! *** ROUTINE tra_mle_trp *** 61 !! 62 !! ** Purpose : Add to the transport the Mixed Layer Eddy induced transport 63 !! 64 !! ** Method : The 3 components of the Mixed Layer Eddy (MLE) induced 65 !! transport are computed as follows : 66 !! zu_mle = dk[ zpsi_uw ] 67 !! zv_mle = dk[ zpsi_vw ] 68 !! zw_mle = - di[ zpsi_uw ] - dj[ zpsi_vw ] 69 !! where zpsi is the MLE streamfunction at uw and vw points (see the doc) 70 !! and added to the input velocity : 71 !! p.n = p.n + z._mle 72 !! 73 !! ** Action : - (pun,pvn,pwn) increased by the mle transport 74 !! CAUTION, the transport is not updated at the last line/raw 75 !! this may be a problem for some advection schemes 76 !! 77 !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 78 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 79 !!---------------------------------------------------------------------- 80 INTEGER , INTENT(in ) :: kt ! ocean time-step index 81 INTEGER , INTENT(in ) :: kit000 ! first time step index 82 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 83 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu ! in : 3 ocean transport components 84 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pv ! out: same 3 transport components 85 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pw ! increased by the MLE induced transport 86 ! 87 INTEGER :: ji, jj, jk ! dummy loop indices 88 INTEGER :: ii, ij, ik, ikmax ! local integers 89 REAL(wp) :: zcuw, zmuw, zc ! local scalar 90 REAL(wp) :: zcvw, zmvw ! - - 91 INTEGER , DIMENSION(jpi,jpj) :: inml_mle 92 REAL(wp), DIMENSION(jpi,jpj) :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH 93 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw 94 !!---------------------------------------------------------------------- 95 ! 96 ! !== MLD used for MLE ==! 97 ! ! compute from the 10m density to deal with the diurnal cycle 98 inml_mle(:,:) = mbkt(:,:) + 1 ! init. to number of ocean w-level (T-level + 1) 99 IF ( nla10 > 0 ) THEN ! avoid case where first level is thicker than 10m 100 DO jk = jpkm1, nlb10, -1 ! from the bottom to nlb10 (10m) 101 DO jj = 1, jpj 102 DO ji = 1, jpi ! index of the w-level at the ML based 103 IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle ) inml_mle(ji,jj) = jk ! Mixed layer 104 END DO 105 END DO 106 END DO 107 ENDIF 108 ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 ) ! max level of the computation 109 ! 110 ! 111 zmld(:,:) = 0._wp !== Horizontal shape of the MLE ==! 112 zbm (:,:) = 0._wp 113 zn2 (:,:) = 0._wp 114 DO jk = 1, ikmax ! MLD and mean buoyancy and N2 over the mixed layer 115 DO jj = 1, jpj 116 DO ji = 1, jpi 117 zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 118 zmld(ji,jj) = zmld(ji,jj) + zc 119 zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0 120 zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 121 END DO 122 END DO 123 END DO 124 125 SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts 126 CASE ( 0 ) != min of the 2 neighbour MLDs 127 DO jj = 1, jpjm1 128 DO ji = 1, fs_jpim1 ! vector opt. 129 zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) ) 130 zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) ) 131 END DO 132 END DO 133 CASE ( 1 ) != average of the 2 neighbour MLDs 134 DO jj = 1, jpjm1 135 DO ji = 1, fs_jpim1 ! vector opt. 136 zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp 137 zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp 138 END DO 139 END DO 140 CASE ( 2 ) != max of the 2 neighbour MLDs 141 DO jj = 1, jpjm1 142 DO ji = 1, fs_jpim1 ! vector opt. 143 zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) ) 144 zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) ) 145 END DO 146 END DO 147 END SELECT 148 ! ! convert density into buoyancy 149 zbm(:,:) = + grav * zbm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 150 ! 151 ! 152 ! !== Magnitude of the MLE stream function ==! 153 ! 154 ! di[bm] Ds 155 ! Psi = Ce H^2 ---------------- e2u mu(z) where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 ) 156 ! e1u Lf fu and the e2u for the "transport" 157 ! (not *e3u as divided by e3u at the end) 158 ! 159 IF( nn_mle == 0 ) THEN ! Fox-Kemper et al. 2010 formulation 160 DO jj = 1, jpjm1 161 DO ji = 1, fs_jpim1 ! vector opt. 162 zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 163 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) & 164 & / ( MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) ) ) 165 ! 166 zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) & 167 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) & 168 & / ( MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) ) ) 169 END DO 170 END DO 171 ! 172 ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 173 DO jj = 1, jpjm1 174 DO ji = 1, fs_jpim1 ! vector opt. 175 zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 176 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) 177 ! 178 zpsim_v(ji,jj) = rc_f * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) & 179 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) 180 END DO 181 END DO 182 ENDIF 183 ! 184 IF( nn_conv == 1 ) THEN ! No MLE in case of convection 185 DO jj = 1, jpjm1 186 DO ji = 1, fs_jpim1 ! vector opt. 187 IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp ) zpsim_u(ji,jj) = 0._wp 188 IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp ) zpsim_v(ji,jj) = 0._wp 189 END DO 190 END DO 191 ENDIF 192 ! 193 ! !== structure function value at uw- and vw-points ==! 194 DO jj = 1, jpjm1 195 DO ji = 1, fs_jpim1 ! vector opt. 196 zhu(ji,jj) = 1._wp / zhu(ji,jj) ! hu --> 1/hu 197 zhv(ji,jj) = 1._wp / zhv(ji,jj) 198 END DO 199 END DO 200 ! 201 zpsi_uw(:,:,:) = 0._wp 202 zpsi_vw(:,:,:) = 0._wp 203 ! 204 DO jk = 2, ikmax ! start from 2 : surface value = 0 205 DO jj = 1, jpjm1 206 DO ji = 1, fs_jpim1 ! vector opt. 207 zcuw = 1._wp - ( gdepw_n(ji+1,jj,jk) + gdepw_n(ji,jj,jk) ) * zhu(ji,jj) 208 zcvw = 1._wp - ( gdepw_n(ji,jj+1,jk) + gdepw_n(ji,jj,jk) ) * zhv(ji,jj) 209 zcuw = zcuw * zcuw 210 zcvw = zcvw * zcvw 211 zmuw = MAX( 0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw ) ) 212 zmvw = MAX( 0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw ) ) 213 ! 214 zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk) 215 zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk) 216 END DO 217 END DO 218 END DO 219 ! 220 ! !== transport increased by the MLE induced transport ==! 221 DO jk = 1, ikmax 222 DO jj = 1, jpjm1 ! CAUTION pu,pv must be defined at row/column i=1 / j=1 223 DO ji = 1, fs_jpim1 ! vector opt. 224 pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 225 pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 226 END DO 227 END DO 228 DO jj = 2, jpjm1 229 DO ji = fs_2, fs_jpim1 ! vector opt. 230 pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk) & 60 SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype ) 61 !!---------------------------------------------------------------------- 62 !! *** ROUTINE tra_mle_trp *** 63 !! 64 !! ** Purpose : Add to the transport the Mixed Layer Eddy induced transport 65 !! 66 !! ** Method : The 3 components of the Mixed Layer Eddy (MLE) induced 67 !! transport are computed as follows : 68 !! zu_mle = dk[ zpsi_uw ] 69 !! zv_mle = dk[ zpsi_vw ] 70 !! zw_mle = - di[ zpsi_uw ] - dj[ zpsi_vw ] 71 !! where zpsi is the MLE streamfunction at uw and vw points (see the doc) 72 !! and added to the input velocity : 73 !! p.n = p.n + z._mle 74 !! 75 !! ** Action : - (pun,pvn,pwn) increased by the mle transport 76 !! CAUTION, the transport is not updated at the last line/raw 77 !! this may be a problem for some advection schemes 78 !! 79 !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 80 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 81 !!---------------------------------------------------------------------- 82 INTEGER , INTENT(in ) :: kt ! ocean time-step index 83 INTEGER , INTENT(in ) :: kit000 ! first time step index 84 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 85 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu ! in : 3 ocean transport components 86 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pv ! out: same 3 transport components 87 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pw ! increased by the MLE induced transport 88 ! 89 INTEGER :: ji, jj, jk ! dummy loop indices 90 INTEGER :: ii, ij, ik, ikmax ! local integers 91 REAL(wp) :: zcuw, zmuw, zc ! local scalar 92 REAL(wp) :: zcvw, zmvw ! - - 93 INTEGER , DIMENSION(jpi,jpj) :: inml_mle 94 REAL(wp), DIMENSION(jpi,jpj) :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH 95 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw 96 !!---------------------------------------------------------------------- 97 ! 98 ! 99 IF(ln_osm_mle) THEN 100 ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 ) ! max level of the computation 101 ! 102 ! 103 SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts 104 CASE ( 0 ) != min of the 2 neighbour MLDs 105 DO jj = 1, jpjm1 106 DO ji = 1, fs_jpim1 ! vector opt. 107 zhu(ji,jj) = MIN( hmle(ji+1,jj), hmle(ji,jj) ) 108 zhv(ji,jj) = MIN( hmle(ji,jj+1), hmle(ji,jj) ) 109 END DO 110 END DO 111 CASE ( 1 ) != average of the 2 neighbour MLDs 112 DO jj = 1, jpjm1 113 DO ji = 1, fs_jpim1 ! vector opt. 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 DO 117 END DO 118 CASE ( 2 ) != max of the 2 neighbour MLDs 119 DO jj = 1, jpjm1 120 DO ji = 1, fs_jpim1 ! vector opt. 121 zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 122 zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) 123 END DO 124 END DO 125 END SELECT 126 IF( nn_mle == 0 ) THEN ! Fox-Kemper et al. 2010 formulation 127 DO jj = 1, jpjm1 128 DO ji = 1, fs_jpim1 ! vector opt. 129 zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj) * e2u(ji,jj) & 130 & * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) & 131 & / ( MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) ) ) 132 ! 133 zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj) * e1v(ji,jj) & 134 & * dbdy_mle(ji,jj) * MIN( 111.e3_wp , e2v(ji,jj) ) & 135 & / ( MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) ) ) 136 END DO 137 END DO 138 ! 139 ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 140 DO jj = 1, jpjm1 141 DO ji = 1, fs_jpim1 ! vector opt. 142 zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2u(ji,jj) & 143 & * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) 144 ! 145 zpsim_v(ji,jj) = rc_f * zhv(ji,jj) * zhv(ji,jj) * e1v(ji,jj) & 146 & * dbdy_mle(ji,jj) * MIN( 111.e3_wp , e2v(ji,jj) ) 147 END DO 148 END DO 149 ENDIF 150 151 ELSE !do not use osn_mle 152 ! !== MLD used for MLE ==! 153 ! ! convert density into buoyancy 154 inml_mle(:,:) = mbkt(:,:) + 1 ! init. to number of ocean w-level (T-level + 1) 155 IF ( nla10 > 0 ) THEN ! avoid case where first level is thicker than 10m 156 DO jk = jpkm1, nlb10, -1 ! from the bottom to nlb10 (10m) 157 DO jj = 1, jpj 158 DO ji = 1, jpi ! index of the w-level at the ML based 159 IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle ) inml_mle(ji,jj) = jk ! Mixed layer 160 END DO 161 END DO 162 END DO 163 ENDIF 164 165 ! 166 ! 167 zmld(:,:) = 0._wp !== Horizontal shape of the MLE ==! 168 zbm (:,:) = 0._wp 169 zn2 (:,:) = 0._wp 170 DO jk = 1, ikmax ! MLD and mean buoyancy and N2 over the mixed layer 171 DO jj = 1, jpj 172 DO ji = 1, jpi 173 zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 174 zmld(ji,jj) = zmld(ji,jj) + zc 175 zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0 176 zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 177 END DO 178 END DO 179 END DO 180 181 SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts 182 CASE ( 0 ) != min of the 2 neighbour MLDs 183 DO jj = 1, jpjm1 184 DO ji = 1, fs_jpim1 ! vector opt. 185 zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) ) 186 zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) ) 187 END DO 188 END DO 189 CASE ( 1 ) != average of the 2 neighbour MLDs 190 DO jj = 1, jpjm1 191 DO ji = 1, fs_jpim1 ! vector opt. 192 zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp 193 zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp 194 END DO 195 END DO 196 CASE ( 2 ) != max of the 2 neighbour MLDs 197 DO jj = 1, jpjm1 198 DO ji = 1, fs_jpim1 ! vector opt. 199 zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) ) 200 zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) ) 201 END DO 202 END DO 203 END SELECT 204 ! ! convert density into buoyancy 205 ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 ) ! max level of the computation 206 zbm(:,:) = + grav * zbm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 207 ! 208 ! 209 ! !== Magnitude of the MLE stream function ==! 210 ! 211 ! di[bm] Ds 212 ! Psi = Ce H^2 ---------------- e2u mu(z) where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 ) 213 ! e1u Lf fu and the e2u for the "transport" 214 ! (not *e3u as divided by e3u at the end) 215 ! 216 IF( nn_mle == 0 ) THEN ! Fox-Kemper et al. 2010 formulation 217 DO jj = 1, jpjm1 218 DO ji = 1, fs_jpim1 ! vector opt. 219 zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 220 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) & 221 & / ( MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) ) ) 222 ! 223 zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) & 224 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) & 225 & / ( MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) ) ) 226 END DO 227 END DO 228 ! 229 ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 230 DO jj = 1, jpjm1 231 DO ji = 1, fs_jpim1 ! vector opt. 232 zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 233 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) 234 ! 235 zpsim_v(ji,jj) = rc_f * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) & 236 & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) 237 END DO 238 END DO 239 ENDIF 240 ! 241 IF( nn_conv == 1 ) THEN ! No MLE in case of convection 242 DO jj = 1, jpjm1 243 DO ji = 1, fs_jpim1 ! vector opt. 244 IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp ) zpsim_u(ji,jj) = 0._wp 245 IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp ) zpsim_v(ji,jj) = 0._wp 246 END DO 247 END DO 248 ENDIF 249 ! 250 ENDIF ! end of lm_osm_mle loop 251 ! !== structure function value at uw- and vw-points ==! 252 DO jj = 1, jpjm1 253 DO ji = 1, fs_jpim1 ! vector opt. 254 zhu(ji,jj) = 1._wp / zhu(ji,jj) ! hu --> 1/hu 255 zhv(ji,jj) = 1._wp / zhv(ji,jj) 256 END DO 257 END DO 258 ! 259 zpsi_uw(:,:,:) = 0._wp 260 zpsi_vw(:,:,:) = 0._wp 261 ! 262 DO jk = 2, ikmax ! start from 2 : surface value = 0 263 DO jj = 1, jpjm1 264 DO ji = 1, fs_jpim1 ! vector opt. 265 zcuw = 1._wp - ( gdepw_n(ji+1,jj,jk) + gdepw_n(ji,jj,jk) ) * zhu(ji,jj) 266 zcvw = 1._wp - ( gdepw_n(ji,jj+1,jk) + gdepw_n(ji,jj,jk) ) * zhv(ji,jj) 267 zcuw = zcuw * zcuw 268 zcvw = zcvw * zcvw 269 zmuw = MAX( 0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw ) ) 270 zmvw = MAX( 0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw ) ) 271 ! 272 zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk) 273 zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk) 274 END DO 275 END DO 276 END DO 277 ! 278 ! !== transport increased by the MLE induced transport ==! 279 DO jk = 1, ikmax 280 DO jj = 1, jpjm1 ! CAUTION pu,pv must be defined at row/column i=1 / j=1 281 DO ji = 1, fs_jpim1 ! vector opt. 282 pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 283 pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 284 END DO 285 END DO 286 DO jj = 2, jpjm1 287 DO ji = fs_2, fs_jpim1 ! vector opt. 288 pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk) & 231 289 & + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) ) 232 END DO 233 END DO 234 END DO 235 236 IF( cdtype == 'TRA') THEN !== outputs ==! 237 ! 238 zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:) ! Lf = N H / f 239 CALL iom_put( "Lf_NHpf" , zLf_NH ) ! Lf = N H / f 240 ! 241 ! divide by cross distance to give streamfunction with dimensions m^2/s 242 DO jk = 1, ikmax+1 243 zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:) 244 zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:) 245 END DO 246 CALL iom_put( "psiu_mle", zpsi_uw ) ! i-mle streamfunction 247 CALL iom_put( "psiv_mle", zpsi_vw ) ! j-mle streamfunction 248 ENDIF 249 ! 250 END SUBROUTINE tra_mle_trp 290 END DO 291 END DO 292 END DO 293 294 IF( cdtype == 'TRA') THEN !== outputs ==! 295 ! 296 IF (ln_osm_mle) THEN 297 zLf_NH(:,:) = SQRT( rb_c * hmle(:,:) ) * r1_ft(:,:) ! Lf = N H / f 298 ELSE 299 zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:) ! Lf = N H / f 300 END IF 301 CALL iom_put( "Lf_NHpf" , zLf_NH ) ! Lf = N H / f 302 ! 303 ! divide by cross distance to give streamfunction with dimensions m^2/s 304 DO jk = 1, ikmax+1 305 zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:) 306 zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:) 307 END DO 308 CALL iom_put( "psiu_mle", zpsi_uw ) ! i-mle streamfunction 309 CALL iom_put( "psiv_mle", zpsi_vw ) ! j-mle streamfunction 310 ENDIF 311 ! 312 END SUBROUTINE tra_mle_trp 251 313 252 314 … … 290 352 WRITE(numout,*) ' Density difference used to define ML for FK rn_rho_c_mle = ', rn_rho_c_mle 291 353 ENDIF 292 ! 293 IF(lwp) THEN 354 355 356 IF( ln_osm_mle .AND. ln_mle ) THEN 357 WRITE(numout,*) 'WARNING: You are running with both OSM-FK and default FK' 358 ! CALL ctl_stop('STOP in zdf_osm_init: Cannot run with both OSM-FK and default FK') 359 END IF 360 361 IF(lwp) THEN 294 362 WRITE(numout,*) 295 363 IF( ln_mle ) THEN
Note: See TracChangeset
for help on using the changeset viewer.