- Timestamp:
- 2020-05-12T14:51:37+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/src/OCE/TRA/tramle.F90
r12658 r12912 21 21 USE lbclnk ! lateral boundary condition / mpp link 22 22 23 ! where OSMOSIS_OBL is used with integrated FK 24 USE zdf_oce, ONLY : ln_zdfosm 25 USE zdfosm, ONLY : ln_osm_mle, hmle, dbdx_mle, dbdy_mle, mld_prof 26 23 27 IMPLICIT NONE 24 28 PRIVATE … … 56 60 CONTAINS 57 61 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) & 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 ! 100 ! 101 IF(ln_osm_mle.and.ln_zdfosm) THEN 102 ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 ) ! max level of the computation 103 ! 104 ! 105 SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts 106 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 113 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 120 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 127 END SELECT 128 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 140 ! 141 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 151 ENDIF 152 153 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 251 ! 252 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) & 231 291 & + 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 292 END DO 293 END DO 294 END DO 295 296 IF( cdtype == 'TRA') THEN !== outputs ==! 297 ! 298 IF (ln_osm_mle.and.ln_zdfosm) THEN 299 zLf_NH(:,:) = SQRT( rb_c * hmle(:,:) ) * r1_ft(:,:) ! Lf = N H / f 300 ELSE 301 zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:) ! Lf = N H / f 302 END IF 303 CALL iom_put( "Lf_NHpf" , zLf_NH ) ! Lf = N H / f 304 ! 305 ! divide by cross distance to give streamfunction with dimensions m^2/s 306 DO jk = 1, ikmax+1 307 zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:) 308 zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:) 309 END DO 310 CALL iom_put( "psiu_mle", zpsi_uw ) ! i-mle streamfunction 311 CALL iom_put( "psiv_mle", zpsi_vw ) ! j-mle streamfunction 312 ENDIF 313 ! 314 END SUBROUTINE tra_mle_trp 251 315 252 316
Note: See TracChangeset
for help on using the changeset viewer.