Changeset 11822 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN
- Timestamp:
- 2019-10-29T11:41:36+01:00 (5 years ago)
- Location:
- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN
- Files:
-
- 9 edited
-
dynadv.F90 (modified) (1 diff)
-
dynhpg.F90 (modified) (7 diffs)
-
dynkeg.F90 (modified) (4 diffs)
-
dynspg.F90 (modified) (1 diff)
-
dynspg_ts.F90 (modified) (28 diffs)
-
dynvor.F90 (modified) (1 diff)
-
dynzdf.F90 (modified) (6 diffs)
-
sshwzv.F90 (modified) (5 diffs)
-
wet_dry.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv.F90
r10893 r11822 108 108 REWIND( numnam_ref ) ! Namelist namdyn_adv in reference namelist : Momentum advection scheme 109 109 READ ( numnam_ref, namdyn_adv, IOSTAT = ios, ERR = 901) 110 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_adv in reference namelist' , lwp)110 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_adv in reference namelist' ) 111 111 REWIND( numnam_cfg ) ! Namelist namdyn_adv in configuration namelist : Momentum advection scheme 112 112 READ ( numnam_cfg, namdyn_adv, IOSTAT = ios, ERR = 902 ) 113 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_adv in configuration namelist' , lwp)113 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_adv in configuration namelist' ) 114 114 IF(lwm) WRITE ( numond, namdyn_adv ) 115 115 -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynhpg.F90
r10946 r11822 37 37 USE trd_oce ! trends: ocean variables 38 38 USE trddyn ! trend manager: dynamics 39 !jcUSE zpshde ! partial step: hor. derivative (zps_hde routine)39 USE zpshde ! partial step: hor. derivative (zps_hde routine) 40 40 ! 41 41 USE in_out_manager ! I/O manager … … 157 157 REWIND( numnam_ref ) ! Namelist namdyn_hpg in reference namelist : Hydrostatic pressure gradient 158 158 READ ( numnam_ref, namdyn_hpg, IOSTAT = ios, ERR = 901) 159 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in reference namelist' , lwp)159 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in reference namelist' ) 160 160 ! 161 161 REWIND( numnam_cfg ) ! Namelist namdyn_hpg in configuration namelist : Hydrostatic pressure gradient 162 162 READ ( numnam_cfg, namdyn_hpg, IOSTAT = ios, ERR = 902 ) 163 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist' , lwp)163 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist' ) 164 164 IF(lwm) WRITE ( numond, namdyn_hpg ) 165 165 ! … … 347 347 REAL(wp) :: zcoef0, zcoef1, zcoef2, zcoef3 ! temporary scalars 348 348 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj 349 REAL(wp), DIMENSION(jpi,jpj) :: zgtsu, zgtsv, zgru, zgrv 349 350 !!---------------------------------------------------------------------- 350 351 ! … … 355 356 ENDIF 356 357 357 ! Partial steps: bottom beforehorizontal gradient of t, s, rd at the last ocean level358 !jc CALL zps_hde ( kt, jpts, ts(:,:,:,:,Kmm), gtsu, gtsv, rhd, gru ,grv )358 ! Partial steps: Compute NOW horizontal gradient of t, s, rd at the last ocean level 359 CALL zps_hde( kt, Kmm, jpts, ts(:,:,:,:,Kmm), zgtsu, zgtsv, rhd, zgru , zgrv ) 359 360 360 361 ! Local constant initialization … … 394 395 END DO 395 396 396 ! partial steps correction at the last level (use gru &grv computed in zpshde.F90)397 ! partial steps correction at the last level (use zgru & zgrv computed in zpshde.F90) 397 398 DO jj = 2, jpjm1 398 399 DO ji = 2, jpim1 … … 404 405 puu (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) - zhpi(ji,jj,iku) ! subtract old value 405 406 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) & ! compute the new one 406 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj)407 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + zgru(ji,jj) ) * r1_e1u(ji,jj) 407 408 puu (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend 408 409 ENDIF … … 410 411 pvv (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) - zhpj(ji,jj,ikv) ! subtract old value 411 412 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) & ! compute the new one 412 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj)413 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + zgrv(ji,jj) ) * r1_e2v(ji,jj) 413 414 pvv (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend 414 415 ENDIF -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynkeg.F90
r10946 r11822 76 76 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 77 77 ! 78 INTEGER :: ji, jj, jk, jb ! dummy loop indices 79 INTEGER :: ii, ifu, ib_bdy ! local integers 80 INTEGER :: ij, ifv, igrd ! - - 81 REAL(wp) :: zu, zv ! local scalars 78 INTEGER :: ji, jj, jk ! dummy loop indices 79 REAL(wp) :: zu, zv ! local scalars 82 80 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhke 83 81 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv … … 99 97 100 98 zhke(:,:,jpk) = 0._wp 101 102 IF (ln_bdy) THEN103 ! Maria Luneva & Fred Wobus: July-2016104 ! compensate for lack of turbulent kinetic energy on liquid bdy points105 DO ib_bdy = 1, nb_bdy106 IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN107 igrd = 2 ! Copying normal velocity into points outside bdy108 DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)109 DO jk = 1, jpkm1110 ii = idx_bdy(ib_bdy)%nbi(jb,igrd)111 ij = idx_bdy(ib_bdy)%nbj(jb,igrd)112 ifu = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) )113 puu(ii-ifu,ij,jk,Kmm) = puu(ii,ij,jk,Kmm) * umask(ii,ij,jk)114 END DO115 END DO116 !117 igrd = 3 ! Copying normal velocity into points outside bdy118 DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)119 DO jk = 1, jpkm1120 ii = idx_bdy(ib_bdy)%nbi(jb,igrd)121 ij = idx_bdy(ib_bdy)%nbj(jb,igrd)122 ifv = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) )123 pvv(ii,ij-ifv,jk,Kmm) = pvv(ii,ij,jk,Kmm) * vmask(ii,ij,jk)124 END DO125 END DO126 ENDIF127 ENDDO128 ENDIF129 99 130 100 SELECT CASE ( kscheme ) !== Horizontal kinetic energy at T-point ==! … … 142 112 END DO 143 113 END DO 144 !145 114 CASE ( nkeg_HW ) !-- Hollingsworth scheme --! 146 115 DO jk = 1, jpkm1 … … 162 131 CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. ) 163 132 ! 164 END SELECT 165 166 IF (ln_bdy) THEN 167 ! restore velocity masks at points outside boundary 168 puu(:,:,:,Kmm) = puu(:,:,:,Kmm) * umask(:,:,:) 169 pvv(:,:,:,Kmm) = pvv(:,:,:,Kmm) * vmask(:,:,:) 170 ENDIF 171 133 END SELECT 172 134 ! 173 135 DO jk = 1, jpkm1 !== grad( KE ) added to the general momentum trends ==! -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynspg.F90
r10946 r11822 205 205 REWIND( numnam_ref ) ! Namelist namdyn_spg in reference namelist : Free surface 206 206 READ ( numnam_ref, namdyn_spg, IOSTAT = ios, ERR = 901) 207 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_spg in reference namelist' , lwp)207 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_spg in reference namelist' ) 208 208 ! 209 209 REWIND( numnam_cfg ) ! Namelist namdyn_spg in configuration namelist : Free surface 210 210 READ ( numnam_cfg, namdyn_spg, IOSTAT = ios, ERR = 902 ) 211 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_spg in configuration namelist' , lwp)211 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_spg in configuration namelist' ) 212 212 IF(lwm) WRITE ( numond, namdyn_spg ) 213 213 ! -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynspg_ts.F90
r11480 r11822 64 64 USE diatmb ! Top,middle,bottom output 65 65 66 USE iom ! to remove 67 66 68 IMPLICIT NONE 67 69 PRIVATE … … 104 106 ! 105 107 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) ) 106 !107 108 IF( ln_dynvor_een .OR. ln_dynvor_eeT ) & 108 & ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 109 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) ) 109 & ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , ftsw(jpi,jpj) , ftse(jpi,jpj), STAT=ierr(2) ) 110 110 ! 111 111 ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(3) ) … … 152 152 LOGICAL :: ll_fw_start ! =T : forward integration 153 153 LOGICAL :: ll_init ! =T : special startup of 2d equations 154 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D 155 INTEGER :: ikbu, iktu, noffset ! local integers 156 INTEGER :: ikbv, iktv ! - - 157 REAL(wp) :: r1_2dt_b, z2dt_bf ! local scalars 158 REAL(wp) :: zx1, zx2, zu_spg, zhura, z1_hu ! - - 159 REAL(wp) :: zy1, zy2, zv_spg, zhvra, z1_hv ! - - 154 INTEGER :: noffset ! local integers : time offset for bdy update 155 REAL(wp) :: r1_2dt_b, z1_hu, z1_hv ! local scalars 160 156 REAL(wp) :: za0, za1, za2, za3 ! - - 161 REAL(wp) :: zmdi, zztmp , z1_ht ! - - 162 REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf 163 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc 164 REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zhdiv 165 REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e, zhtp2_e 166 REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zhvst_e 157 REAL(wp) :: zmdi, zztmp, zldg ! - - 158 REAL(wp) :: zhu_bck, zhv_bck, zhdiv ! - - 159 REAL(wp) :: zun_save, zvn_save ! - - 160 REAL(wp), DIMENSION(jpi,jpj) :: zu_trd, zu_frc, zu_spg, zssh_frc 161 REAL(wp), DIMENSION(jpi,jpj) :: zv_trd, zv_frc, zv_spg 162 REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhtp2_e 163 REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zsshp2_e 167 164 REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points 165 REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV ! fluxes 168 166 ! 169 167 REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True. … … 185 183 zwdramp = r_rn_wdmin1 ! simplest ramp 186 184 ! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp 187 ! ! reciprocal of baroclinic time step 188 IF( kt == nit000 .AND. neuler == 0 ) THEN ; z2dt_bf = rdt 189 ELSE ; z2dt_bf = 2.0_wp * rdt 190 ENDIF 191 r1_2dt_b = 1.0_wp / z2dt_bf 185 ! ! inverse of baroclinic time step 186 IF( kt == nit000 .AND. neuler == 0 ) THEN ; r1_2dt_b = 1._wp / ( rdt ) 187 ELSE ; r1_2dt_b = 1._wp / ( 2._wp * rdt ) 188 ENDIF 192 189 ! 193 190 ll_init = ln_bt_av ! if no time averaging, then no specific restart … … 213 210 ll_fw_start =.FALSE. 214 211 ENDIF 215 ! 216 ! Set averaging weights and cycle length: 212 ! ! Set averaging weights and cycle length: 217 213 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 218 214 ! 219 ENDIF220 !221 IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities)222 DO jj = 2, jpjm1223 DO ji = fs_2, fs_jpim1 ! vector opt.224 zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )225 zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )226 END DO227 END DO228 ELSE ! bottom friction only229 DO jj = 2, jpjm1230 DO ji = fs_2, fs_jpim1 ! vector opt.231 zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )232 zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )233 END DO234 END DO235 ENDIF236 !237 ! Set arrays to remove/compute coriolis trend.238 ! Do it once at kt=nit000 if volume is fixed, else at each long time step.239 ! Note that these arrays are also used during barotropic loop. These are however frozen240 ! although they should be updated in the variable volume case. Not a big approximation.241 ! To remove this approximation, copy lines below inside barotropic loop242 ! and update depths at T-F points (ht and zhf resp.) at each barotropic time step243 !244 IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN245 !246 SELECT CASE( nvor_scheme )247 CASE( np_EEN ) != EEN scheme using e3f (energy & enstrophy scheme)248 SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point249 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)250 DO jj = 1, jpjm1251 DO ji = 1, jpim1252 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + &253 & ht(ji ,jj ) + ht(ji+1,jj ) ) * 0.25_wp254 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)255 END DO256 END DO257 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)258 DO jj = 1, jpjm1259 DO ji = 1, jpim1260 zwz(ji,jj) = ( ht (ji ,jj+1) + ht (ji+1,jj+1) &261 & + ht (ji ,jj ) + ht (ji+1,jj ) ) &262 & / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) &263 & + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) )264 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)265 END DO266 END DO267 END SELECT268 CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )269 !270 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp271 DO jj = 2, jpj272 DO ji = 2, jpi273 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)274 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )275 ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)276 ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )277 END DO278 END DO279 !280 CASE( np_EET ) != EEN scheme using e3t (energy conserving scheme)281 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp282 DO jj = 2, jpj283 DO ji = 2, jpi284 z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) )285 ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht286 ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht287 ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht288 ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht289 END DO290 END DO291 !292 CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT !293 !294 zwz(:,:) = 0._wp295 zhf(:,:) = 0._wp296 297 !!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed298 !!gm A priori a better value should be something like :299 !!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)300 !!gm divided by the sum of the corresponding mask301 !!gm302 !!303 IF( .NOT.ln_sco ) THEN304 305 !!gm agree the JC comment : this should be done in a much clear way306 307 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case308 ! Set it to zero for the time being309 ! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level310 ! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth311 ! ENDIF312 ! zhf(:,:) = gdepw_0(:,:,jk+1)313 !314 ELSE315 !316 !zhf(:,:) = hbatf(:,:)317 DO jj = 1, jpjm1318 DO ji = 1, jpim1319 zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) &320 & + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) &321 & / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) &322 & + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp )323 END DO324 END DO325 ENDIF326 !327 DO jj = 1, jpjm1328 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))329 END DO330 !331 DO jk = 1, jpkm1332 DO jj = 1, jpjm1333 zhf(:,jj) = zhf(:,jj) + e3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)334 END DO335 END DO336 CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )337 ! JC: TBC. hf should be greater than 0338 DO jj = 1, jpj339 DO ji = 1, jpi340 IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) ! zhf is actually hf here but it saves an array341 END DO342 END DO343 zwz(:,:) = ff_f(:,:) * zwz(:,:)344 END SELECT345 215 ENDIF 346 216 ! … … 351 221 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 352 222 ENDIF 223 ! 353 224 354 225 ! ----------------------------------------------------------------------------- … … 357 228 ! 358 229 ! 359 ! !* e3*d/dt(Ua) (Vertically integrated) 360 ! ! -------------------------------------------------- 361 zu_frc(:,:) = 0._wp 362 zv_frc(:,:) = 0._wp 363 ! 364 DO jk = 1, jpkm1 365 zu_frc(:,:) = zu_frc(:,:) + e3u(:,:,jk,Kmm) * puu(:,:,jk,Krhs) * umask(:,:,jk) 366 zv_frc(:,:) = zv_frc(:,:) + e3v(:,:,jk,Kmm) * pvv(:,:,jk,Krhs) * vmask(:,:,jk) 367 END DO 368 ! 369 zu_frc(:,:) = zu_frc(:,:) * r1_hu(:,:,Kmm) 370 zv_frc(:,:) = zv_frc(:,:) * r1_hv(:,:,Kmm) 371 ! 372 ! 373 ! !* baroclinic momentum trend (remove the vertical mean trend) 374 DO jk = 1, jpkm1 ! ----------------------------------------------------------- 375 DO jj = 2, jpjm1 376 DO ji = fs_2, fs_jpim1 ! vector opt. 377 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - zu_frc(ji,jj) * umask(ji,jj,jk) 378 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - zv_frc(ji,jj) * vmask(ji,jj,jk) 379 END DO 380 END DO 230 ! != zu_frc = 1/H e3*d/dt(Ua) =! (Vertical mean of Ua, the 3D trends) 231 ! ! --------------------------- ! 232 zu_frc(:,:) = SUM( e3u(:,:,:,Kmm) * uu(:,:,:,Krhs) * umask(:,:,:) , DIM=3 ) * r1_hu(:,:,Kmm) 233 zv_frc(:,:) = SUM( e3v(:,:,:,Kmm) * vv(:,:,:,Krhs) * vmask(:,:,:) , DIM=3 ) * r1_hv(:,:,Kmm) 234 ! 235 ! 236 ! != U(Krhs) => baroclinic trend =! (remove its vertical mean) 237 DO jk = 1, jpkm1 ! ----------------------------- ! 238 uu(:,:,jk,Krhs) = ( uu(:,:,jk,Krhs) - zu_frc(:,:) ) * umask(:,:,jk) 239 vv(:,:,jk,Krhs) = ( vv(:,:,jk,Krhs) - zv_frc(:,:) ) * vmask(:,:,jk) 381 240 END DO 382 241 … … 384 243 !!gm Is it correct to do so ? I think so... 385 244 386 387 ! !* barotropic Coriolis trends (vorticity scheme dependent) 388 ! ! -------------------------------------------------------- 389 ! 390 zwx(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:) ! now fluxes 391 zwy(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:) 392 ! 393 SELECT CASE( nvor_scheme ) 394 CASE( np_ENT ) ! enstrophy conserving scheme (f-point) 395 DO jj = 2, jpjm1 396 DO ji = 2, jpim1 ! vector opt. 397 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * r1_hu(ji,jj,Kmm) & 398 & * ( e1e2t(ji+1,jj)*ht(ji+1,jj)*ff_t(ji+1,jj) * ( pvv_b(ji+1,jj,Kmm) + pvv_b(ji+1,jj-1,Kmm) ) & 399 & + e1e2t(ji ,jj)*ht(ji ,jj)*ff_t(ji ,jj) * ( pvv_b(ji ,jj,Kmm) + pvv_b(ji ,jj-1,Kmm) ) ) 400 ! 401 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * r1_hv(ji,jj,Kmm) & 402 & * ( e1e2t(ji,jj+1)*ht(ji,jj+1)*ff_t(ji,jj+1) * ( puu_b(ji,jj+1,Kmm) + puu_b(ji-1,jj+1,Kmm) ) & 403 & + e1e2t(ji,jj )*ht(ji,jj )*ff_t(ji,jj ) * ( puu_b(ji,jj ,Kmm) + puu_b(ji-1,jj ,Kmm) ) ) 404 END DO 405 END DO 406 ! 407 CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX 408 DO jj = 2, jpjm1 409 DO ji = fs_2, fs_jpim1 ! vector opt. 410 zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 411 zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 412 zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 413 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 414 ! energy conserving formulation for planetary vorticity term 415 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 416 zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 417 END DO 418 END DO 419 ! 420 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 421 DO jj = 2, jpjm1 422 DO ji = fs_2, fs_jpim1 ! vector opt. 423 zy1 = r1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 424 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 425 zx1 = - r1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 426 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 427 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 428 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 429 END DO 430 END DO 431 ! 432 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 433 DO jj = 2, jpjm1 434 DO ji = fs_2, fs_jpim1 ! vector opt. 435 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 436 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 437 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 438 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 439 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 440 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 441 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 442 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 443 END DO 444 END DO 445 ! 446 END SELECT 447 ! 448 ! !* Right-Hand-Side of the barotropic momentum equation 449 ! ! ---------------------------------------------------- 450 IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient 451 IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters 245 ! != remove 2D Coriolis and pressure gradient trends =! 246 ! ! ------------------------------------------------- ! 247 ! 248 IF( kt == nit000 .OR. .NOT. ln_linssh ) CALL dyn_cor_2D_init( Kmm ) ! Set zwz, the barotropic Coriolis force coefficient 249 ! ! recompute zwz = f/depth at every time step for (.NOT.ln_linssh) as the water colomn height changes 250 ! 251 ! !* 2D Coriolis trends 252 zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:) ! now fluxes 253 zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:) ! NB: FULL domain : put a value in last row and column 254 ! 255 CALL dyn_cor_2d( hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV, & ! <<== in 256 & zu_trd, zv_trd ) ! ==>> out 257 ! 258 IF( .NOT.ln_linssh ) THEN !* surface pressure gradient (variable volume only) 259 ! 260 IF( ln_wd_il ) THEN ! W/D : limiter applied to spgspg 261 CALL wad_spg( pssh(:,:,Kmm), zcpx, zcpy ) ! Calculating W/D gravity filters, zcpx and zcpy 452 262 DO jj = 2, jpjm1 453 DO ji = 2, jpim1 454 ll_tmp1 = MIN( pssh(ji,jj,Kmm) , pssh(ji+1,jj,Kmm) ) > & 455 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 456 & MAX( pssh(ji,jj,Kmm) + ht_0(ji,jj) , pssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) & 457 & > rn_wdmin1 + rn_wdmin2 458 ll_tmp2 = ( ABS( pssh(ji+1,jj,Kmm) - pssh(ji ,jj,Kmm)) > 1.E-12 ).AND.( & 459 & MAX( pssh(ji,jj,Kmm) , pssh(ji+1,jj,Kmm) ) > & 460 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 461 IF(ll_tmp1) THEN 462 zcpx(ji,jj) = 1.0_wp 463 ELSEIF(ll_tmp2) THEN 464 ! no worries about pssh(ji+1,jj,Kmm) - pssh(ji ,jj,Kmm) = 0, it won't happen ! here 465 zcpx(ji,jj) = ABS( (pssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - pssh(ji,jj,Kmm) - ht_0(ji,jj)) & 466 & / (pssh(ji+1,jj,Kmm) - pssh(ji ,jj,Kmm)) ) 467 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 468 ELSE 469 zcpx(ji,jj) = 0._wp 470 ENDIF 471 ! 472 ll_tmp1 = MIN( pssh(ji,jj,Kmm) , pssh(ji,jj+1,Kmm) ) > & 473 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 474 & MAX( pssh(ji,jj,Kmm) + ht_0(ji,jj) , pssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) & 475 & > rn_wdmin1 + rn_wdmin2 476 ll_tmp2 = ( ABS( pssh(ji,jj,Kmm) - pssh(ji,jj+1,Kmm)) > 1.E-12 ).AND.( & 477 & MAX( pssh(ji,jj,Kmm) , pssh(ji,jj+1,Kmm) ) > & 478 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 479 480 IF(ll_tmp1) THEN 481 zcpy(ji,jj) = 1.0_wp 482 ELSE IF(ll_tmp2) THEN 483 ! no worries about pssh(ji,jj+1,Kmm) - pssh(ji,jj ,Kmm) = 0, it won't happen ! here 484 zcpy(ji,jj) = ABS( (pssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - pssh(ji,jj,Kmm) - ht_0(ji,jj)) & 485 & / (pssh(ji,jj+1,Kmm) - pssh(ji,jj ,Kmm)) ) 486 zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) ) 487 ELSE 488 zcpy(ji,jj) = 0._wp 489 ENDIF 490 END DO 491 END DO 492 ! 493 DO jj = 2, jpjm1 494 DO ji = 2, jpim1 263 DO ji = 2, jpim1 ! SPG with the application of W/D gravity filters 495 264 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj ,Kmm) - pssh(ji ,jj ,Kmm) ) & 496 265 & * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth … … 499 268 END DO 500 269 END DO 501 ! 502 ELSE 503 ! 270 ELSE ! now suface pressure gradient 504 271 DO jj = 2, jpjm1 505 272 DO ji = fs_2, fs_jpim1 ! vector opt. … … 519 286 END DO 520 287 ! 521 ! ! Add bottom stress contribution from baroclinic velocities: 522 IF (ln_bt_fw) THEN 523 DO jj = 2, jpjm1 524 DO ji = fs_2, fs_jpim1 ! vector opt. 525 ikbu = mbku(ji,jj) 526 ikbv = mbkv(ji,jj) 527 zwx(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm) ! NOW bottom baroclinic velocities 528 zwy(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm) 529 END DO 530 END DO 531 ELSE 532 DO jj = 2, jpjm1 533 DO ji = fs_2, fs_jpim1 ! vector opt. 534 ikbu = mbku(ji,jj) 535 ikbv = mbkv(ji,jj) 536 zwx(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb) ! BEFORE bottom baroclinic velocities 537 zwy(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb) 538 END DO 539 END DO 540 ENDIF 541 ! 542 ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 543 IF( ln_wd_il ) THEN 544 zztmp = -1._wp / rdtbt 545 DO jj = 2, jpjm1 546 DO ji = fs_2, fs_jpim1 ! vector opt. 547 zu_frc(ji,jj) = zu_frc(ji,jj) + & 548 & MAX(r1_hu(ji,jj,Kmm) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ), zztmp ) * zwx(ji,jj) * wdrampu(ji,jj) 549 zv_frc(ji,jj) = zv_frc(ji,jj) + & 550 & MAX(r1_hv(ji,jj,Kmm) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ), zztmp ) * zwy(ji,jj) * wdrampv(ji,jj) 551 END DO 552 END DO 553 ELSE 554 DO jj = 2, jpjm1 555 DO ji = fs_2, fs_jpim1 ! vector opt. 556 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zwx(ji,jj) 557 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zwy(ji,jj) 558 END DO 559 END DO 560 END IF 561 ! 562 IF( ln_isfcav ) THEN ! Add TOP stress contribution from baroclinic velocities: 563 IF( ln_bt_fw ) THEN 564 DO jj = 2, jpjm1 288 ! != Add bottom stress contribution from baroclinic velocities =! 289 ! ! ----------------------------------------------------------- ! 290 CALL dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, zu_frc, zv_frc, zCdU_u, zCdU_v ) ! also provide the barotropic drag coefficients 291 ! != Add atmospheric pressure forcing =! 292 ! ! ---------------------------------- ! 293 IF( ln_apr_dyn ) THEN 294 IF( ln_bt_fw ) THEN ! FORWARD integration: use kt+1/2 pressure (NOW+1/2) 295 DO jj = 2, jpjm1 565 296 DO ji = fs_2, fs_jpim1 ! vector opt. 566 iktu = miku(ji,jj) 567 iktv = mikv(ji,jj) 568 zwx(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm) ! NOW top baroclinic velocities 569 zwy(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm) 570 END DO 571 END DO 572 ELSE 573 DO jj = 2, jpjm1 297 zu_frc(ji,jj) = zu_frc(ji,jj) + grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 298 zv_frc(ji,jj) = zv_frc(ji,jj) + grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) 299 END DO 300 END DO 301 ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW) 302 zztmp = grav * r1_2 303 DO jj = 2, jpjm1 574 304 DO ji = fs_2, fs_jpim1 ! vector opt. 575 iktu = miku(ji,jj) 576 iktv = mikv(ji,jj) 577 zwx(ji,jj) = puu(ji,jj,iktu,Kbb) - puu_b(ji,jj,Kbb) ! BEFORE top baroclinic velocities 578 zwy(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb) 579 END DO 580 END DO 581 ENDIF 582 ! 583 ! Note that the "unclipped" top friction parameter is used even with explicit drag 584 DO jj = 2, jpjm1 585 DO ji = fs_2, fs_jpim1 ! vector opt. 586 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zwx(ji,jj) 587 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zwy(ji,jj) 588 END DO 589 END DO 590 ENDIF 591 ! 305 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 306 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) 307 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) & 308 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj) 309 END DO 310 END DO 311 ENDIF 312 ENDIF 313 ! 314 ! != Add atmospheric pressure forcing =! 315 ! ! ---------------------------------- ! 592 316 IF( ln_bt_fw ) THEN ! Add wind forcing 593 317 DO jj = 2, jpjm1 … … 607 331 ENDIF 608 332 ! 609 IF( ln_apr_dyn ) THEN ! Add atm pressure forcing 610 IF( ln_bt_fw ) THEN 611 DO jj = 2, jpjm1 612 DO ji = fs_2, fs_jpim1 ! vector opt. 613 zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 614 zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) 615 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 616 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg 617 END DO 618 END DO 619 ELSE 620 zztmp = grav * r1_2 621 DO jj = 2, jpjm1 622 DO ji = fs_2, fs_jpim1 ! vector opt. 623 zu_spg = zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 624 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) 625 zv_spg = zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) & 626 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj) 627 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 628 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg 629 END DO 630 END DO 631 ENDIF 632 ENDIF 633 ! !* Right-Hand-Side of the barotropic ssh equation 634 ! ! ----------------------------------------------- 635 ! ! Surface net water flux and rivers 636 IF (ln_bt_fw) THEN 637 zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) 638 ELSE 333 ! !----------------! 334 ! !== sssh_frc ==! Right-Hand-Side of the barotropic ssh equation (over the FULL domain) 335 ! !----------------! 336 ! != Net water flux forcing applied to a water column =! 337 ! ! --------------------------------------------------- ! 338 IF (ln_bt_fw) THEN ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2) 339 zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) 340 ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW) 639 341 zztmp = r1_rau0 * r1_2 640 zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) & 641 & + fwfisf(:,:) + fwfisf_b(:,:) ) 642 ENDIF 643 ! 644 IF( ln_sdw ) THEN ! Stokes drift divergence added if necessary 342 zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) + fwfisf(:,:) + fwfisf_b(:,:) ) 343 ENDIF 344 ! != Add Stokes drift divergence =! (if exist) 345 IF( ln_sdw ) THEN ! ----------------------------- ! 645 346 zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:) 646 347 ENDIF 647 348 ! 648 349 #if defined key_asminc 649 ! ! Include the IAU weighted SSH increment 350 ! != Add the IAU weighted SSH increment =! 351 ! ! ------------------------------------ ! 650 352 IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 651 353 zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:) 652 354 ENDIF 653 355 #endif 654 ! ! *Fill boundary data arrays for AGRIF356 ! != Fill boundary data arrays for AGRIF 655 357 ! ! ------------------------------------ 656 358 #if defined key_agrif … … 674 376 vb_e (:,:) = 0._wp 675 377 ENDIF 676 378 ! 379 IF( ln_linssh ) THEN ! mid-step ocean depth is fixed (hup2_e=hu_n=hu_0) 380 zhup2_e(:,:) = hu(:,:,Kmm) 381 zhvp2_e(:,:) = hv(:,:,Kmm) 382 zhtp2_e(:,:) = ht(:,:) 383 ENDIF 677 384 ! 678 385 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields … … 696 403 ENDIF 697 404 ! 698 !699 !700 405 ! Initialize sums: 701 406 puu_b (:,:,Kaa) = 0._wp ! After barotropic velocities (or transport if flux form) … … 717 422 ! 718 423 l_full_nf_update = jn == icycle ! false: disable full North fold update (performances) for jn = 1 to icycle-1 719 ! ! ------------------ 720 ! !* Update the forcing (BDY and tides) 721 ! ! ------------------ 722 ! Update only tidal forcing at open boundaries 723 IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 ) 724 IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset ) 725 ! 726 ! Set extrapolation coefficients for predictor step: 424 ! 425 ! !== Update the forcing ==! (BDY and tides) 426 ! 427 IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, kt_offset= noffset+1 ) 428 IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jn, kt_offset= noffset ) 429 ! 430 ! !== extrapolation at mid-step ==! (jn+1/2) 431 ! 432 ! !* Set extrapolation coefficients for predictor step: 727 433 IF ((jn<3).AND.ll_init) THEN ! Forward 728 434 za1 = 1._wp … … 734 440 za3 = 0.281105_wp ! za3 = bet 735 441 ENDIF 736 737 ! Extrapolate barotropic velocities at step jit+0.5: 442 ! 443 ! !* Extrapolate barotropic velocities at mid-step (jn+1/2) 444 !-- m+1/2 m m-1 m-2 --! 445 !-- u = (3/2+beta) u -(1/2+2beta) u + beta u --! 446 !-------------------------------------------------------------------------! 738 447 ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 739 448 va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) … … 742 451 ! ! ------------------ 743 452 ! Extrapolate Sea Level at step jit+0.5: 453 !-- m+1/2 m m-1 m-2 --! 454 !-- ssh = (3/2+beta) ssh -(1/2+2beta) ssh + beta ssh --! 455 !--------------------------------------------------------------------------------! 744 456 zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 745 457 746 ! set wetting & drying mask at tracer points for this barotropic sub-step 747 IF ( ln_wd_dl ) THEN 748 ! 749 IF ( ln_wd_dl_rmp ) THEN 750 DO jj = 1, jpj 751 DO ji = 1, jpi ! vector opt. 752 IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 753 ! IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 754 ztwdmask(ji,jj) = 1._wp 755 ELSE IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 756 ztwdmask(ji,jj) = (tanh(50._wp*( ( zsshp2_e(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )*r_rn_wdmin1)) ) 757 ELSE 758 ztwdmask(ji,jj) = 0._wp 759 END IF 760 END DO 761 END DO 762 ELSE 763 DO jj = 1, jpj 764 DO ji = 1, jpi ! vector opt. 765 IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 766 ztwdmask(ji,jj) = 1._wp 767 ELSE 768 ztwdmask(ji,jj) = 0._wp 769 ENDIF 770 END DO 771 END DO 772 ENDIF 773 ! 774 ENDIF 458 ! set wetting & drying mask at tracer points for this barotropic mid-step 459 IF( ln_wd_dl ) CALL wad_tmsk( zsshp2_e, ztwdmask ) 775 460 ! 776 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points 777 DO ji = 2, fs_jpim1 ! Vector opt. 778 zwx(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 779 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 780 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 781 zwy(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 782 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 783 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) 784 END DO 785 END DO 786 CALL lbc_lnk_multi( 'dynspg_ts', zwx, 'U', 1._wp, zwy, 'V', 1._wp ) 461 ! ! ocean t-depth at mid-step 462 zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:) 787 463 ! 788 zhup2_e(:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points 789 zhvp2_e(:,:) = hv_0(:,:) + zwy(:,:) 790 zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:) 791 ELSE 792 zhup2_e(:,:) = hu(:,:,Kmm) 793 zhvp2_e(:,:) = hv(:,:,Kmm) 794 zhtp2_e(:,:) = ht(:,:) 795 ENDIF 796 ! !* after ssh 797 ! ! ----------- 798 ! 799 ! Enforce volume conservation at open boundaries: 464 ! ! ocean u- and v-depth at mid-step (separate DO-loops remove the need of a lbc_lnk) 465 DO jj = 1, jpj 466 DO ji = 1, jpim1 ! not jpi-column 467 zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj) & 468 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 469 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 470 END DO 471 END DO 472 DO jj = 1, jpjm1 ! not jpj-row 473 DO ji = 1, jpi 474 zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj) & 475 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 476 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 477 END DO 478 END DO 479 ! 480 ENDIF 481 ! 482 ! !== after SSH ==! (jn+1) 483 ! 484 ! ! update (ua_e,va_e) to enforce volume conservation at open boundaries 485 ! ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d 800 486 IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e ) 801 487 ! 802 zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5 803 zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 488 ! ! resulting flux at mid-step (not over the full domain) 489 zhU(1:jpim1,1:jpj ) = e2u(1:jpim1,1:jpj ) * ua_e(1:jpim1,1:jpj ) * zhup2_e(1:jpim1,1:jpj ) ! not jpi-column 490 zhV(1:jpi ,1:jpjm1) = e1v(1:jpi ,1:jpjm1) * va_e(1:jpi ,1:jpjm1) * zhvp2_e(1:jpi ,1:jpjm1) ! not jpj-row 804 491 ! 805 492 #if defined key_agrif … … 808 495 IF((nbondi == -1).OR.(nbondi == 2)) THEN 809 496 DO jj = 1, jpj 810 z wx(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj)811 z wy(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj)497 zhU(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj) 498 zhV(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj) 812 499 END DO 813 500 ENDIF 814 501 IF((nbondi == 1).OR.(nbondi == 2)) THEN 815 502 DO jj=1,jpj 816 z wx(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj)817 z wy(nlci-nbghostcells :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells :nlci-1,jj)503 zhU(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj) 504 zhV(nlci-nbghostcells :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells :nlci-1,jj) 818 505 END DO 819 506 ENDIF 820 507 IF((nbondj == -1).OR.(nbondj == 2)) THEN 821 508 DO ji=1,jpi 822 z wy(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1)823 z wx(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1)509 zhV(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1) 510 zhU(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1) 824 511 END DO 825 512 ENDIF 826 513 IF((nbondj == 1).OR.(nbondj == 2)) THEN 827 514 DO ji=1,jpi 828 z wy(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2)829 z wx(ji,nlcj-nbghostcells :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells :nlcj-1)515 zhV(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2) 516 zhU(ji,nlcj-nbghostcells :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells :nlcj-1) 830 517 END DO 831 518 ENDIF 832 519 ENDIF 833 520 #endif 834 IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 835 836 IF ( ln_wd_dl ) THEN 837 ! 838 ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells 839 ! 840 DO jj = 1, jpjm1 841 DO ji = 1, jpim1 842 IF ( zwx(ji,jj) > 0.0 ) THEN 843 zuwdmask(ji, jj) = ztwdmask(ji ,jj) 844 ELSE 845 zuwdmask(ji, jj) = ztwdmask(ji+1,jj) 846 END IF 847 zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj) 848 un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj) 849 850 IF ( zwy(ji,jj) > 0.0 ) THEN 851 zvwdmask(ji, jj) = ztwdmask(ji, jj ) 852 ELSE 853 zvwdmask(ji, jj) = ztwdmask(ji, jj+1) 854 END IF 855 zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj) 856 vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj) 857 END DO 858 END DO 521 IF( ln_wd_il ) CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rdtbt) !!gm wad_lmt_bt use of lbc_lnk on zhU, zhV 522 523 IF( ln_wd_dl ) THEN ! un_e and vn_e are set to zero at faces where 524 ! ! the direction of the flow is from dry cells 525 CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask ) ! not jpi colomn for U, not jpj row for V 859 526 ! 860 527 ENDIF 861 862 ! Sum over sub-time-steps to compute advective velocities 863 za2 = wgtbtp2(jn) 864 un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 865 vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 866 867 ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc = True) 528 ! 529 ! 530 ! Compute Sea Level at step jit+1 531 !-- m+1 m m+1/2 --! 532 !-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --! 533 !-------------------------------------------------------------------------! 534 DO jj = 2, jpjm1 ! INNER domain 535 DO ji = 2, jpim1 536 zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj) 537 ssha_e(ji,jj) = ( sshn_e(ji,jj) - rdtbt * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj) 538 END DO 539 END DO 540 ! 541 CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp ) 542 ! 543 ! ! Sum over sub-time-steps to compute advective velocities 544 za2 = wgtbtp2(jn) ! zhU, zhV hold fluxes extrapolated at jn+0.5 545 un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:) 546 vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:) 547 ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True) 868 548 IF ( ln_wd_dl_bc ) THEN 869 zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:) 870 zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:) 871 END IF 872 873 ! Set next sea level: 874 DO jj = 2, jpjm1 875 DO ji = fs_2, fs_jpim1 ! vector opt. 876 zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) & 877 & + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1e2t(ji,jj) 878 END DO 879 END DO 880 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 881 882 CALL lbc_lnk( 'dynspg_ts', ssha_e, 'T', 1._wp ) 883 549 zuwdav2(1:jpim1,1:jpj ) = zuwdav2(1:jpim1,1:jpj ) + za2 * zuwdmask(1:jpim1,1:jpj ) ! not jpi-column 550 zvwdav2(1:jpi ,1:jpjm1) = zvwdav2(1:jpi ,1:jpjm1) + za2 * zvwdmask(1:jpi ,1:jpjm1) ! not jpj-row 551 END IF 552 ! 884 553 ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T) 885 554 IF( ln_bdy ) CALL bdy_ssh( ssha_e ) … … 890 559 ! Sea Surface Height at u-,v-points (vvl case only) 891 560 IF( .NOT.ln_linssh ) THEN 892 DO jj = 2, jpjm1 561 DO jj = 2, jpjm1 ! INNER domain, will be extended to whole domain later 893 562 DO ji = 2, jpim1 ! NO Vector Opt. 894 563 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & … … 900 569 END DO 901 570 END DO 902 CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )903 571 ENDIF 904 ! 905 ! Half-step back interpolation of SSH for surface pressure computation: 906 !---------------------------------------------------------------------- 907 IF ((jn==1).AND.ll_init) THEN 908 za0=1._wp ! Forward-backward 909 za1=0._wp 910 za2=0._wp 911 za3=0._wp 912 ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12 913 za0= 1.0833333333333_wp ! za0 = 1-gam-eps 914 za1=-0.1666666666666_wp ! za1 = gam 915 za2= 0.0833333333333_wp ! za2 = eps 916 za3= 0._wp 917 ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 918 IF (rn_bt_alpha==0._wp) THEN 919 za0=0.614_wp ! za0 = 1/2 + gam + 2*eps 920 za1=0.285_wp ! za1 = 1/2 - 2*gam - 3*eps 921 za2=0.088_wp ! za2 = gam 922 za3=0.013_wp ! za3 = eps 923 ELSE 924 zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha 925 zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha 926 za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon 927 za1 = 1._wp - za0 - zgamma - zepsilon 928 za2 = zgamma 929 za3 = zepsilon 930 ENDIF 931 ENDIF 932 ! 572 ! 573 ! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2 574 !-- m+1/2 m+1 m m-1 m-2 --! 575 !-- ssh' = za0 * ssh + za1 * ssh + za2 * ssh + za3 * ssh --! 576 !------------------------------------------------------------------------------------------! 577 CALL ts_bck_interp( jn, ll_init, za0, za1, za2, za3 ) ! coeficients of the interpolation 933 578 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 934 579 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 935 936 IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters 937 DO jj = 2, jpjm1 938 DO ji = 2, jpim1 939 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 940 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 941 & MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) & 942 & > rn_wdmin1 + rn_wdmin2 943 ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji+1,jj)) > 1.E-12 ).AND.( & 944 & MAX( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 945 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 946 947 IF(ll_tmp1) THEN 948 zcpx(ji,jj) = 1.0_wp 949 ELSE IF(ll_tmp2) THEN 950 ! no worries about zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj) = 0, it won't happen ! here 951 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 952 & / (zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj)) ) 953 ELSE 954 zcpx(ji,jj) = 0._wp 955 ENDIF 956 ! 957 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 958 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 959 & MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) & 960 & > rn_wdmin1 + rn_wdmin2 961 ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji,jj+1)) > 1.E-12 ).AND.( & 962 & MAX( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 963 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 964 965 IF(ll_tmp1) THEN 966 zcpy(ji,jj) = 1.0_wp 967 ELSEIF(ll_tmp2) THEN 968 ! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here 969 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 970 & / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj )) ) 971 ELSE 972 zcpy(ji,jj) = 0._wp 973 ENDIF 974 END DO 975 END DO 976 ENDIF 977 ! 978 ! Compute associated depths at U and V points: 979 IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form 980 ! 981 DO jj = 2, jpjm1 982 DO ji = 2, jpim1 983 zx1 = r1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) & 984 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) & 985 & + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) ) 986 zy1 = r1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) & 987 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) & 988 & + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) ) 989 zhust_e(ji,jj) = hu_0(ji,jj) + zx1 990 zhvst_e(ji,jj) = hv_0(ji,jj) + zy1 991 END DO 992 END DO 993 ! 580 ! 581 ! ! Surface pressure gradient 582 zldg = ( 1._wp - rn_scal_load ) * grav ! local factor 583 DO jj = 2, jpjm1 584 DO ji = 2, jpim1 585 zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 586 zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 587 END DO 588 END DO 589 IF( ln_wd_il ) THEN ! W/D : gravity filters applied on pressure gradient 590 CALL wad_spg( zsshp2_e, zcpx, zcpy ) ! Calculating W/D gravity filters 591 zu_spg(2:jpim1,2:jpjm1) = zu_spg(2:jpim1,2:jpjm1) * zcpx(2:jpim1,2:jpjm1) 592 zv_spg(2:jpim1,2:jpjm1) = zv_spg(2:jpim1,2:jpjm1) * zcpy(2:jpim1,2:jpjm1) 994 593 ENDIF 995 594 ! … … 997 596 ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated 998 597 ! at each time step. We however keep them constant here for optimization. 999 ! Recall that zwx and zwy arrays hold fluxes at this stage: 1000 ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5 1001 ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 1002 ! 1003 SELECT CASE( nvor_scheme ) 1004 CASE( np_ENT ) ! energy conserving scheme (t-point) 1005 DO jj = 2, jpjm1 1006 DO ji = 2, jpim1 ! vector opt. 1007 1008 z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zhup2_e(ji,jj) + 1._wp - ssumask(ji,jj) ) 1009 z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zhvp2_e(ji,jj) + 1._wp - ssvmask(ji,jj) ) 1010 1011 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu & 1012 & * ( e1e2t(ji+1,jj)*zhtp2_e(ji+1,jj)*ff_t(ji+1,jj) * ( va_e(ji+1,jj) + va_e(ji+1,jj-1) ) & 1013 & + e1e2t(ji ,jj)*zhtp2_e(ji ,jj)*ff_t(ji ,jj) * ( va_e(ji ,jj) + va_e(ji ,jj-1) ) ) 1014 ! 1015 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv & 1016 & * ( e1e2t(ji,jj+1)*zhtp2_e(ji,jj+1)*ff_t(ji,jj+1) * ( ua_e(ji,jj+1) + ua_e(ji-1,jj+1) ) & 1017 & + e1e2t(ji,jj )*zhtp2_e(ji,jj )*ff_t(ji,jj ) * ( ua_e(ji,jj ) + ua_e(ji-1,jj ) ) ) 1018 END DO 1019 END DO 1020 ! 1021 CASE( np_ENE, np_MIX ) ! energy conserving scheme (f-point) 1022 DO jj = 2, jpjm1 1023 DO ji = fs_2, fs_jpim1 ! vector opt. 1024 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 1025 zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 1026 zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 1027 zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 1028 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 1029 zv_trd(ji,jj) =-r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 1030 END DO 1031 END DO 1032 ! 1033 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 1034 DO jj = 2, jpjm1 1035 DO ji = fs_2, fs_jpim1 ! vector opt. 1036 zy1 = r1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 1037 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 1038 zx1 = - r1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 1039 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 1040 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 1041 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 1042 END DO 1043 END DO 1044 ! 1045 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 1046 DO jj = 2, jpjm1 1047 DO ji = fs_2, fs_jpim1 ! vector opt. 1048 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 1049 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 1050 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 1051 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 1052 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 1053 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 1054 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 1055 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 1056 END DO 1057 END DO 1058 ! 1059 END SELECT 598 ! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated) 599 CALL dyn_cor_2d( zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd ) 1060 600 ! 1061 601 ! Add tidal astronomical forcing if defined … … 1063 603 DO jj = 2, jpjm1 1064 604 DO ji = fs_2, fs_jpim1 ! vector opt. 1065 zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 1066 zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 1067 zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg 1068 zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg 605 zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 606 zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 1069 607 END DO 1070 608 END DO … … 1080 618 END DO 1081 619 END DO 1082 ENDIF 1083 ! 1084 ! Surface pressure trend: 1085 IF( ln_wd_il ) THEN 1086 DO jj = 2, jpjm1 1087 DO ji = 2, jpim1 1088 ! Add surface pressure gradient 1089 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 1090 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 1091 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg * zcpx(ji,jj) 1092 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg * zcpy(ji,jj) 1093 END DO 1094 END DO 1095 ELSE 1096 DO jj = 2, jpjm1 1097 DO ji = fs_2, fs_jpim1 ! vector opt. 1098 ! Add surface pressure gradient 1099 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 1100 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 1101 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg 1102 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg 1103 END DO 1104 END DO 1105 END IF 1106 620 ENDIF 1107 621 ! 1108 622 ! Set next velocities: 623 ! Compute barotropic speeds at step jit+1 (h : total height of the water colomn) 624 !-- VECTOR FORM 625 !-- m+1 m / m+1/2 \ --! 626 !-- u = u + delta_t' * \ (1-r)*g * grad_x( ssh') - f * k vect u + frc / --! 627 !-- --! 628 !-- FLUX FORM --! 629 !-- m+1 __1__ / m m / m+1/2 m+1/2 m+1/2 n \ \ --! 630 !-- u = m+1 | h * u + delta_t' * \ h * (1-r)*g * grad_x( ssh') - h * f * k vect u + h * frc / | --! 631 !-- h \ / --! 632 !------------------------------------------------------------------------------------------------------------------------! 1109 633 IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form 1110 634 DO jj = 2, jpjm1 1111 635 DO ji = fs_2, fs_jpim1 ! vector opt. 1112 636 ua_e(ji,jj) = ( un_e(ji,jj) & 1113 & + rdtbt * ( zwx(ji,jj) &637 & + rdtbt * ( zu_spg(ji,jj) & 1114 638 & + zu_trd(ji,jj) & 1115 639 & + zu_frc(ji,jj) ) & … … 1117 641 1118 642 va_e(ji,jj) = ( vn_e(ji,jj) & 1119 & + rdtbt * ( zwy(ji,jj) &643 & + rdtbt * ( zv_spg(ji,jj) & 1120 644 & + zv_trd(ji,jj) & 1121 645 & + zv_frc(ji,jj) ) & 1122 646 & ) * ssvmask(ji,jj) 1123 1124 647 END DO 1125 648 END DO … … 1127 650 ELSE !* Flux form 1128 651 DO jj = 2, jpjm1 1129 DO ji = fs_2, fs_jpim1 ! vector opt. 1130 1131 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 1132 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 1133 1134 zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj)) 1135 zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj)) 1136 1137 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) & 1138 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 1139 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 1140 & + hu(ji,jj,Kmm) * zu_frc(ji,jj) ) & 1141 & ) * zhura 1142 1143 va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) & 1144 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 1145 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & 1146 & + hv(ji,jj,Kmm) * zv_frc(ji,jj) ) & 1147 & ) * zhvra 652 DO ji = 2, jpim1 653 ! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2 654 ! ! backward interpolated depth used in spg terms at jn+1/2 655 zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 656 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 657 zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 658 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 659 ! ! inverse depth at jn+1 660 z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) ) 661 z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) ) 662 ! 663 ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) & 664 & + rdtbt * ( zhu_bck * zu_spg (ji,jj) & ! 665 & + zhup2_e(ji,jj) * zu_trd (ji,jj) & ! 666 & + hu(ji,jj,Kmm) * zu_frc (ji,jj) ) ) * z1_hu 667 ! 668 va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) & 669 & + rdtbt * ( zhv_bck * zv_spg (ji,jj) & ! 670 & + zhvp2_e(ji,jj) * zv_trd (ji,jj) & ! 671 & + hv(ji,jj,Kmm) * zv_frc (ji,jj) ) ) * z1_hv 1148 672 END DO 1149 673 END DO … … 1158 682 END DO 1159 683 ENDIF 1160 1161 1162 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 1163 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 1164 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 1165 hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 1166 hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) 1167 ! 1168 ENDIF 1169 ! !* domain lateral boundary 1170 CALL lbc_lnk_multi( 'dynspg_ts', ua_e, 'U', -1._wp, va_e , 'V', -1._wp ) 684 685 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 686 hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1) 687 hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) ) 688 hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1) 689 hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) ) 690 CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp & 691 & , hu_e , 'U', -1._wp, hv_e , 'V', -1._wp & 692 & , hur_e, 'U', -1._wp, hvr_e, 'V', -1._wp ) 693 ELSE 694 CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp ) 695 ENDIF 696 ! 1171 697 ! 1172 698 ! ! open boundaries … … 1216 742 ! Set advection velocity correction: 1217 743 IF (ln_bt_fw) THEN 1218 zwx(:,:) = un_adv(:,:)1219 zwy(:,:) = vn_adv(:,:)1220 744 IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN 1221 un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) ) 1222 vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) ) 1223 ! 1224 ! Update corrective fluxes for next time step: 1225 un_bf(:,:) = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:)) 1226 vn_bf(:,:) = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:)) 745 DO jj = 1, jpj 746 DO ji = 1, jpi 747 zun_save = un_adv(ji,jj) 748 zvn_save = vn_adv(ji,jj) 749 ! ! apply the previously computed correction 750 un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - atfp * un_bf(ji,jj) ) 751 vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - atfp * vn_bf(ji,jj) ) 752 ! ! Update corrective fluxes for next time step 753 un_bf(ji,jj) = atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) ) 754 vn_bf(ji,jj) = atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) ) 755 ! ! Save integrated transport for next computation 756 ub2_b(ji,jj) = zun_save 757 vb2_b(ji,jj) = zvn_save 758 END DO 759 END DO 1227 760 ELSE 1228 un_bf(:,:) = 0._wp 1229 vn_bf(:,:) = 0._wp 1230 END IF 1231 ! Save integrated transport for next computation 1232 ub2_b(:,:) = zwx(:,:) 1233 vb2_b(:,:) = zwy(:,:) 761 un_bf(:,:) = 0._wp ! corrective fluxes for next time step set to zero 762 vn_bf(:,:) = 0._wp 763 ub2_b(:,:) = un_adv(:,:) ! Save integrated transport for next computation 764 vb2_b(:,:) = vn_adv(:,:) 765 END IF 1234 766 ENDIF 1235 767 … … 1307 839 ! 1308 840 IF( ln_diatmb ) THEN 1309 CALL iom_put( "baro_u" , uu_b(:,:,Kmm)*ssumask(:,:)+zmdi*(1.-ssumask(:,:) ) ) ! Barotropic U Velocity1310 CALL iom_put( "baro_v" , vv_b(:,:,Kmm)*ssvmask(:,:)+zmdi*(1.-ssvmask(:,:) ) ) ! Barotropic V Velocity841 CALL iom_put( "baro_u" , puu_b(:,:,Kmm)*ssumask(:,:)+zmdi*(1.-ssumask(:,:) ) ) ! Barotropic U Velocity 842 CALL iom_put( "baro_v" , pvv_b(:,:,Kmm)*ssvmask(:,:)+zmdi*(1.-ssvmask(:,:) ) ) ! Barotropic V Velocity 1311 843 ENDIF 1312 844 ! … … 1582 1114 END SUBROUTINE dyn_spg_ts_init 1583 1115 1116 1117 SUBROUTINE dyn_cor_2D_init( Kmm ) 1118 !!--------------------------------------------------------------------- 1119 !! *** ROUTINE dyn_cor_2D_init *** 1120 !! 1121 !! ** Purpose : Set time splitting options 1122 !! Set arrays to remove/compute coriolis trend. 1123 !! Do it once during initialization if volume is fixed, else at each long time step. 1124 !! Note that these arrays are also used during barotropic loop. These are however frozen 1125 !! although they should be updated in the variable volume case. Not a big approximation. 1126 !! To remove this approximation, copy lines below inside barotropic loop 1127 !! and update depths at T-F points (ht and zhf resp.) at each barotropic time step 1128 !! 1129 !! Compute zwz = f / ( height of the water colomn ) 1130 !!---------------------------------------------------------------------- 1131 INTEGER, INTENT(in) :: Kmm ! Time index 1132 INTEGER :: ji ,jj, jk ! dummy loop indices 1133 REAL(wp) :: z1_ht 1134 REAL(wp), DIMENSION(jpi,jpj) :: zhf 1135 !!---------------------------------------------------------------------- 1136 ! 1137 SELECT CASE( nvor_scheme ) 1138 CASE( np_EEN ) != EEN scheme using e3f (energy & enstrophy scheme) 1139 SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point 1140 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 1141 DO jj = 1, jpjm1 1142 DO ji = 1, jpim1 1143 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 1144 & ht(ji ,jj ) + ht(ji+1,jj ) ) * 0.25_wp 1145 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1146 END DO 1147 END DO 1148 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 1149 DO jj = 1, jpjm1 1150 DO ji = 1, jpim1 1151 zwz(ji,jj) = ( ht (ji ,jj+1) + ht (ji+1,jj+1) & 1152 & + ht (ji ,jj ) + ht (ji+1,jj ) ) & 1153 & / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) & 1154 & + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) ) 1155 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1156 END DO 1157 END DO 1158 END SELECT 1159 CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp ) 1160 ! 1161 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1162 DO jj = 2, jpj 1163 DO ji = 2, jpi 1164 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 1165 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) 1166 ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) 1167 ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) 1168 END DO 1169 END DO 1170 ! 1171 CASE( np_EET ) != EEN scheme using e3t (energy conserving scheme) 1172 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1173 DO jj = 2, jpj 1174 DO ji = 2, jpi 1175 z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) ) 1176 ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht 1177 ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht 1178 ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht 1179 ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht 1180 END DO 1181 END DO 1182 ! 1183 CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT ! 1184 ! 1185 zwz(:,:) = 0._wp 1186 zhf(:,:) = 0._wp 1187 1188 !!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed 1189 !!gm A priori a better value should be something like : 1190 !!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1) 1191 !!gm divided by the sum of the corresponding mask 1192 !!gm 1193 !! 1194 IF( .NOT.ln_sco ) THEN 1195 1196 !!gm agree the JC comment : this should be done in a much clear way 1197 1198 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 1199 ! Set it to zero for the time being 1200 ! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level 1201 ! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth 1202 ! ENDIF 1203 ! zhf(:,:) = gdepw_0(:,:,jk+1) 1204 ! 1205 ELSE 1206 ! 1207 !zhf(:,:) = hbatf(:,:) 1208 DO jj = 1, jpjm1 1209 DO ji = 1, jpim1 1210 zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) & 1211 & + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) & 1212 & / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) & 1213 & + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp ) 1214 END DO 1215 END DO 1216 ENDIF 1217 ! 1218 DO jj = 1, jpjm1 1219 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 1220 END DO 1221 ! 1222 DO jk = 1, jpkm1 1223 DO jj = 1, jpjm1 1224 zhf(:,jj) = zhf(:,jj) + e3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 1225 END DO 1226 END DO 1227 CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp ) 1228 ! JC: TBC. hf should be greater than 0 1229 DO jj = 1, jpj 1230 DO ji = 1, jpi 1231 IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) 1232 END DO 1233 END DO 1234 zwz(:,:) = ff_f(:,:) * zwz(:,:) 1235 END SELECT 1236 1237 END SUBROUTINE dyn_cor_2d_init 1238 1239 1240 1241 SUBROUTINE dyn_cor_2d( phu, phv, punb, pvnb, zhU, zhV, zu_trd, zv_trd ) 1242 !!--------------------------------------------------------------------- 1243 !! *** ROUTINE dyn_cor_2d *** 1244 !! 1245 !! ** Purpose : Compute u and v coriolis trends 1246 !!---------------------------------------------------------------------- 1247 INTEGER :: ji ,jj ! dummy loop indices 1248 REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! - - 1249 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: phu, phv, punb, pvnb, zhU, zhV 1250 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd 1251 !!---------------------------------------------------------------------- 1252 SELECT CASE( nvor_scheme ) 1253 CASE( np_ENT ) ! enstrophy conserving scheme (f-point) 1254 DO jj = 2, jpjm1 1255 DO ji = 2, jpim1 1256 z1_hu = ssumask(ji,jj) / ( phu(ji,jj) + 1._wp - ssumask(ji,jj) ) 1257 z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) ) 1258 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu & 1259 & * ( e1e2t(ji+1,jj)*ht(ji+1,jj)*ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) ) & 1260 & + e1e2t(ji ,jj)*ht(ji ,jj)*ff_t(ji ,jj) * ( pvnb(ji ,jj) + pvnb(ji ,jj-1) ) ) 1261 ! 1262 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv & 1263 & * ( e1e2t(ji,jj+1)*ht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) ) & 1264 & + e1e2t(ji,jj )*ht(ji,jj )*ff_t(ji,jj ) * ( punb(ji,jj ) + punb(ji-1,jj ) ) ) 1265 END DO 1266 END DO 1267 ! 1268 CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX 1269 DO jj = 2, jpjm1 1270 DO ji = fs_2, fs_jpim1 ! vector opt. 1271 zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj) 1272 zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) 1273 zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj) 1274 zx2 = ( zhU(ji ,jj) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj) 1275 ! energy conserving formulation for planetary vorticity term 1276 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 1277 zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 1278 END DO 1279 END DO 1280 ! 1281 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 1282 DO jj = 2, jpjm1 1283 DO ji = fs_2, fs_jpim1 ! vector opt. 1284 zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) & 1285 & + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) 1286 zx1 = - r1_8 * ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) & 1287 & + zhU(ji ,jj ) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj) 1288 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 1289 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 1290 END DO 1291 END DO 1292 ! 1293 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 1294 DO jj = 2, jpjm1 1295 DO ji = fs_2, fs_jpim1 ! vector opt. 1296 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) & 1297 & + ftnw(ji+1,jj) * zhV(ji+1,jj ) & 1298 & + ftse(ji,jj ) * zhV(ji ,jj-1) & 1299 & + ftsw(ji+1,jj) * zhV(ji+1,jj-1) ) 1300 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) & 1301 & + ftse(ji,jj+1) * zhU(ji ,jj+1) & 1302 & + ftnw(ji,jj ) * zhU(ji-1,jj ) & 1303 & + ftne(ji,jj ) * zhU(ji ,jj ) ) 1304 END DO 1305 END DO 1306 ! 1307 END SELECT 1308 ! 1309 END SUBROUTINE dyn_cor_2D 1310 1311 1312 SUBROUTINE wad_tmsk( pssh, ptmsk ) 1313 !!---------------------------------------------------------------------- 1314 !! *** ROUTINE wad_lmt *** 1315 !! 1316 !! ** Purpose : set wetting & drying mask at tracer points 1317 !! for the current barotropic sub-step 1318 !! 1319 !! ** Method : ??? 1320 !! 1321 !! ** Action : ptmsk : wetting & drying t-mask 1322 !!---------------------------------------------------------------------- 1323 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pssh ! 1324 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: ptmsk ! 1325 ! 1326 INTEGER :: ji, jj ! dummy loop indices 1327 !!---------------------------------------------------------------------- 1328 ! 1329 IF( ln_wd_dl_rmp ) THEN 1330 DO jj = 1, jpj 1331 DO ji = 1, jpi 1332 IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 1333 ! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 1334 ptmsk(ji,jj) = 1._wp 1335 ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 1336 ptmsk(ji,jj) = TANH( 50._wp*( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )*r_rn_wdmin1) ) 1337 ELSE 1338 ptmsk(ji,jj) = 0._wp 1339 ENDIF 1340 END DO 1341 END DO 1342 ELSE 1343 DO jj = 1, jpj 1344 DO ji = 1, jpi 1345 IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp 1346 ELSE ; ptmsk(ji,jj) = 0._wp 1347 ENDIF 1348 END DO 1349 END DO 1350 ENDIF 1351 ! 1352 END SUBROUTINE wad_tmsk 1353 1354 1355 SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk ) 1356 !!---------------------------------------------------------------------- 1357 !! *** ROUTINE wad_lmt *** 1358 !! 1359 !! ** Purpose : set wetting & drying mask at tracer points 1360 !! for the current barotropic sub-step 1361 !! 1362 !! ** Method : ??? 1363 !! 1364 !! ** Action : ptmsk : wetting & drying t-mask 1365 !!---------------------------------------------------------------------- 1366 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pTmsk ! W & D t-mask 1367 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phU, phV, pu, pv ! ocean velocities and transports 1368 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pUmsk, pVmsk ! W & D u- and v-mask 1369 ! 1370 INTEGER :: ji, jj ! dummy loop indices 1371 !!---------------------------------------------------------------------- 1372 ! 1373 DO jj = 1, jpj 1374 DO ji = 1, jpim1 ! not jpi-column 1375 IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj) 1376 ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj) 1377 ENDIF 1378 phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj) 1379 pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj) 1380 END DO 1381 END DO 1382 ! 1383 DO jj = 1, jpjm1 ! not jpj-row 1384 DO ji = 1, jpi 1385 IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj ) 1386 ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1) 1387 ENDIF 1388 phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 1389 pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj) 1390 END DO 1391 END DO 1392 ! 1393 END SUBROUTINE wad_Umsk 1394 1395 1396 SUBROUTINE wad_spg( pshn, zcpx, zcpy ) 1397 !!--------------------------------------------------------------------- 1398 !! *** ROUTINE wad_sp *** 1399 !! 1400 !! ** Purpose : 1401 !!---------------------------------------------------------------------- 1402 INTEGER :: ji ,jj ! dummy loop indices 1403 LOGICAL :: ll_tmp1, ll_tmp2 1404 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pshn 1405 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy 1406 !!---------------------------------------------------------------------- 1407 DO jj = 2, jpjm1 1408 DO ji = 2, jpim1 1409 ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji+1,jj) ) > & 1410 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 1411 & MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji+1,jj) + ht_0(ji+1,jj) ) & 1412 & > rn_wdmin1 + rn_wdmin2 1413 ll_tmp2 = ( ABS( pshn(ji+1,jj) - pshn(ji ,jj)) > 1.E-12 ).AND.( & 1414 & MAX( pshn(ji,jj) , pshn(ji+1,jj) ) > & 1415 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 1416 IF(ll_tmp1) THEN 1417 zcpx(ji,jj) = 1.0_wp 1418 ELSEIF(ll_tmp2) THEN 1419 ! no worries about pshn(ji+1,jj) - pshn(ji ,jj) = 0, it won't happen ! here 1420 zcpx(ji,jj) = ABS( (pshn(ji+1,jj) + ht_0(ji+1,jj) - pshn(ji,jj) - ht_0(ji,jj)) & 1421 & / (pshn(ji+1,jj) - pshn(ji ,jj)) ) 1422 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 1423 ELSE 1424 zcpx(ji,jj) = 0._wp 1425 ENDIF 1426 ! 1427 ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji,jj+1) ) > & 1428 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 1429 & MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji,jj+1) + ht_0(ji,jj+1) ) & 1430 & > rn_wdmin1 + rn_wdmin2 1431 ll_tmp2 = ( ABS( pshn(ji,jj) - pshn(ji,jj+1)) > 1.E-12 ).AND.( & 1432 & MAX( pshn(ji,jj) , pshn(ji,jj+1) ) > & 1433 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1434 1435 IF(ll_tmp1) THEN 1436 zcpy(ji,jj) = 1.0_wp 1437 ELSE IF(ll_tmp2) THEN 1438 ! no worries about pshn(ji,jj+1) - pshn(ji,jj ) = 0, it won't happen ! here 1439 zcpy(ji,jj) = ABS( (pshn(ji,jj+1) + ht_0(ji,jj+1) - pshn(ji,jj) - ht_0(ji,jj)) & 1440 & / (pshn(ji,jj+1) - pshn(ji,jj )) ) 1441 zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) ) 1442 ELSE 1443 zcpy(ji,jj) = 0._wp 1444 ENDIF 1445 END DO 1446 END DO 1447 1448 END SUBROUTINE wad_spg 1449 1450 1451 1452 SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v ) 1453 !!---------------------------------------------------------------------- 1454 !! *** ROUTINE dyn_drg_init *** 1455 !! 1456 !! ** Purpose : - add the baroclinic top/bottom drag contribution to 1457 !! the baroclinic part of the barotropic RHS 1458 !! - compute the barotropic drag coefficients 1459 !! 1460 !! ** Method : computation done over the INNER domain only 1461 !!---------------------------------------------------------------------- 1462 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 1463 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in ) :: puu, pvv ! ocean velocities and RHS of momentum equation 1464 REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(in ) :: puu_b, pvv_b ! barotropic velocities at main time levels 1465 REAL(wp), DIMENSION(jpi,jpj) , INTENT(inout) :: pu_RHSi, pv_RHSi ! baroclinic part of the barotropic RHS 1466 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients 1467 ! 1468 INTEGER :: ji, jj ! dummy loop indices 1469 INTEGER :: ikbu, ikbv, iktu, iktv 1470 REAL(wp) :: zztmp 1471 REAL(wp), DIMENSION(jpi,jpj) :: zu_i, zv_i 1472 !!---------------------------------------------------------------------- 1473 ! 1474 ! !== Set the barotropic drag coef. ==! 1475 ! 1476 IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities) 1477 1478 DO jj = 2, jpjm1 1479 DO ji = 2, jpim1 ! INNER domain 1480 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 1481 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) 1482 END DO 1483 END DO 1484 ELSE ! bottom friction only 1485 DO jj = 2, jpjm1 1486 DO ji = 2, jpim1 ! INNER domain 1487 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 1488 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) 1489 END DO 1490 END DO 1491 ENDIF 1492 ! 1493 ! !== BOTTOM stress contribution from baroclinic velocities ==! 1494 ! 1495 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities 1496 1497 DO jj = 2, jpjm1 1498 DO ji = 2, jpim1 ! INNER domain 1499 ikbu = mbku(ji,jj) 1500 ikbv = mbkv(ji,jj) 1501 zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm) 1502 zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm) 1503 END DO 1504 END DO 1505 ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities 1506 1507 DO jj = 2, jpjm1 1508 DO ji = 2, jpim1 ! INNER domain 1509 ikbu = mbku(ji,jj) 1510 ikbv = mbkv(ji,jj) 1511 zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb) 1512 zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb) 1513 END DO 1514 END DO 1515 ENDIF 1516 ! 1517 IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please ! 1518 zztmp = -1._wp / rdtbt 1519 DO jj = 2, jpjm1 1520 DO ji = 2, jpim1 ! INNER domain 1521 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( & 1522 & r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) 1523 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( & 1524 & r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp ) 1525 END DO 1526 END DO 1527 ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation) 1528 1529 DO jj = 2, jpjm1 1530 DO ji = 2, jpim1 ! INNER domain 1531 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj) 1532 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj) 1533 END DO 1534 END DO 1535 END IF 1536 ! 1537 ! !== TOP stress contribution from baroclinic velocities ==! (no W/D case) 1538 ! 1539 IF( ln_isfcav ) THEN 1540 ! 1541 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity 1542 1543 DO jj = 2, jpjm1 1544 DO ji = 2, jpim1 ! INNER domain 1545 iktu = miku(ji,jj) 1546 iktv = mikv(ji,jj) 1547 zu_i(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm) 1548 zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm) 1549 END DO 1550 END DO 1551 ELSE ! CENTRED integration: use BEFORE top baroclinic velocity 1552 1553 DO jj = 2, jpjm1 1554 DO ji = 2, jpim1 ! INNER domain 1555 iktu = miku(ji,jj) 1556 iktv = mikv(ji,jj) 1557 zu_i(ji,jj) = puu(ji,jj,iktu,Kbb) - puu_b(ji,jj,Kbb) 1558 zv_i(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb) 1559 END DO 1560 END DO 1561 ENDIF 1562 ! 1563 ! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation) 1564 1565 DO jj = 2, jpjm1 1566 DO ji = 2, jpim1 ! INNER domain 1567 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj) 1568 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj) 1569 END DO 1570 END DO 1571 ! 1572 ENDIF 1573 ! 1574 END SUBROUTINE dyn_drg_init 1575 1576 SUBROUTINE ts_bck_interp( jn, ll_init, & ! <== in 1577 & za0, za1, za2, za3 ) ! ==> out 1578 !!---------------------------------------------------------------------- 1579 INTEGER ,INTENT(in ) :: jn ! index of sub time step 1580 LOGICAL ,INTENT(in ) :: ll_init ! 1581 REAL(wp),INTENT( out) :: za0, za1, za2, za3 ! Half-step back interpolation coefficient 1582 ! 1583 REAL(wp) :: zepsilon, zgamma ! - - 1584 !!---------------------------------------------------------------------- 1585 ! ! set Half-step back interpolation coefficient 1586 IF ( jn==1 .AND. ll_init ) THEN !* Forward-backward 1587 za0 = 1._wp 1588 za1 = 0._wp 1589 za2 = 0._wp 1590 za3 = 0._wp 1591 ELSEIF( jn==2 .AND. ll_init ) THEN !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12 1592 za0 = 1.0833333333333_wp ! za0 = 1-gam-eps 1593 za1 =-0.1666666666666_wp ! za1 = gam 1594 za2 = 0.0833333333333_wp ! za2 = eps 1595 za3 = 0._wp 1596 ELSE !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 1597 IF( rn_bt_alpha == 0._wp ) THEN ! Time diffusion 1598 za0 = 0.614_wp ! za0 = 1/2 + gam + 2*eps 1599 za1 = 0.285_wp ! za1 = 1/2 - 2*gam - 3*eps 1600 za2 = 0.088_wp ! za2 = gam 1601 za3 = 0.013_wp ! za3 = eps 1602 ELSE ! no time diffusion 1603 zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha 1604 zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha 1605 za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon 1606 za1 = 1._wp - za0 - zgamma - zepsilon 1607 za2 = zgamma 1608 za3 = zepsilon 1609 ENDIF 1610 ENDIF 1611 END SUBROUTINE ts_bck_interp 1612 1613 1584 1614 !!====================================================================== 1585 1615 END MODULE dynspg_ts -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynvor.F90
r10946 r11822 858 858 REWIND( numnam_ref ) ! Namelist namdyn_vor in reference namelist : Vorticity scheme options 859 859 READ ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901) 860 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in reference namelist' , lwp)860 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in reference namelist' ) 861 861 REWIND( numnam_cfg ) ! Namelist namdyn_vor in configuration namelist : Vorticity scheme options 862 862 READ ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 ) 863 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist' , lwp)863 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist' ) 864 864 IF(lwm) WRITE ( numond, namdyn_vor ) 865 865 ! -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzdf.F90
r10946 r11822 172 172 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 173 173 & / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 174 zWui = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) )175 zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )174 zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) / ze3ua 175 zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 176 176 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 177 177 zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) … … 187 187 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 188 188 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 189 zWui = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) )190 zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )189 zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) / ze3ua 190 zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 191 191 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 192 192 zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) … … 201 201 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa) 202 202 zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2) 203 zWus = 0.5_wp * ( wi(ji ,jj,2) + wi(ji+1,jj,2) )203 zWus = ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) / ze3ua 204 204 zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp ) 205 205 zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) ) … … 338 338 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 339 339 & / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 340 zWvi = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * wvmask(ji,jj,jk )341 zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1)340 zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) / ze3va 341 zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 342 342 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp ) 343 343 zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp ) … … 353 353 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 354 354 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 355 zWvi = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * wvmask(ji,jj,jk )356 zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1)355 zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) / ze3va 356 zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 357 357 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp ) 358 358 zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp ) … … 367 367 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa) 368 368 zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2) 369 zWvs = 0.5_wp * ( wi(ji,jj ,2) + wi(ji,jj+1,2) )369 zWvs = ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) / ze3va 370 370 zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp ) 371 371 zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) ) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/sshwzv.F90
r11480 r11822 9 9 !! - ! 2010-09 (D.Storkey and E.O'Dea) bug fixes for BDY module 10 10 !! 3.3 ! 2011-10 (M. Leclair) split former ssh_wzv routine and remove all vvl related work 11 !! 4.0 ! 2018-12 (A. Coward) add mixed implicit/explicit advection 11 12 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) Rename ssh_nxt -> ssh_atf. Now only does time filtering. 12 13 !!---------------------------------------------------------------------- … … 278 279 !! : wi : now vertical velocity (for implicit treatment) 279 280 !! 280 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 281 !! Reference : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent 282 !! implicit scheme for vertical advection in oceanic modeling. 283 !! Ocean Modelling, 91, 38-69. 281 284 !!---------------------------------------------------------------------- 282 285 INTEGER, INTENT(in) :: kt ! time step … … 284 287 ! 285 288 INTEGER :: ji, jj, jk ! dummy loop indices 286 REAL(wp) :: zCu, zcff, z1_e3 w! local scalars289 REAL(wp) :: zCu, zcff, z1_e3t ! local scalars 287 290 REAL(wp) , PARAMETER :: Cu_min = 0.15_wp ! local parameters 288 REAL(wp) , PARAMETER :: Cu_max = 0. 27! local parameters291 REAL(wp) , PARAMETER :: Cu_max = 0.30_wp ! local parameters 289 292 REAL(wp) , PARAMETER :: Cu_cut = 2._wp*Cu_max - Cu_min ! local parameters 290 293 REAL(wp) , PARAMETER :: Fcu = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters … … 297 300 IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity ' 298 301 IF(lwp) WRITE(numout,*) '~~~~~ ' 299 ! 300 Cu_adv(:,:,jpk) = 0._wp ! bottom value : Cu_adv=0 (set once for all) 301 ENDIF 302 ! 303 DO jk = 1, jpkm1 ! calculate Courant numbers 304 DO jj = 2, jpjm1 305 DO ji = 2, fs_jpim1 ! vector opt. 306 z1_e3w = 1._wp / e3w(ji,jj,jk,Kmm) 307 Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) & 308 & + ( MAX( e2u(ji ,jj)*e3uw(ji ,jj,jk,Kmm)*uu(ji ,jj,jk,Kmm), 0._wp ) - & 309 & MIN( e2u(ji-1,jj)*e3uw(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) ) & 310 & * r1_e1e2t(ji,jj) & 311 & + ( MAX( e1v(ji,jj )*e3vw(ji,jj ,jk,Kmm)*vv(ji,jj ,jk,Kmm), 0._wp ) - & 312 & MIN( e1v(ji,jj-1)*e3vw(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm), 0._wp ) ) & 313 & * r1_e1e2t(ji,jj) & 314 & ) * z1_e3w 302 wi(:,:,:) = 0._wp 303 ENDIF 304 ! 305 ! Calculate Courant numbers 306 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 307 DO jk = 1, jpkm1 308 DO jj = 2, jpjm1 309 DO ji = 2, fs_jpim1 ! vector opt. 310 z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 311 ! 2*rdt and not r2dt (for restartability) 312 Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) & 313 & + ( MAX( e2u(ji ,jj)*e3u(ji ,jj,jk,Kmm)*uu(ji ,jj,jk,Kmm) + un_td(ji ,jj,jk), 0._wp ) - & 314 & MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) ) & 315 & * r1_e1e2t(ji,jj) & 316 & + ( MAX( e1v(ji,jj )*e3v(ji,jj ,jk,Kmm)*vv(ji,jj ,jk,Kmm) + vn_td(ji,jj ,jk), 0._wp ) - & 317 & MIN( e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm) + vn_td(ji,jj-1,jk), 0._wp ) ) & 318 & * r1_e1e2t(ji,jj) & 319 & ) * z1_e3t 320 END DO 315 321 END DO 316 322 END DO 317 END DO 323 ELSE 324 DO jk = 1, jpkm1 325 DO jj = 2, jpjm1 326 DO ji = 2, fs_jpim1 ! vector opt. 327 z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 328 ! 2*rdt and not r2dt (for restartability) 329 Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) & 330 & + ( MAX( e2u(ji ,jj)*e3u(ji ,jj,jk,Kmm)*uu(ji ,jj,jk,Kmm), 0._wp ) - & 331 & MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) ) & 332 & * r1_e1e2t(ji,jj) & 333 & + ( MAX( e1v(ji,jj )*e3v(ji,jj ,jk,Kmm)*vv(ji,jj ,jk,Kmm), 0._wp ) - & 334 & MIN( e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm), 0._wp ) ) & 335 & * r1_e1e2t(ji,jj) & 336 & ) * z1_e3t 337 END DO 338 END DO 339 END DO 340 ENDIF 341 CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1. ) 318 342 ! 319 343 CALL iom_put("Courant",Cu_adv) 320 344 ! 321 wi(:,:,:) = 0._wp ! Includes top and bottom values set to zero322 345 IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN ! Quick check if any breaches anywhere 323 DO jk = 1, jpkm1! or scan Courant criterion and partition324 DO jj = 2, jpjm1! w where necessary325 DO ji = 2, fs_jpim1 ! vector opt.346 DO jk = jpkm1, 2, -1 ! or scan Courant criterion and partition 347 DO jj = 1, jpj ! w where necessary 348 DO ji = 1, jpi 326 349 ! 327 zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) ) 350 zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) ) 351 ! alt: 352 ! IF ( wn(ji,jj,jk) > 0._wp ) THEN 353 ! zCu = Cu_adv(ji,jj,jk) 354 ! ELSE 355 ! zCu = Cu_adv(ji,jj,jk-1) 356 ! ENDIF 328 357 ! 329 IF( zCu < Cu_min ) THEN!<-- Fully explicit358 IF( zCu <= Cu_min ) THEN !<-- Fully explicit 330 359 zcff = 0._wp 331 360 ELSEIF( zCu < Cu_cut ) THEN !<-- Mixed explicit … … 340 369 ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk) 341 370 ! 342 Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient 371 Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient below and in stp_ctl 343 372 END DO 344 373 END DO 345 374 END DO 375 Cu_adv(:,:,1) = 0._wp 346 376 ELSE 347 377 ! Fully explicit everywhere 348 Cu_adv = 0.0_wp ! Reuse array to output coefficient 378 Cu_adv(:,:,:) = 0._wp ! Reuse array to output coefficient below and in stp_ctl 379 wi (:,:,:) = 0._wp 349 380 ENDIF 350 381 CALL iom_put("wimp",wi) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/wet_dry.F90
r11027 r11822 81 81 REWIND( numnam_ref ) ! Namelist namwad in reference namelist : Parameters for Wetting/Drying 82 82 READ ( numnam_ref, namwad, IOSTAT = ios, ERR = 905) 83 905 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in reference namelist' , .TRUE.)83 905 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in reference namelist' ) 84 84 REWIND( numnam_cfg ) ! Namelist namwad in configuration namelist : Parameters for Wetting/Drying 85 85 READ ( numnam_cfg, namwad, IOSTAT = ios, ERR = 906) 86 906 IF( ios > 0 ) CALL ctl_nam ( ios , 'namwad in configuration namelist' , .TRUE.)86 906 IF( ios > 0 ) CALL ctl_nam ( ios , 'namwad in configuration namelist' ) 87 87 IF(lwm) WRITE ( numond, namwad ) 88 88 !
Note: See TracChangeset
for help on using the changeset viewer.
