- Timestamp:
- 2020-09-14T17:40:34+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@13382 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/DYN/dynzdf.F90
r11281 r13463 37 37 38 38 !! * Substitutions 39 # include "vectopt_loop_substitute.h90" 39 # include "do_loop_substitute.h90" 40 # include "domzgr_substitute.h90" 40 41 !!---------------------------------------------------------------------- 41 42 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 45 46 CONTAINS 46 47 47 SUBROUTINE dyn_zdf( kt )48 SUBROUTINE dyn_zdf( kt, Kbb, Kmm, Krhs, puu, pvv, Kaa ) 48 49 !!---------------------------------------------------------------------- 49 50 !! *** ROUTINE dyn_zdf *** … … 54 55 !! 55 56 !! ** Method : - Leap-Frog time stepping on all trends but the vertical mixing 56 !! u a = ub + 2*dt * uavector form or linear free surf.57 !! u a = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_aotherwise57 !! u(after) = u(before) + 2*dt * u(rhs) vector form or linear free surf. 58 !! u(after) = ( e3u_b*u(before) + 2*dt * e3u_n*u(rhs) ) / e3u_after otherwise 58 59 !! - update the after velocity with the implicit vertical mixing. 59 60 !! This requires to solver the following system: 60 !! u a = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_adk[ua] ]61 !! u(after) = u(after) + 1/e3u_after dk+1[ mi(avm) / e3uw_after dk[ua] ] 61 62 !! with the following surface/top/bottom boundary condition: 62 63 !! surface: wind stress input (averaged over kt-1/2 & kt+1/2) 63 64 !! top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90) 64 65 !! 65 !! ** Action : ( ua,va) after velocity66 !! ** Action : (puu(:,:,:,Kaa),pvv(:,:,:,Kaa)) after velocity 66 67 !!--------------------------------------------------------------------- 67 INTEGER, INTENT(in) :: kt ! ocean time-step index 68 INTEGER , INTENT( in ) :: kt ! ocean time-step index 69 INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs, Kaa ! ocean time level indices 70 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 68 71 ! 69 72 INTEGER :: ji, jj, jk ! dummy loop indices … … 90 93 ENDIF 91 94 ENDIF 92 ! !* set time step93 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! = rdt (restart with Euler time stepping)94 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2. * rdt ! = 2 rdt (leapfrog)95 ENDIF96 !97 95 ! !* explicit top/bottom drag case 98 IF( .NOT.ln_drgimp ) CALL zdf_drg_exp( kt, ub, vb, ua, va ) ! add top/bottom friction trend to (ua,va)96 IF( .NOT.ln_drgimp ) CALL zdf_drg_exp( kt, Kmm, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add top/bottom friction trend to (puu(Kaa),pvv(Kaa)) 99 97 ! 100 98 ! 101 99 IF( l_trddyn ) THEN !* temporary save of ta and sa trends 102 100 ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 103 ztrdu(:,:,:) = ua(:,:,:)104 ztrdv(:,:,:) = va(:,:,:)105 ENDIF 106 ! 107 ! !== RHS: Leap-Frog time stepping on all trends but the vertical mixing ==! (put in ua,va)101 ztrdu(:,:,:) = puu(:,:,:,Krhs) 102 ztrdv(:,:,:) = pvv(:,:,:,Krhs) 103 ENDIF 104 ! 105 ! !== RHS: Leap-Frog time stepping on all trends but the vertical mixing ==! (put in puu(:,:,:,Kaa),pvv(:,:,:,Kaa)) 108 106 ! 109 107 ! ! time stepping except vertical diffusion 110 108 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 111 DO jk = 1, jpkm1112 ua(:,:,jk) = ( ub(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk)113 va(:,:,jk) = ( vb(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk)114 END DO109 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 110 puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kbb) + rDt * puu(ji,jj,jk,Krhs) ) * umask(ji,jj,jk) 111 pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kbb) + rDt * pvv(ji,jj,jk,Krhs) ) * vmask(ji,jj,jk) 112 END_3D 115 113 ELSE ! applied on thickness weighted velocity 116 DO jk = 1, jpkm1 117 ua(:,:,jk) = ( e3u_b(:,:,jk) * ub(:,:,jk) & 118 & + r2dt * e3u_n(:,:,jk) * ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk) 119 va(:,:,jk) = ( e3v_b(:,:,jk) * vb(:,:,jk) & 120 & + r2dt * e3v_n(:,:,jk) * va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk) 121 END DO 114 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 115 puu(ji,jj,jk,Kaa) = ( e3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb ) & 116 & + rDt * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Krhs) ) & 117 & / e3u(ji,jj,jk,Kaa) * umask(ji,jj,jk) 118 pvv(ji,jj,jk,Kaa) = ( e3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb ) & 119 & + rDt * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Krhs) ) & 120 & / e3v(ji,jj,jk,Kaa) * vmask(ji,jj,jk) 121 END_3D 122 122 ENDIF 123 123 ! ! add top/bottom friction … … 125 125 ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does 126 126 ! not lead to the effective stress seen over the whole barotropic loop. 127 ! G. Madec : in linear free surface, e3u _a = e3u_n = e3u_0, so systematic use of e3u_a127 ! G. Madec : in linear free surface, e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa) 128 128 IF( ln_drgimp .AND. ln_dynspg_ts ) THEN 129 DO jk = 1, jpkm1! remove barotropic velocities130 ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk)131 va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk)132 END DO133 DO jj = 2, jpjm1 ! Add bottom/top stress due to barotropic component only134 DO ji = fs_2, fs_jpim1 ! vector opt.135 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points136 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points)137 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)138 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)139 ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua140 va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va141 END DO142 END DO129 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! remove barotropic velocities 130 puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - uu_b(ji,jj,Kaa) ) * umask(ji,jj,jk) 131 pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - vv_b(ji,jj,Kaa) ) * vmask(ji,jj,jk) 132 END_3D 133 DO_2D( 0, 0, 0, 0 ) 134 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 135 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 136 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) & 137 & + r_vvl * e3u(ji,jj,iku,Kaa) 138 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) & 139 & + r_vvl * e3v(ji,jj,ikv,Kaa) 140 puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua 141 pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va 142 END_2D 143 143 IF( ln_isfcav ) THEN ! Ocean cavities (ISF) 144 DO jj = 2, jpjm1145 DO ji = fs_2, fs_jpim1 ! vector opt.146 iku = miku(ji,jj) ! top ocean level at u- and v-points147 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points)148 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)149 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)150 ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua151 va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va152 END DO153 END DO144 DO_2D( 0, 0, 0, 0 ) 145 iku = miku(ji,jj) ! top ocean level at u- and v-points 146 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 147 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) & 148 & + r_vvl * e3u(ji,jj,iku,Kaa) 149 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) & 150 & + r_vvl * e3v(ji,jj,ikv,Kaa) 151 puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua 152 pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va 153 END_2D 154 154 END IF 155 155 ENDIF … … 158 158 ! 159 159 ! !* Matrix construction 160 zdt = r 2dt * 0.5160 zdt = rDt * 0.5 161 161 IF( ln_zad_Aimp ) THEN !! 162 162 SELECT CASE( nldf_dyn ) 163 163 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 164 DO jk = 1, jpkm1 165 DO jj = 2, jpjm1 166 DO ji = fs_2, fs_jpim1 ! vector opt. 167 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 168 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 169 & / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 170 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 171 & / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 172 zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) / ze3ua 173 zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 174 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 175 zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) 176 zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 177 END DO 178 END DO 179 END DO 164 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 165 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) & 166 & + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 167 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 168 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 169 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 170 & / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 171 zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) / ze3ua 172 zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 173 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 174 zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) 175 zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 176 END_3D 180 177 CASE DEFAULT ! iso-level lateral mixing 181 DO jk = 1, jpkm1 182 DO jj = 2, jpjm1 183 DO ji = fs_2, fs_jpim1 ! vector opt. 184 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 185 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 186 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 187 zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) / ze3ua 188 zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 189 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 190 zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) 191 zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 192 END DO 193 END DO 194 END DO 178 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 179 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) & ! after scale factor at U-point 180 & + r_vvl * e3u(ji,jj,jk,Kaa) 181 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) & 182 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 183 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) & 184 & / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 185 zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) / ze3ua 186 zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 187 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 188 zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) 189 zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 190 END_3D 195 191 END SELECT 196 DO jj = 2, jpjm1 !* Surface boundary conditions197 DO ji = fs_2, fs_jpim1 ! vector opt.198 zwi(ji,jj,1) = 0._wp199 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)200 zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) / ( ze3ua * e3uw_n(ji,jj,2) ) * wumask(ji,jj,2)201 zWus = ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) / ze3ua202 zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp )203 zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ))204 END DO205 END DO192 DO_2D( 0, 0, 0, 0 ) 193 zwi(ji,jj,1) = 0._wp 194 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) & 195 & + r_vvl * e3u(ji,jj,1,Kaa) 196 zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) & 197 & / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2) 198 zWus = ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) / ze3ua 199 zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp ) 200 zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) ) 201 END_2D 206 202 ELSE 207 203 SELECT CASE( nldf_dyn ) 208 204 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 209 DO jk = 1, jpkm1 210 DO jj = 2, jpjm1 211 DO ji = fs_2, fs_jpim1 ! vector opt. 212 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 213 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 214 & / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 215 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 216 & / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 217 zwi(ji,jj,jk) = zzwi 218 zws(ji,jj,jk) = zzws 219 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 220 END DO 221 END DO 222 END DO 205 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 206 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) & 207 & + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 208 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 209 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 210 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 211 & / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 212 zwi(ji,jj,jk) = zzwi 213 zws(ji,jj,jk) = zzws 214 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 215 END_3D 223 216 CASE DEFAULT ! iso-level lateral mixing 224 DO jk = 1, jpkm1 225 DO jj = 2, jpjm1 226 DO ji = fs_2, fs_jpim1 ! vector opt. 227 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 228 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 229 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 230 zwi(ji,jj,jk) = zzwi 231 zws(ji,jj,jk) = zzws 232 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 233 END DO 234 END DO 235 END DO 217 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 218 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) & 219 & + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 220 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) & 221 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 222 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) & 223 & / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 224 zwi(ji,jj,jk) = zzwi 225 zws(ji,jj,jk) = zzws 226 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 227 END_3D 236 228 END SELECT 237 DO jj = 2, jpjm1 !* Surface boundary conditions 238 DO ji = fs_2, fs_jpim1 ! vector opt. 239 zwi(ji,jj,1) = 0._wp 240 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 241 END DO 242 END DO 229 DO_2D( 0, 0, 0, 0 ) 230 zwi(ji,jj,1) = 0._wp 231 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 232 END_2D 243 233 ENDIF 244 234 ! … … 251 241 ! 252 242 IF ( ln_drgimp ) THEN ! implicit bottom friction 253 DO jj = 2, jpjm1 254 DO ji = 2, jpim1 255 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 256 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point 257 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 258 END DO 259 END DO 243 DO_2D( 0, 0, 0, 0 ) 244 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 245 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) & 246 & + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point 247 zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 248 END_2D 260 249 IF ( ln_isfcav ) THEN ! top friction (always implicit) 261 DO jj = 2, jpjm1 262 DO ji = 2, jpim1 263 !!gm top Cd is masked (=0 outside cavities) no need of test on mik>=2 ==>> it has been suppressed 264 iku = miku(ji,jj) ! ocean top level at u- and v-points 265 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point 266 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 267 END DO 268 END DO 250 DO_2D( 0, 0, 0, 0 ) 251 !!gm top Cd is masked (=0 outside cavities) no need of test on mik>=2 ==>> it has been suppressed 252 iku = miku(ji,jj) ! ocean top level at u- and v-points 253 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) & 254 & + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point 255 zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 256 END_2D 269 257 END IF 270 258 ENDIF … … 282 270 ! m is decomposed in the product of an upper and a lower triangular matrix 283 271 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 284 ! The solution (the after velocity) is in ua272 ! The solution (the after velocity) is in puu(:,:,:,Kaa) 285 273 !----------------------------------------------------------------------- 286 274 ! 287 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 288 DO jj = 2, jpjm1 289 DO ji = fs_2, fs_jpim1 ! vector opt. 290 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 291 END DO 292 END DO 293 END DO 294 ! 295 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 296 DO ji = fs_2, fs_jpim1 ! vector opt. 297 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 298 ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 299 & / ( ze3ua * rau0 ) * umask(ji,jj,1) 300 END DO 301 END DO 302 DO jk = 2, jpkm1 303 DO jj = 2, jpjm1 304 DO ji = fs_2, fs_jpim1 305 ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 306 END DO 307 END DO 308 END DO 309 ! 310 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==! 311 DO ji = fs_2, fs_jpim1 ! vector opt. 312 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 313 END DO 314 END DO 315 DO jk = jpk-2, 1, -1 316 DO jj = 2, jpjm1 317 DO ji = fs_2, fs_jpim1 318 ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 319 END DO 320 END DO 321 END DO 275 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 276 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 277 END_3D 278 ! 279 DO_2D( 0, 0, 0, 0 ) 280 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) & 281 & + r_vvl * e3u(ji,jj,1,Kaa) 282 puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 283 & / ( ze3ua * rho0 ) * umask(ji,jj,1) 284 END_2D 285 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 286 puu(ji,jj,jk,Kaa) = puu(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * puu(ji,jj,jk-1,Kaa) 287 END_3D 288 ! 289 DO_2D( 0, 0, 0, 0 ) 290 puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1) 291 END_2D 292 DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 ) 293 puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jj,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk) 294 END_3D 322 295 ! 323 296 ! !== Vertical diffusion on v ==! 324 297 ! 325 298 ! !* Matrix construction 326 zdt = r 2dt * 0.5299 zdt = rDt * 0.5 327 300 IF( ln_zad_Aimp ) THEN !! 328 301 SELECT CASE( nldf_dyn ) 329 302 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzv) 330 DO jk = 1, jpkm1 331 DO jj = 2, jpjm1 332 DO ji = fs_2, fs_jpim1 ! vector opt. 333 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point 334 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 335 & / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 336 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 337 & / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 338 zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) / ze3va 339 zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 340 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp ) 341 zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp ) 342 zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 343 END DO 344 END DO 345 END DO 303 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 304 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) & 305 & + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 306 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 307 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 308 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 309 & / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 310 zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) / ze3va 311 zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 312 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp ) 313 zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp ) 314 zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 315 END_3D 346 316 CASE DEFAULT ! iso-level lateral mixing 347 DO jk = 1, jpkm1 348 DO jj = 2, jpjm1 349 DO ji = fs_2, fs_jpim1 ! vector opt. 350 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point 351 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 352 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 353 zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) / ze3va 354 zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 355 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp ) 356 zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp ) 357 zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 358 END DO 359 END DO 360 END DO 317 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 318 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) & 319 & + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 320 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) & 321 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 322 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) & 323 & / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 324 zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) / ze3va 325 zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 326 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp ) 327 zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp ) 328 zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 329 END_3D 361 330 END SELECT 362 DO jj = 2, jpjm1 !* Surface boundary conditions363 DO ji = fs_2, fs_jpim1 ! vector opt.364 zwi(ji,jj,1) = 0._wp365 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)366 zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw_n(ji,jj,2) ) * wvmask(ji,jj,2)367 zWvs = ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) / ze3va368 zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp )369 zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ))370 END DO371 END DO331 DO_2D( 0, 0, 0, 0 ) 332 zwi(ji,jj,1) = 0._wp 333 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) & 334 & + r_vvl * e3v(ji,jj,1,Kaa) 335 zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) & 336 & / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2) 337 zWvs = ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) / ze3va 338 zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp ) 339 zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) ) 340 END_2D 372 341 ELSE 373 342 SELECT CASE( nldf_dyn ) 374 343 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 375 DO jk = 1, jpkm1 376 DO jj = 2, jpjm1 377 DO ji = fs_2, fs_jpim1 ! vector opt. 378 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point 379 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 380 & / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 381 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 382 & / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 383 zwi(ji,jj,jk) = zzwi 384 zws(ji,jj,jk) = zzws 385 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 386 END DO 387 END DO 388 END DO 344 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 345 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) & 346 & + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 347 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 348 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 349 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 350 & / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 351 zwi(ji,jj,jk) = zzwi 352 zws(ji,jj,jk) = zzws 353 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 354 END_3D 389 355 CASE DEFAULT ! iso-level lateral mixing 390 DO jk = 1, jpkm1 391 DO jj = 2, jpjm1 392 DO ji = fs_2, fs_jpim1 ! vector opt. 393 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point 394 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 395 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 396 zwi(ji,jj,jk) = zzwi 397 zws(ji,jj,jk) = zzws 398 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 399 END DO 400 END DO 401 END DO 356 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 357 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) & 358 & + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 359 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) & 360 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 361 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) & 362 & / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 363 zwi(ji,jj,jk) = zzwi 364 zws(ji,jj,jk) = zzws 365 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 366 END_3D 402 367 END SELECT 403 DO jj = 2, jpjm1 !* Surface boundary conditions 404 DO ji = fs_2, fs_jpim1 ! vector opt. 405 zwi(ji,jj,1) = 0._wp 406 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 407 END DO 408 END DO 368 DO_2D( 0, 0, 0, 0 ) 369 zwi(ji,jj,1) = 0._wp 370 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 371 END_2D 409 372 ENDIF 410 373 ! … … 416 379 ! 417 380 IF( ln_drgimp ) THEN 418 DO jj = 2, jpjm1 419 DO ji = 2, jpim1 420 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 421 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point 422 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va 423 END DO 424 END DO 381 DO_2D( 0, 0, 0, 0 ) 382 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 383 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) & 384 & + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point 385 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va 386 END_2D 425 387 IF ( ln_isfcav ) THEN 426 DO jj = 2, jpjm1 427 DO ji = 2, jpim1 428 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 429 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point 430 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va 431 END DO 432 END DO 388 DO_2D( 0, 0, 0, 0 ) 389 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 390 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) & 391 & + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point 392 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va 393 END_2D 433 394 ENDIF 434 395 ENDIF … … 449 410 !----------------------------------------------------------------------- 450 411 ! 451 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 452 DO jj = 2, jpjm1 453 DO ji = fs_2, fs_jpim1 ! vector opt. 454 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 455 END DO 456 END DO 457 END DO 458 ! 459 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 460 DO ji = fs_2, fs_jpim1 ! vector opt. 461 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 462 va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 463 & / ( ze3va * rau0 ) * vmask(ji,jj,1) 464 END DO 465 END DO 466 DO jk = 2, jpkm1 467 DO jj = 2, jpjm1 468 DO ji = fs_2, fs_jpim1 ! vector opt. 469 va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 470 END DO 471 END DO 472 END DO 473 ! 474 DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==! 475 DO ji = fs_2, fs_jpim1 ! vector opt. 476 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 477 END DO 478 END DO 479 DO jk = jpk-2, 1, -1 480 DO jj = 2, jpjm1 481 DO ji = fs_2, fs_jpim1 482 va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 483 END DO 484 END DO 485 END DO 412 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 413 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 414 END_3D 415 ! 416 DO_2D( 0, 0, 0, 0 ) 417 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) & 418 & + r_vvl * e3v(ji,jj,1,Kaa) 419 pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 420 & / ( ze3va * rho0 ) * vmask(ji,jj,1) 421 END_2D 422 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 423 pvv(ji,jj,jk,Kaa) = pvv(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * pvv(ji,jj,jk-1,Kaa) 424 END_3D 425 ! 426 DO_2D( 0, 0, 0, 0 ) 427 pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1) 428 END_2D 429 DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 ) 430 pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jj,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk) 431 END_3D 486 432 ! 487 433 IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics 488 ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:)489 ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:)490 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )434 ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) / rDt - ztrdu(:,:,:) 435 ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) / rDt - ztrdv(:,:,:) 436 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt, Kmm ) 491 437 DEALLOCATE( ztrdu, ztrdv ) 492 438 ENDIF 493 439 ! ! print mean trends (used for debugging) 494 IF( ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf - Ua: ', mask1=umask, &495 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )440 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' zdf - Ua: ', mask1=umask, & 441 & tab3d_2=pvv(:,:,:,Kaa), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 496 442 ! 497 443 IF( ln_timing ) CALL timing_stop('dyn_zdf')
Note: See TracChangeset
for help on using the changeset viewer.