Changeset 11380 for NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/DYN/dynspg_ts.F90
- Timestamp:
- 2019-07-31T15:56:02+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/DYN/dynspg_ts.F90
r11372 r11380 64 64 USE diatmb ! Top,middle,bottom output 65 65 66 USE iom ! to remove67 66 68 67 IMPLICIT NONE … … 78 77 REAL(wp),SAVE :: rdtbt ! Barotropic time step 79 78 ! 80 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 ! 1st & 2nd weights used in time filtering of barotropic fields 81 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff_f/h at F points 82 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter 83 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme) 79 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 ! 1st & 2nd weights used in time filtering of barotropic fields 80 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff_f/h at F points 81 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter 82 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme) 83 84 !! Arrays at barotropic time step: ! befbefore! before ! now ! after ! 85 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ubb_e , ub_e , un_e , ua_e !: u-external velocity 86 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vbb_e , vb_e , vn_e , va_e !: v-external velocity 87 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshbb_e, sshb_e, sshn_e, ssha_e !: external ssh 88 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_e !: external u-depth 89 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hv_e !: external v-depth 90 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hur_e !: inverse of u-depth 91 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hvr_e !: inverse of v-depth 92 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub2_b , vb2_b !: Half step fluxes (ln_bt_fw=T) 93 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_bf , vn_bf !: Asselin filtered half step fluxes (ln_bt_fw=T) 94 95 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_0_xtd , hu_0_xtd , hv_0_xtd 96 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_n_xtd , hv_n_xtd 97 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_e1e2t_xtd, r1_e1e2u_xtd, r1_e1e2v_xtd 98 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: e1e2t_xtd 99 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: e2u_xtd , e1v_xtd 100 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_e1u_xtd, r1_e2v_xtd 101 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssmask_xtd, ssumask_xtd, ssvmask_xtd 102 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ff_t_xtd ! used in ENT scheme 103 104 #if defined key_agrif 105 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub2_i_b, vb2_i_b !: Half step time integrated fluxes 106 #endif 84 107 85 108 REAL(wp) :: r1_12 = 1._wp / 12._wp ! local ratios … … 101 124 !! *** routine dyn_spg_ts_alloc *** 102 125 !!---------------------------------------------------------------------- 103 INTEGER :: ierr(3) 126 INTEGER :: ierr(6) 127 INTEGER :: idbi, idei, idbj, idej ! lower/upper bounds of extended arrays 104 128 !!---------------------------------------------------------------------- 105 129 ierr(:) = 0 106 ! 107 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) ) 108 IF( ln_dynvor_een .OR. ln_dynvor_eeT ) & 109 & ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , ftsw(jpi,jpj) , ftse(jpi,jpj), STAT=ierr(2) ) 130 idbi = 1 - nn_hlts ; idbj = 1 - nn_hlts 131 idei = jpi + nn_hlts ; idej = jpj + nn_hlts 132 ! 133 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(idbi:idei,idbj:idej), STAT=ierr(1) ) 134 IF( ln_dynvor_een .OR. ln_dynvor_eeT ) THEN 135 ALLOCATE( ftnw(idbi:idei,idbj:idej) , ftne(idbi:idei,idbj:idej) , ftsw(idbi:idei,idbj:idej) , ftse(idbi:idei,idbj:idej) & 136 & , STAT=ierr(2) ) 137 ELSEIF( ln_dynvor_enT ) THEN 138 ALLOCATE( ff_t_xtd(idbi:idei,idbj:idej), STAT=ierr(2) ) 139 END IF 110 140 ! 111 141 ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(3) ) 142 ALLOCATE( ssha_e(idbi:idei,idbj:idej), sshn_e(idbi:idei,idbj:idej), sshb_e(idbi:idei,idbj:idej), sshbb_e(idbi:idei,idbj:idej) & 143 & , ua_e(idbi:idei,idbj:idej), un_e(idbi:idei,idbj:idej), ub_e(idbi:idei,idbj:idej), ubb_e(idbi:idei,idbj:idej) & 144 & , va_e(idbi:idei,idbj:idej), vn_e(idbi:idei,idbj:idej), vb_e(idbi:idei,idbj:idej), vbb_e(idbi:idei,idbj:idej) & 145 & , hu_e(idbi:idei,idbj:idej), hur_e(idbi:idei,idbj:idej), hv_e(idbi:idei,idbj:idej), hvr_e(idbi:idei,idbj:idej) & 146 & , STAT=ierr(4) ) 147 ! 148 ALLOCATE( ub2_b(jpi,jpj), vb2_b(jpi,jpj), un_bf(jpi,jpj), vn_bf(jpi,jpj) , STAT=ierr(5) ) 149 150 #if defined key_agrif 151 ALLOCATE( ub2_i_b(jpi,jpj), vb2_i_b(jpi,jpj) , STAT=ierr(5) ) 152 #endif 153 ALLOCATE( ht_0_xtd (idbi:idei,idbj:idej), hu_0_xtd (idbi:idei,idbj:idej), hv_0_xtd (idbi:idei,idbj:idej) & 154 & , r1_e1e2t_xtd(idbi:idei,idbj:idej), r1_e1e2u_xtd(idbi:idei,idbj:idej), r1_e1e2v_xtd(idbi:idei,idbj:idej) & 155 & , ssmask_xtd (idbi:idei,idbj:idej), ssumask_xtd (idbi:idei,idbj:idej), ssvmask_xtd (idbi:idei,idbj:idej) & 156 & , e1e2t_xtd (idbi:idei,idbj:idej), e2u_xtd (idbi:idei,idbj:idej), e1v_xtd (idbi:idei,idbj:idej) & 157 & , r1_e1u_xtd (idbi:idei,idbj:idej), r1_e2v_xtd (idbi:idei,idbj:idej) & 158 & , STAT=ierr(6) ) 112 159 ! 113 160 dyn_spg_ts_alloc = MAXVAL( ierr(:) ) … … 146 193 INTEGER, INTENT(in) :: kt ! ocean time-step index 147 194 ! 148 INTEGER :: ji, jj, jk, j n! dummy loop indices195 INTEGER :: ji, jj, jk, jm ! dummy loop indices 149 196 LOGICAL :: ll_fw_start ! =T : forward integration 150 197 LOGICAL :: ll_init ! =T : special startup of 2d equations … … 155 202 REAL(wp) :: zhu_bck, zhv_bck, zhdiv ! - - 156 203 REAL(wp) :: zun_save, zvn_save ! - - 157 REAL(wp), DIMENSION(jpi,jpj) :: zu_trd, zu_frc, zu_spg, zssh_frc 158 REAL(wp), DIMENSION(jpi,jpj) :: zv_trd, zv_frc, zv_spg 159 REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhtp2_e 160 REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zsshp2_e 161 REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points 162 REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV ! fluxes 204 ! 205 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zu_trd, zssh_frc 206 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zv_trd 207 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zhU , zhV 208 ! 209 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: hu_n_xtd, hv_n_xtd 210 ! 211 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zu_spg , zv_spg 212 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zu_frc , zv_frc 213 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsshu_a, zhup2_e, zhtp2_e 214 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsshv_a, zhvp2_e, zsshp2_e 215 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zCdU_u , zCdU_v ! top/bottom stress at u- & v-points 216 163 217 ! 164 218 REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True. … … 170 224 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points 171 225 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2 ! averages over the sub-steps of zuwdmask and zvwdmask 172 !!---------------------------------------------------------------------- 173 ! 226 227 228 INTEGER :: idbi, idei, idbj, idej ! lower/upper bounds of extended arrays 229 INTEGER :: ixtd ! number of halos over which the solution is currently correct 230 INTEGER :: ibi, iei, ibj, iej ! lower and upper bounds over which the solution is currently correct 231 !!---------------------------------------------------------------------- 232 ! 233 234 174 235 IF( ln_wd_il ) ALLOCATE( zcpx(jpi,jpj), zcpy(jpi,jpj) ) 175 236 ! !* Allocate temporary arrays 176 IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj)) 237 IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj) ) 238 239 idbi = 1 - nn_hlts ; idbj = 1 - nn_hlts 240 idei = jpi + nn_hlts ; idej = jpj + nn_hlts 241 ! ! allocate local arrays 242 ALLOCATE( zu_spg (idbi:idei,idbj:idej), zv_spg (idbi:idei,idbj:idej) & 243 & , zsshu_a (idbi:idei,idbj:idej), zsshv_a (idbi:idei,idbj:idej) & 244 & , zhup2_e (idbi:idei,idbj:idej), zhvp2_e (idbi:idei,idbj:idej) & 245 & , zCdU_u (idbi:idei,idbj:idej), zCdU_v (idbi:idei,idbj:idej) & 246 & , zhtp2_e (idbi:idei,idbj:idej), zsshp2_e(idbi:idei,idbj:idej) & 247 & , zu_trd (idbi:idei,idbj:idej), zu_frc (idbi:idei,idbj:idej) & 248 & , zv_trd (idbi:idei,idbj:idej), zv_frc (idbi:idei,idbj:idej) & 249 & , zhU (idbi:idei,idbj:idej), zhV (idbi:idei,idbj:idej) & 250 & , zssh_frc(idbi:idei,idbj:idej) ) 251 ! ! allocate redundant arrays 252 ALLOCATE( hu_n_xtd(idbi:idei,idbj:idej), hv_n_xtd(idbi:idei,idbj:idej) ) 177 253 ! 178 254 zmdi=1.e+20 ! missing data indicator for masking … … 227 303 ! != zu_frc = 1/H e3*d/dt(Ua) =! (Vertical mean of Ua, the 3D trends) 228 304 ! ! --------------------------- ! 229 zu_frc( :,:) = SUM( e3u_n(:,:,:) * ua(:,:,:) * umask(:,:,:) , DIM=3 ) * r1_hu_n(:,:)230 zv_frc( :,:) = SUM( e3v_n(:,:,:) * va(:,:,:) * vmask(:,:,:) , DIM=3 ) * r1_hv_n(:,:)305 zu_frc(1:jpi,1:jpj) = SUM( e3u_n(:,:,:) * ua(:,:,:) * umask(:,:,:) , DIM=3 ) * r1_hu_n(:,:) 306 zv_frc(1:jpi,1:jpj) = SUM( e3v_n(:,:,:) * va(:,:,:) * vmask(:,:,:) , DIM=3 ) * r1_hv_n(:,:) 231 307 ! 232 308 ! 233 309 ! != Ua => baroclinic trend =! (remove its vertical mean) 234 310 DO jk = 1, jpkm1 ! ------------------------ ! 235 ua(:,:,jk) = ( ua(:,:,jk) - zu_frc( :,:) ) * umask(:,:,jk)236 va(:,:,jk) = ( va(:,:,jk) - zv_frc( :,:) ) * vmask(:,:,jk)311 ua(:,:,jk) = ( ua(:,:,jk) - zu_frc(1:jpi,1:jpj) ) * umask(:,:,jk) 312 va(:,:,jk) = ( va(:,:,jk) - zv_frc(1:jpi,1:jpj) ) * vmask(:,:,jk) 237 313 END DO 238 314 … … 243 319 ! ! ------------------------------------------------- ! 244 320 ! 245 IF( kt == nit000 .OR. .NOT. ln_linssh ) CALL dyn_cor_2D_init ! Set zwz, the barotropic Coriolis force coefficient 246 ! ! recompute zwz = f/depth at every time step for (.NOT.ln_linssh) as the water colomn height changes 321 ! Set zwz, the barotropic Coriolis force coefficient 322 ! recompute zwz = f/depth at every time step for (.NOT.ln_linssh) as the water colomn height changes 323 IF( kt == nit000 .OR. .NOT. ln_linssh ) CALL dyn_cor_2d_init( idbi, idei, idbj, idej ) 247 324 ! 248 325 ! !* 2D Coriolis trends 249 zhU(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:) ! now fluxes 250 zhV(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) ! NB: FULL domain : put a value in last row and column 251 ! 252 CALL dyn_cor_2d( hu_n, hv_n, un_b, vn_b, zhU, zhV, & ! <<== in 253 & zu_trd, zv_trd ) ! ==>> out 326 zhU(1:jpi,1:jpj) = un_b(:,:) * hu_n(:,:) * e2u(:,:) ! now fluxes 327 zhV(1:jpi,1:jpj) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) ! NB: FULL domain : put a value in last row and column 328 ! 329 ! ! ht_n, hu_n, hv_n, un_b, vn_b are of size 1:jpi 1:jpj 330 ! ! zhU, zhV, zu_trd, zv_trd are of size idbi:idei idbj:idej 331 CALL dyn_cor_2d( ht_n, hu_n, hv_n, un_b, vn_b, zhU, zhV, 1 , jpi , 1 , jpj & ! <<== in 332 & , zu_trd, zv_trd, idbi, idei, idbj, idej ) ! ==>> out 254 333 ! 255 334 IF( .NOT.ln_linssh ) THEN !* surface pressure gradient (variable volume only) … … 285 364 ! != Add bottom stress contribution from baroclinic velocities =! 286 365 ! ! ----------------------------------------------------------- ! 287 CALL dyn_drg_init( zu_frc, zv_frc, zCdU_u, zCdU_v ) ! also provide the barotropic drag coefficients 366 CALL drg_init( idbi, idei, idbj, idej, zu_frc, zv_frc, zCdU_u, zCdU_v ) ! also provide the barotropic drag coefficients 367 ! ! arrays are computed on inner domain 288 368 ! 289 369 ! != Add atmospheric pressure forcing =! … … 335 415 ! ! --------------------------------------------------- ! 336 416 IF (ln_bt_fw) THEN ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2) 337 zssh_frc( :,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )417 zssh_frc(1:jpi,1:jpj) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) 338 418 ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW) 339 419 zztmp = r1_rau0 * r1_2 340 zssh_frc( :,:) = zztmp * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) + fwfisf(:,:) + fwfisf_b(:,:) )420 zssh_frc(1:jpi,1:jpj) = zztmp * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) + fwfisf(:,:) + fwfisf_b(:,:) ) 341 421 ENDIF 342 422 ! != Add Stokes drift divergence =! (if exist) 343 423 IF( ln_sdw ) THEN ! ----------------------------- ! 344 zssh_frc( :,:) = zssh_frc(:,:) + div_sd(:,:)424 zssh_frc(1:jpi,1:jpj) = zssh_frc(1:jpi,1:jpj) + div_sd(:,:) 345 425 ENDIF 346 426 ! … … 349 429 ! ! ------------------------------------ ! 350 430 IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 351 zssh_frc( :,:) = zssh_frc(:,:) - ssh_iau(:,:)431 zssh_frc(1:jpi,1:jpj) = zssh_frc(1:jpi,1:jpj) - ssh_iau(:,:) 352 432 ENDIF 353 433 #endif … … 357 437 IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt ) 358 438 #endif 439 440 359 441 ! 360 442 ! ----------------------------------------------------------------------- … … 366 448 ! ! ==================== ! 367 449 ! Initialize barotropic variables: 368 IF( ll_init ) THEN450 IF( ll_init ) THEN 369 451 sshbb_e(:,:) = 0._wp 370 452 ubb_e (:,:) = 0._wp … … 376 458 ! 377 459 IF( ln_linssh ) THEN ! mid-step ocean depth is fixed (hup2_e=hu_n=hu_0) 378 zh up2_e(:,:) = hu_n(:,:)379 zh vp2_e(:,:) = hv_n(:,:)380 zh tp2_e(:,:) = ht_n(:,:)460 zhtp2_e(1:jpi,1:jpj) = ht_n(1:jpi,1:jpj) 461 zhup2_e(1:jpi,1:jpj) = hu_n(1:jpi,1:jpj) 462 zhvp2_e(1:jpi,1:jpj) = hv_n(1:jpi,1:jpj) 381 463 ENDIF 382 464 ! 383 465 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 384 sshn_e( :,:) = sshn(:,:)385 un_e ( :,:) = un_b(:,:)386 vn_e ( :,:) = vn_b(:,:)387 ! 388 hu_e ( :,:) = hu_n(:,:)389 hv_e ( :,:) = hv_n(:,:)390 hur_e ( :,:) = r1_hu_n(:,:)391 hvr_e ( :,:) = r1_hv_n(:,:)466 sshn_e(1:jpi,1:jpj) = sshn(1:jpi,1:jpj) 467 un_e (1:jpi,1:jpj) = un_b(1:jpi,1:jpj) 468 vn_e (1:jpi,1:jpj) = vn_b(1:jpi,1:jpj) 469 ! 470 hu_e (1:jpi,1:jpj) = hu_n(1:jpi,1:jpj) 471 hv_e (1:jpi,1:jpj) = hv_n(1:jpi,1:jpj) 472 hur_e (1:jpi,1:jpj) = r1_hu_n(1:jpi,1:jpj) 473 hvr_e (1:jpi,1:jpj) = r1_hv_n(1:jpi,1:jpj) 392 474 ELSE ! CENTRED integration: start from BEFORE fields 393 sshn_e(:,:) = sshb(:,:) 394 un_e (:,:) = ub_b(:,:) 395 vn_e (:,:) = vb_b(:,:) 396 ! 397 hu_e (:,:) = hu_b(:,:) 398 hv_e (:,:) = hv_b(:,:) 399 hur_e (:,:) = r1_hu_b(:,:) 400 hvr_e (:,:) = r1_hv_b(:,:) 401 ENDIF 475 sshn_e(1:jpi,1:jpj) = sshb(1:jpi,1:jpj) 476 un_e (1:jpi,1:jpj) = ub_b(1:jpi,1:jpj) 477 vn_e (1:jpi,1:jpj) = vb_b(1:jpi,1:jpj) 478 ! 479 hu_e (1:jpi,1:jpj) = hu_b(1:jpi,1:jpj) 480 hv_e (1:jpi,1:jpj) = hv_b(1:jpi,1:jpj) 481 hur_e (1:jpi,1:jpj) = r1_hu_b(1:jpi,1:jpj) 482 hvr_e (1:jpi,1:jpj) = r1_hv_b(1:jpi,1:jpj) 483 ENDIF 484 ! 485 hu_n_xtd(1:jpi,1:jpj) = hu_n(1:jpi,1:jpj) 486 hv_n_xtd(1:jpi,1:jpj) = hv_n(1:jpi,1:jpj) 487 ! 488 ! 489 ! ! Extend arrays 490 ! ! -------------- 491 ! 492 IF( ln_linssh ) THEN 493 CALL lbc_lnk_multi( 'dynspg_ts', hu_n_xtd, 'U', -1._wp, hv_n_xtd, 'V', -1._wp & 494 & , zCdU_u , 'U', -1._wp, zCdU_v , 'V', -1._wp & 495 & , zu_frc , 'U', -1._wp, zv_frc , 'V', -1._wp & 496 & , un_e , 'U', -1._wp, vn_e , 'V', -1._wp & 497 & , hu_e , 'U', -1._wp, hv_e , 'V', -1._wp & 498 & , hur_e , 'U', -1._wp, hvr_e , 'V', -1._wp & 499 & , zhtp2_e , 'T', 1._wp, zhup2_e, 'U', -1._wp, zhvp2_e , 'V', -1._wp & 500 & , zssh_frc, 'T', 1._wp, sshn_e , 'T', 1._wp, khlcom = nn_hls+nn_hlts ) 501 ELSE 502 CALL lbc_lnk_multi( 'dynspg_ts', hu_n_xtd, 'U', -1._wp, hv_n_xtd, 'V', -1._wp & 503 & , zCdU_u , 'U', -1._wp, zCdU_v , 'V', -1._wp & 504 & , zu_frc , 'U', -1._wp, zv_frc , 'V', -1._wp & 505 & , un_e , 'U', -1._wp, vn_e , 'V', -1._wp & 506 & , hu_e , 'U', -1._wp, hv_e , 'V', -1._wp & 507 & , hur_e , 'U', -1._wp, hvr_e , 'V', -1._wp & 508 & , zssh_frc, 'T', 1._wp, sshn_e , 'T', 1._wp, khlcom = nn_hls+nn_hlts ) 509 END IF 402 510 ! 403 511 ! Initialize sums: … … 413 521 zuwdav2 (:,:) = 0._wp 414 522 zvwdav2 (:,:) = 0._wp 415 END IF 416 523 END IF 524 525 ixtd = nn_hls + nn_hlts ! solution is now correct over the whole domain (interior + regular halos + time splitting halos) 526 ibi = 1 - nn_hlts ; ibj = 1 - nn_hlts 527 iei = jpi + nn_hlts ; iej = jpj + nn_hlts 417 528 ! ! ==================== ! 418 DO j n= 1, icycle ! sub-time-step loop !529 DO jm = 1, icycle ! sub-time-step loop ! 419 530 ! ! ==================== ! 420 531 ! 421 l_full_nf_update = j n == icycle ! false: disable full North fold update (performances) for jn= 1 to icycle-1532 l_full_nf_update = jm == icycle ! false: disable full North fold update (performances) for jm = 1 to icycle-1 422 533 ! 423 534 ! !== Update the forcing ==! (BDY and tides) 424 535 ! 425 IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=j n, kt_offset= noffset+1 )426 IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=j n, kt_offset= noffset )427 ! 428 ! !== extrapolation at mid-step ==! (j n+1/2)536 IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jm, kt_offset= noffset+1 ) 537 IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jm, kt_offset= noffset ) 538 ! 539 ! !== extrapolation at mid-step ==! (jm+1/2) 429 540 ! 430 541 ! !* Set extrapolation coefficients for predictor step: 431 IF ((jn<3).AND.ll_init) THEN ! Forward542 IF( (jm<3) .AND. ll_init ) THEN ! Forward 432 543 za1 = 1._wp 433 544 za2 = 0._wp … … 439 550 ENDIF 440 551 ! 441 ! !* Extrapolate barotropic velocities at mid-step (j n+1/2)552 ! !* Extrapolate barotropic velocities at mid-step (jm+1/2) 442 553 !-- m+1/2 m m-1 m-2 --! 443 554 !-- u = (3/2+beta) u -(1/2+2beta) u + beta u --! 444 555 !-------------------------------------------------------------------------! 445 ua_e( :,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)446 va_e( :,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)556 ua_e(ibi:iei,ibj:iej) = za1 * un_e(ibi:iei,ibj:iej) + za2 * ub_e(ibi:iei,ibj:iej) + za3 * ubb_e(ibi:iei,ibj:iej) 557 va_e(ibi:iei,ibj:iej) = za1 * vn_e(ibi:iei,ibj:iej) + za2 * vb_e(ibi:iei,ibj:iej) + za3 * vbb_e(ibi:iei,ibj:iej) 447 558 448 559 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 449 560 ! ! ------------------ 450 ! Extrapolate Sea Level at step jit+0.5:561 ! !* Extrapolate Sea Level at mid-step (jm+1/2) 451 562 !-- m+1/2 m m-1 m-2 --! 452 563 !-- ssh = (3/2+beta) ssh -(1/2+2beta) ssh + beta ssh --! 453 564 !--------------------------------------------------------------------------------! 454 zsshp2_e( :,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)455 565 zsshp2_e(ibi:iei,ibj:iej) = za1 * sshn_e (ibi:iei,ibj:iej) + za2 * sshb_e(ibi:iei,ibj:iej) & 566 & + za3 * sshbb_e(ibi:iei,ibj:iej) 456 567 ! set wetting & drying mask at tracer points for this barotropic mid-step 457 568 IF( ln_wd_dl ) CALL wad_tmsk( zsshp2_e, ztwdmask ) 458 569 ! 459 570 ! ! ocean t-depth at mid-step 460 zhtp2_e( :,:) = ht_0(:,:) + zsshp2_e(:,:)571 zhtp2_e(ibi:iei,ibj:iej) = ht_0_xtd(ibi:iei,ibj:iej) + zsshp2_e(ibi:iei,ibj:iej) 461 572 ! 462 ! ! ocean u- and v-depth at mid-step (separate DO-loops remove the need of a lbc_lnk) 463 DO jj = 1, jpj 464 DO ji = 1, jpim1 ! not jpi-column 465 zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj) & 466 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 467 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 468 END DO 469 END DO 470 DO jj = 1, jpjm1 ! not jpj-row 471 DO ji = 1, jpi 472 zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj) & 473 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 474 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 573 ! ! ocean u- and v-depth at mid-step 574 DO jj = ibj, iej-1 ! not last column, not last row 575 DO ji = ibi, iei-1 576 zhup2_e(ji,jj) = hu_0_xtd(ji,jj) + r1_2 * r1_e1e2u_xtd(ji,jj) & 577 & * ( e1e2t_xtd(ji ,jj) * zsshp2_e(ji ,jj) & 578 & + e1e2t_xtd(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask_xtd(ji,jj) 579 zhvp2_e(ji,jj) = hv_0_xtd(ji,jj) + r1_2 * r1_e1e2v_xtd(ji,jj) & 580 & * ( e1e2t_xtd(ji,jj ) * zsshp2_e(ji,jj ) & 581 & + e1e2t_xtd(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask_xtd(ji,jj) 475 582 END DO 476 583 END DO … … 478 585 ENDIF 479 586 ! 480 ! !== after SSH ==! (j n+1)587 ! !== after SSH ==! (jm+1) 481 588 ! 482 589 ! ! update (ua_e,va_e) to enforce volume conservation at open boundaries 483 590 ! ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d 484 IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, j n, ua_e, va_e, zhup2_e, zhvp2_e )591 IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jm, ua_e, va_e, zhup2_e, zhvp2_e ) 485 592 ! 486 593 ! ! resulting flux at mid-step (not over the full domain) 487 zhU( 1:jpim1,1:jpj ) = e2u(1:jpim1,1:jpj ) * ua_e(1:jpim1,1:jpj ) * zhup2_e(1:jpim1,1:jpj ) ! not jpi-column488 zhV( 1:jpi ,1:jpjm1) = e1v(1:jpi ,1:jpjm1) * va_e(1:jpi ,1:jpjm1) * zhvp2_e(1:jpi ,1:jpjm1) ! not jpj-row594 zhU(ibi:iei-1,ibj:iej-1) = e2u_xtd(ibi:iei-1,ibj:iej-1) * ua_e(ibi:iei-1,ibj:iej-1) * zhup2_e(ibi:iei-1,ibj:iej-1) 595 zhV(ibi:iei-1,ibj:iej-1) = e1v_xtd(ibi:iei-1,ibj:iej-1) * va_e(ibi:iei-1,ibj:iej-1) * zhvp2_e(ibi:iei-1,ibj:iej-1) 489 596 ! 490 597 #if defined key_agrif … … 524 631 ! 525 632 ENDIF 526 !527 !528 ! Compute Sea Level at step jit+1529 !-- m+1 m m+1/2 --!530 !-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --!531 !-------------------------------------------------------------------------!532 DO jj = 2, jpjm1 ! INNER domain533 DO ji = 2, jpim1534 zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj)535 ssha_e(ji,jj) = ( sshn_e(ji,jj) - rdtbt * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj)536 END DO537 END DO538 !539 CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp )540 !541 ! ! Sum over sub-time-steps to compute advective velocities542 za2 = wgtbtp2(jn) ! zhU, zhV hold fluxes extrapolated at jn+0.5543 un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:)544 vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:)545 633 ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True) 546 634 IF ( ln_wd_dl_bc ) THEN … … 548 636 zvwdav2(1:jpi ,1:jpjm1) = zvwdav2(1:jpi ,1:jpjm1) + za2 * zvwdmask(1:jpi ,1:jpjm1) ! not jpj-row 549 637 END IF 638 ! 639 ! Sum over sub-time-steps to compute advective velocities (only correct on interior domain) 640 za2 = wgtbtp2(jm) ! zhU, zhV hold fluxes extrapolated at jm+1/2 641 un_adv(1:jpi,1:jpj) = un_adv(1:jpi,1:jpj) + za2 * zhU(1:jpi,1:jpj) * r1_e2u(1:jpi,1:jpj) 642 vn_adv(1:jpi,1:jpj) = vn_adv(1:jpi,1:jpj) + za2 * zhV(1:jpi,1:jpj) * r1_e1v(1:jpi,1:jpj) 643 ! 644 ! 645 ! Compute Sea Level at step jit+1 646 !-- m+1 m m+1/2 --! 647 !-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --! 648 !-------------------------------------------------------------------------! 649 ! correct domain reduction 650 ixtd = ixtd - 1 651 ibi = ibi + 1 ; ibj = ibj + 1 652 iei = iei - 1 ; iej = iej - 1 653 DO jj = ibj, iej 654 DO ji = ibi, iei 655 zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t_xtd(ji,jj) 656 ssha_e(ji,jj) = ( sshn_e(ji,jj) - rdtbt * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask_xtd(ji,jj) 657 END DO 658 END DO 659 ! 660 IF( nn_hlts == 0 ) THEN 661 CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp ) 662 ixtd = nn_hls + nn_hlts ! solution is now correct over the whole domain 663 ibi = 1 - nn_hlts ; ibj = 1 - nn_hlts 664 iei = jpi + nn_hlts ; iej = jpj + nn_hlts 665 END IF 666 550 667 ! 551 668 ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T) 552 IF( ln_bdy ) CALL bdy_ssh( ssha_e ) 669 IF( ln_bdy ) THEN 670 CALL swap_bdyptr ! bdy treatment is now done on extended domain 671 CALL bdy_ssh( ssha_e, idbi, idei, idbj, idej, ldcomall=.true., pmask=ssmask_xtd, khlcom=nn_hls+nn_hlts ) 672 CALL swap_bdyptr ! bdy treatment is now done on regular domain 673 END IF 674 553 675 #if defined key_agrif 554 IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( j n)676 IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jm ) 555 677 #endif 556 !557 ! Sea Surface Height at u-,v-points (vvl case only)558 IF( .NOT.ln_linssh ) THEN559 DO jj = 2, jpjm1 ! INNER domain, will be extended to whole domain later560 DO ji = 2, jpim1 ! NO Vector Opt.561 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &562 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &563 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) )564 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &565 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &566 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) )567 END DO568 END DO569 ENDIF570 678 ! 571 679 ! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2 572 !-- m+1/2 m+1 m m-1 m-2 --! 573 !-- ssh' = za0 * ssh + za1 * ssh + za2 * ssh + za3 * ssh --! 574 !------------------------------------------------------------------------------------------! 575 CALL ts_bck_interp( jn, ll_init, za0, za1, za2, za3 ) ! coeficients of the interpolation 576 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 577 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 680 !-- m+1/2 m+1 m m-1 m-2 --! 681 !-- ssh' = za0 * ssh + za1 * ssh + za2 * ssh + za3 * ssh --! 682 !-----------------------------------------------------------------------------------------! 683 CALL ts_bck_interp( jm, ll_init, za0, za1, za2, za3 ) ! coeficients of the interpolation 684 zsshp2_e(ibi:iei,ibj:iej) = za0 * ssha_e(ibi:iei,ibj:iej) + za1 * sshn_e (ibi:iei,ibj:iej) & 685 & + za2 * sshb_e(ibi:iei,ibj:iej) + za3 * sshbb_e(ibi:iei,ibj:iej) 686 ! 687 ! 688 ! Sea Surface Height at u-,v-points (vvl case only) 689 IF( .NOT.ln_linssh ) THEN 690 DO jj = ibj, iej-1 691 DO ji = ibi, iei-1 692 zsshu_a(ji,jj) = r1_2 * ssumask_xtd(ji,jj) * r1_e1e2u_xtd(ji,jj) & 693 & * ( e1e2t_xtd(ji ,jj ) * ssha_e(ji ,jj ) & 694 & + e1e2t_xtd(ji+1,jj ) * ssha_e(ji+1,jj ) ) 695 zsshv_a(ji,jj) = r1_2 * ssvmask_xtd(ji,jj) * r1_e1e2v_xtd(ji,jj) & 696 & * ( e1e2t_xtd(ji ,jj ) * ssha_e(ji ,jj ) & 697 & + e1e2t_xtd(ji ,jj+1) * ssha_e(ji ,jj+1) ) 698 END DO 699 END DO 700 ENDIF 578 701 ! 579 702 ! ! Surface pressure gradient 580 703 zldg = ( 1._wp - rn_scal_load ) * grav ! local factor 581 DO jj = 2, jpjm1582 DO ji = 2, jpim1583 zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u (ji,jj)584 zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v (ji,jj)704 DO jj = ibj, iej-1 705 DO ji = ibi, iei-1 706 zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u_xtd(ji,jj) 707 zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v_xtd(ji,jj) 585 708 END DO 586 709 END DO … … 592 715 ! 593 716 ! Add Coriolis trend: 594 ! zwz array belowor triads normally depend on sea level with ln_linssh=F and should be updated717 ! - zwz array used in dyn_cor_2d or triads normally depend on sea level with ln_linssh=F and should be updated 595 718 ! at each time step. We however keep them constant here for optimization. 596 ! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated) 597 CALL dyn_cor_2d( zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd ) 719 ! - Recall that zhU and zhV hold fluxes at jm+1/2 (extrapolated not backward interpolated) 720 ! - zu_trd_xtd and zv_trd_xtd are only correct on (ibi+1:iei-1,ibj+1:iej-1) 721 ! NOTE : input flux arguments have to be correct (ibi:iei,ibj:iej) -> a lbc call between input arguments computation 722 ! and this call without fluxes (typically after ssh at step m+1 computation) would not yield correct results 723 CALL dyn_cor_2d( zhtp2_e, zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV , idbi, idei, idbj, idej & 724 & , zu_trd, zv_trd , idbi, idei, idbj, idej ) 598 725 ! 599 726 ! Add tidal astronomical forcing if defined 727 ! pot_astro is correct on 1:jpi,1:jpj 600 728 IF ( ln_tide .AND. ln_tide_pot ) THEN 601 729 DO jj = 2, jpjm1 … … 610 738 !jth do implicitly instead 611 739 IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs 612 DO jj = 2, jpjm1613 DO ji = fs_2, fs_jpim1 ! vector opt.740 DO jj = ibj, iej 741 DO ji = ibi, iei 614 742 zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj) 615 743 zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj) … … 624 752 !-- u = u + delta_t' * \ (1-r)*g * grad_x( ssh') - f * k vect u + frc / --! 625 753 !-- --! 626 !-- FLUX FORM--!754 !-- FLUX FORM --! 627 755 !-- m+1 __1__ / m m / m+1/2 m+1/2 m+1/2 n \ \ --! 628 756 !-- u = m+1 | h * u + delta_t' * \ h * (1-r)*g * grad_x( ssh') - h * f * k vect u + h * frc / | --! 629 757 !-- h \ / --! 630 758 !------------------------------------------------------------------------------------------------------------------------! 759 ! correct domain reduction 760 ixtd = ixtd - 1 761 ibi = ibi + 1 ; ibj = ibj + 1 762 iei = iei - 1 ; iej = iej - 1 763 ! 631 764 IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form 632 DO jj = 2, jpjm1633 DO ji = fs_2, fs_jpim1 ! vector opt.765 DO jj = ibj, iej 766 DO ji = ibi, iei 634 767 ua_e(ji,jj) = ( un_e(ji,jj) & 635 768 & + rdtbt * ( zu_spg(ji,jj) & 636 769 & + zu_trd(ji,jj) & 637 770 & + zu_frc(ji,jj) ) & 638 & ) * ssumask (ji,jj)771 & ) * ssumask_xtd(ji,jj) 639 772 640 773 va_e(ji,jj) = ( vn_e(ji,jj) & … … 642 775 & + zv_trd(ji,jj) & 643 776 & + zv_frc(ji,jj) ) & 644 & ) * ssvmask (ji,jj)777 & ) * ssvmask_xtd(ji,jj) 645 778 END DO 646 779 END DO 647 780 ! 648 781 ELSE !* Flux form 649 DO jj = 2, jpjm1 650 DO ji = 2, jpim1 651 ! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2 652 ! ! backward interpolated depth used in spg terms at jn+1/2 653 zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 654 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 655 zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 656 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 657 ! ! inverse depth at jn+1 658 z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) ) 659 z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) ) 782 DO jj = ibj, iej 783 DO ji = ibi, iei 784 ! ! hu_e, hv_e hold depth at jm, zhup2_e, zhvp2_e hold extrapolated depth at jm+1/2 785 ! ! backward interpolated depth used in spg terms at jm+1/2 786 zhu_bck = hu_0_xtd(ji,jj) + r1_2*r1_e1e2u_xtd(ji,jj) * & 787 & ( e1e2t_xtd(ji ,jj) * zsshp2_e(ji ,jj) & 788 & + e1e2t_xtd(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask_xtd(ji,jj) 789 zhv_bck = hv_0_xtd(ji,jj) + r1_2*r1_e1e2v_xtd(ji,jj) * & 790 & ( e1e2t_xtd(ji,jj ) * zsshp2_e(ji,jj ) & 791 & + e1e2t_xtd(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask_xtd(ji,jj) 792 ! ! inverse depth at jm+1 793 z1_hu = ssumask_xtd(ji,jj) / ( hu_0_xtd(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask_xtd(ji,jj) ) 794 z1_hv = ssvmask_xtd(ji,jj) / ( hv_0_xtd(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask_xtd(ji,jj) ) 660 795 ! 661 ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj)&662 & + rdtbt * ( zhu_bck * zu_spg(ji,jj) & !663 & + zhup2_e(ji,jj) * zu_trd 664 & + hu_n (ji,jj) * zu_frc(ji,jj) ) ) * z1_hu796 ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) & 797 & + rdtbt * ( zhu_bck * zu_spg(ji,jj) & ! 798 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & ! 799 & + hu_n_xtd (ji,jj) * zu_frc(ji,jj) ) ) * z1_hu 665 800 ! 666 va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj)&667 & + rdtbt * ( zhv_bck * zv_spg(ji,jj) & !668 & + zhvp2_e(ji,jj) * zv_trd 669 & + hv_n (ji,jj) * zv_frc(ji,jj) ) ) * z1_hv801 va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) & 802 & + rdtbt * ( zhv_bck * zv_spg(ji,jj) & ! 803 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & ! 804 & + hv_n_xtd (ji,jj) * zv_frc(ji,jj) ) ) * z1_hv 670 805 END DO 671 806 END DO … … 681 816 ENDIF 682 817 683 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 684 hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1) 685 hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) ) 686 hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1) 687 hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) ) 688 CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp & 689 & , hu_e , 'U', -1._wp, hv_e , 'V', -1._wp & 690 & , hur_e, 'U', -1._wp, hvr_e, 'V', -1._wp ) 691 ELSE 692 CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp ) 693 ENDIF 694 ! 695 ! 696 ! ! open boundaries 697 IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 818 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) and inverse depth 819 hu_e (ibi:iei,ibj:iej) = hu_0_xtd(ibi:iei,ibj:iej) + zsshu_a(ibi:iei,ibj:iej) 820 hv_e (ibi:iei,ibj:iej) = hv_0_xtd(ibi:iei,ibj:iej) + zsshv_a(ibi:iei,ibj:iej) 821 hur_e(ibi:iei,ibj:iej) = ssumask_xtd(ibi:iei,ibj:iej) / ( hu_e(ibi:iei,ibj:iej) + 1._wp - ssumask_xtd(ibi:iei,ibj:iej) ) 822 hvr_e(ibi:iei,ibj:iej) = ssvmask_xtd(ibi:iei,ibj:iej) / ( hv_e(ibi:iei,ibj:iej) + 1._wp - ssvmask_xtd(ibi:iei,ibj:iej) ) 823 ENDIF 824 825 IF( ixtd == 0 ) THEN 826 IF( .NOT. ln_linssh ) THEN 827 CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp & ! after 828 & , un_e , 'U', -1._wp, vn_e , 'V', -1._wp & ! now 829 & , ub_e , 'U', -1._wp, vb_e , 'V', -1._wp & ! before 830 & , ubb_e , 'U', -1._wp, vbb_e , 'V', -1._wp & ! before before 831 & , ssha_e, 'T', 1._wp, sshn_e , 'T', 1._wp & ! after, now 832 & , sshb_e, 'T', 1._wp, sshbb_e, 'T', 1._wp & ! before, before before 833 & , hu_e , 'U', -1._wp, hv_e , 'V', -1._wp & 834 & , hur_e , 'U', -1._wp, hvr_e , 'V', -1._wp & 835 & , khlcom = nn_hls+nn_hlts ) 836 ELSE 837 CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp & ! after 838 & , un_e , 'U', -1._wp, vn_e , 'V', -1._wp & ! now 839 & , ub_e , 'U', -1._wp, vb_e , 'V', -1._wp & ! before 840 & , ubb_e , 'U', -1._wp, vbb_e , 'V', -1._wp & ! before before 841 & , ssha_e, 'T', 1._wp, sshn_e , 'T', 1._wp & ! after, now 842 & , sshb_e, 'T', 1._wp, sshbb_e, 'T', 1._wp & ! before, before before 843 & , khlcom = nn_hls+nn_hlts ) 844 END IF 845 ixtd = nn_hls + nn_hlts ! solution is now correct over the whole domain 846 ibi = 1 - nn_hlts ; ibj = 1 - nn_hlts 847 iei = jpi + nn_hlts ; iej = jpj + nn_hlts 848 END IF 849 ! 850 ! 851 ! ! open boundaries 852 ! ! bdy treatment is here done on regular domain (nn_hlts forced to 1 if ln_bdy or ln_tides) 853 IF( ln_bdy ) CALL bdy_dyn2d( jm, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e, idbi, idei, idbj, idej & 854 & , ldcomall=.true., pumask=ssumask_xtd, pvmask=ssvmask_xtd, khlcom=nn_hls+nn_hlts ) 698 855 #if defined key_agrif 699 IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( j n) ! Agrif856 IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jm ) ! Agrif 700 857 #endif 701 858 ! !* Swap … … 715 872 ! !* Sum over whole bt loop 716 873 ! ! ---------------------- 717 za1 = wgtbtp1(j n)874 za1 = wgtbtp1(jm) 718 875 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities 719 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:)720 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:)876 ua_b(1:jpi,1:jpj) = ua_b(1:jpi,1:jpj) + za1 * ua_e(1:jpi,1:jpj) 877 va_b(1:jpi,1:jpj) = va_b(1:jpi,1:jpj) + za1 * va_e(1:jpi,1:jpj) 721 878 ELSE ! Sum transports 722 879 IF ( .NOT.ln_wd_dl ) THEN 723 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:)724 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:)880 ua_b(1:jpi,1:jpj) = ua_b(1:jpi,1:jpj) + za1 * ua_e(1:jpi,1:jpj) * hu_e(1:jpi,1:jpj) 881 va_b(1:jpi,1:jpj) = va_b(1:jpi,1:jpj) + za1 * va_e(1:jpi,1:jpj) * hv_e(1:jpi,1:jpj) 725 882 ELSE 726 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:)727 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:)883 ua_b(1:jpi,1:jpj) = ua_b(1:jpi,1:jpj) + za1 * ua_e(1:jpi,1:jpj) * hu_e(1:jpi,1:jpj) * zuwdmask(1:jpi,1:jpj) 884 va_b(1:jpi,1:jpj) = va_b(1:jpi,1:jpj) + za1 * va_e(1:jpi,1:jpj) * hv_e(1:jpi,1:jpj) * zvwdmask(1:jpi,1:jpj) 728 885 END IF 729 886 ENDIF 730 887 ! ! Sum sea level 731 ssha( :,:) = ssha(:,:) + za1 * ssha_e(:,:)888 ssha(1:jpi,1:jpj) = ssha(1:jpi,1:jpj) + za1 * ssha_e(1:jpi,1:jpj) 732 889 733 890 ! ! ==================== ! … … 737 894 ! Phase 3. update the general trend with the barotropic trend 738 895 ! ----------------------------------------------------------------------------- 896 ! Correction on regular halos 897 CALL lbc_lnk_multi( 'dynspg_ts', un_adv, 'U', -1._wp, vn_adv, 'V', -1._wp & 898 & , ua_b , 'U', -1._wp, va_b , 'V', -1._wp & 899 & , ssha , 'T', -1._wp ) 739 900 ! 740 901 ! Set advection velocity correction: … … 783 944 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 784 945 END DO 785 END DO 786 CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions946 END DO ! Boundary conditions 947 CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp, khlcom=nn_hls+nn_hlts ) ! change array used? 787 948 ! 788 949 DO jk=1,jpkm1 … … 791 952 END DO 792 953 ! Save barotropic velocities not transport: 793 ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a( :,:) + 1._wp - ssumask(:,:) )794 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a( :,:) + 1._wp - ssvmask(:,:) )954 ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(1:jpi,1:jpj) + 1._wp - ssumask(:,:) ) 955 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(1:jpi,1:jpj) + 1._wp - ssvmask(:,:) ) 795 956 ENDIF 796 957 … … 841 1002 ENDIF 842 1003 ! 1004 1005 ! deallocate temporary arrays 1006 DEALLOCATE( zu_trd , zv_trd & 1007 & , zu_frc , zv_frc & 1008 & , zu_spg , zv_spg & 1009 & , zsshu_a , zsshv_a & 1010 & , zhup2_e , zhvp2_e & 1011 & , zCdU_u , zCdU_v & 1012 & , zhU , zhV & 1013 & , zssh_frc, zsshp2_e & 1014 & , zhtp2_e & 1015 & , hu_n_xtd, hv_n_xtd ) 1016 ! 843 1017 END SUBROUTINE dyn_spg_ts 844 1018 … … 850 1024 !! ** Purpose : Set time-splitting weights for temporal averaging (or not) 851 1025 !!---------------------------------------------------------------------- 852 LOGICAL, INTENT(in ) :: ll_av ! temporal averaging=.true.853 LOGICAL, INTENT(in ) :: ll_fw ! forward time splitting =.true.854 INTEGER, INTENT(inout) :: jpit! cycle length1026 LOGICAL, INTENT(in ) :: ll_av ! temporal averaging=.true. 1027 LOGICAL, INTENT(in ) :: ll_fw ! forward time splitting =.true. 1028 INTEGER, INTENT(inout) :: jpit ! cycle length 855 1029 REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) :: zwgt1, & ! Primary weights 856 1030 zwgt2 ! Secondary weights 857 1031 858 INTEGER :: jic, j n, ji ! temporary integers1032 INTEGER :: jic, jm, ji ! temporary integers 859 1033 REAL(wp) :: za1, za2 860 1034 !!---------------------------------------------------------------------- … … 880 1054 881 1055 CASE( 1 ) ! Boxcar, width = nn_baro 882 DO j n= 1, 3*nn_baro883 za1 = ABS(float(j n-jic))/float(nn_baro)1056 DO jm = 1, 3*nn_baro 1057 za1 = ABS(float(jm-jic))/float(nn_baro) 884 1058 IF (za1 < 0.5_wp) THEN 885 zwgt1(j n) = 1._wp886 jpit = j n1059 zwgt1(jm) = 1._wp 1060 jpit = jm 887 1061 ENDIF 888 1062 ENDDO 889 1063 890 1064 CASE( 2 ) ! Boxcar, width = 2 * nn_baro 891 DO j n= 1, 3*nn_baro892 za1 = ABS(float(j n-jic))/float(nn_baro)1065 DO jm = 1, 3*nn_baro 1066 za1 = ABS(float(jm-jic))/float(nn_baro) 893 1067 IF (za1 < 1._wp) THEN 894 zwgt1(j n) = 1._wp895 jpit = j n1068 zwgt1(jm) = 1._wp 1069 jpit = jm 896 1070 ENDIF 897 1071 ENDDO … … 905 1079 906 1080 ! Set secondary weights 907 DO j n= 1, jpit908 DO ji = j n, jpit909 zwgt2(j n) = zwgt2(jn) + zwgt1(ji)1081 DO jm = 1, jpit 1082 DO ji = jm, jpit 1083 zwgt2(jm) = zwgt2(jm) + zwgt1(ji) 910 1084 END DO 911 1085 END DO … … 914 1088 za1 = 1._wp / SUM(zwgt1(1:jpit)) 915 1089 za2 = 1._wp / SUM(zwgt2(1:jpit)) 916 DO j n= 1, jpit917 zwgt1(j n) = zwgt1(jn) * za1918 zwgt2(j n) = zwgt2(jn) * za21090 DO jm = 1, jpit 1091 zwgt1(jm) = zwgt1(jm) * za1 1092 zwgt2(jm) = zwgt2(jm) * za2 919 1093 END DO 920 1094 ! … … 1111 1285 ENDIF 1112 1286 ! 1287 ! initialize extended scale factors 1288 ht_0_xtd (1:jpi,1:jpj) = ht_0 (1:jpi,1:jpj) 1289 hu_0_xtd (1:jpi,1:jpj) = hu_0 (1:jpi,1:jpj) 1290 hv_0_xtd (1:jpi,1:jpj) = hv_0 (1:jpi,1:jpj) 1291 r1_e1e2t_xtd(1:jpi,1:jpj) = r1_e1e2t(1:jpi,1:jpj) 1292 r1_e1e2u_xtd(1:jpi,1:jpj) = r1_e1e2u(1:jpi,1:jpj) 1293 r1_e1e2v_xtd(1:jpi,1:jpj) = r1_e1e2v(1:jpi,1:jpj) 1294 e1e2t_xtd (1:jpi,1:jpj) = e1e2t (1:jpi,1:jpj) 1295 ssmask_xtd (1:jpi,1:jpj) = ssmask (1:jpi,1:jpj) 1296 ssumask_xtd (1:jpi,1:jpj) = ssumask (1:jpi,1:jpj) 1297 ssvmask_xtd (1:jpi,1:jpj) = ssvmask (1:jpi,1:jpj) 1298 e2u_xtd (1:jpi,1:jpj) = e2u (1:jpi,1:jpj) 1299 e1v_xtd (1:jpi,1:jpj) = e1v (1:jpi,1:jpj) 1300 r1_e1u_xtd (1:jpi,1:jpj) = r1_e1u (1:jpi,1:jpj) 1301 r1_e2v_xtd (1:jpi,1:jpj) = r1_e2v (1:jpi,1:jpj) 1302 ! 1303 CALL lbc_lnk_multi( 'dynspg_ts', ht_0_xtd , 'T', 1._wp, hu_0_xtd , 'U', -1._wp, hv_0_xtd , 'V', -1._wp & 1304 & , r1_e1e2t_xtd, 'T', 1._wp, r1_e1e2u_xtd, 'U', -1._wp, r1_e1e2v_xtd, 'V', -1._wp & 1305 & , ssmask_xtd , 'T', 1._wp, ssumask_xtd , 'U', -1._wp, ssvmask_xtd , 'V', -1._wp & 1306 & , e1e2t_xtd , 'T', 1._wp, e2u_xtd , 'U', -1._wp, e1v_xtd , 'V', -1._wp & 1307 & , r1_e1u_xtd , 'U', -1._wp, r1_e2v_xtd , 'V', -1._wp & 1308 & , khlcom = nn_hls+nn_hlts ) 1309 IF( ln_dynvor_enT ) THEN 1310 ff_t_xtd (1:jpi,1:jpj) = ff_t (1:jpi,1:jpj) 1311 CALL lbc_lnk_multi( 'dynspg_ts', ff_t_xtd , 'F', -1._wp, khlcom = nn_hls+nn_hlts ) 1312 END IF 1313 1314 ! 1113 1315 END SUBROUTINE dyn_spg_ts_init 1114 1316 1115 1317 1116 SUBROUTINE dyn_cor_2d_init 1318 SUBROUTINE dyn_cor_2d_init( kdbi, kdei, kdbj, kdej ) 1117 1319 !!--------------------------------------------------------------------- 1118 1320 !! *** ROUTINE dyn_cor_2d_init *** … … 1128 1330 !! Compute zwz = f / ( height of the water colomn ) 1129 1331 !!---------------------------------------------------------------------- 1332 INTEGER , INTENT(in ) :: kdbi, kdei, kdbj, kdej 1130 1333 INTEGER :: ji ,jj, jk ! dummy loop indices 1131 1334 REAL(wp) :: z1_ht … … 1137 1340 SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point 1138 1341 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 1139 DO jj = 1, jpj m11140 DO ji = 1, jpi m11141 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + 1142 & 1342 DO jj = 1, jpj-1 1343 DO ji = 1, jpi-1 1344 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 1345 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp 1143 1346 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1144 1347 END DO 1145 1348 END DO 1146 1349 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 1147 DO jj = 1, jpj m11148 DO ji = 1, jpi m11149 zwz(ji,jj) = ( ht_n (ji ,jj+1) + ht_n (ji+1,jj+1) &1350 DO jj = 1, jpj-1 1351 DO ji = 1, jpi-1 1352 zwz(ji,jj) = ( ht_n (ji ,jj+1) + ht_n (ji+1,jj+1) & 1150 1353 & + ht_n (ji ,jj ) + ht_n (ji+1,jj ) ) & 1151 1354 & / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) & … … 1155 1358 END DO 1156 1359 END SELECT 1157 CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp ) 1158 ! 1159 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1160 DO jj = 2, jpj 1161 DO ji = 2, jpi 1360 ! 1361 DO jj = 2, jpj-1 1362 DO ji = 2, jpi-1 1162 1363 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 1163 1364 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) … … 1166 1367 END DO 1167 1368 END DO 1369 CALL lbc_lnk_multi( 'dynspg_ts', ftne, 'F', 1._wp, ftnw, 'F', 1._wp & 1370 & , ftse, 'F', 1._wp, ftsw, 'F', 1._wp, khlcom = nn_hls+nn_hlts ) 1168 1371 ! 1169 1372 CASE( np_EET ) != EEN scheme using e3t (energy conserving scheme) 1170 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp1171 1373 DO jj = 2, jpj 1172 1374 DO ji = 2, jpi … … 1178 1380 END DO 1179 1381 END DO 1382 CALL lbc_lnk_multi( 'dynspg_ts', ftne, 'F', 1._wp, ftnw, 'F', 1._wp & 1383 & , ftse, 'F', 1._wp, ftsw, 'F', 1._wp, khlcom = nn_hls+nn_hlts ) 1180 1384 ! 1181 1385 CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT ! … … 1223 1427 END DO 1224 1428 END DO 1225 CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp ) 1226 ! JC: TBC. hf should be greater than 0 1227 DO jj = 1, jpj 1228 DO ji = 1, jpi 1429 ! JC: TBC. hf should be greater than 0 1430 DO jj = 2, jpjm1 1431 DO ji = 2, jpim1 1229 1432 IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) 1230 1433 END DO 1231 1434 END DO 1232 zwz(:,:) = ff_f(:,:) * zwz(:,:) 1435 zwz(2:jpim1,2:jpjm1) = ff_f(2:jpim1,2:jpjm1) * zwz(2:jpim1,2:jpjm1) 1436 CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp, khlcom = nn_hls+nn_hlts ) 1233 1437 END SELECT 1234 1438 … … 1237 1441 1238 1442 1239 SUBROUTINE dyn_cor_2d( hu_n, hv_n, un_b, vn_b, zhU, zhV, zu_trd, zv_trd ) 1443 SUBROUTINE dyn_cor_2d( phgtt, phgtu, phgtv, pun, pvn, phU, phV, kdbi , kdei , kdbj , kdej , pu_trd, pv_trd & 1444 & , kdbi2, kdei2, kdbj2, kdej2 ) 1240 1445 !!--------------------------------------------------------------------- 1241 1446 !! *** ROUTINE dyn_cor_2d *** 1242 1447 !! 1243 1448 !! ** Purpose : Compute u and v coriolis trends 1244 !!---------------------------------------------------------------------- 1245 INTEGER :: ji ,jj ! dummy loop indices 1246 REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! - - 1247 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: hu_n, hv_n, un_b, vn_b, zhU, zhV 1248 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd 1449 !! 1450 !! kdXX2 are useful in the initialisation where some arrays are not over the whole domain 1451 !! and some are 1452 !!---------------------------------------------------------------------- 1453 REAL(wp), DIMENSION(kdbi :kdei ,kdbj :kdej ), INTENT(in ) :: phgtt, phgtu, phgtv, pun, pvn ! height, speed 1454 INTEGER , INTENT(in ) :: kdbi , kdei , kdbj , kdej ! arrays size 1455 REAL(wp), DIMENSION(kdbi2:kdei2,kdbj2:kdej2), INTENT(in ) :: phU, phV ! flux 1456 REAL(wp), DIMENSION(kdbi2:kdei2,kdbj2:kdej2), INTENT( out) :: pu_trd, pv_trd 1457 INTEGER , INTENT(in ) :: kdbi2, kdei2, kdbj2, kdej2 ! arrays size 1458 INTEGER :: ji, jj ! dummy loop indices 1459 REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! local integer 1249 1460 !!---------------------------------------------------------------------- 1250 1461 SELECT CASE( nvor_scheme ) 1251 1462 CASE( np_ENT ) ! enstrophy conserving scheme (f-point) 1252 DO jj = 2, jpjm11253 DO ji = 2, jpim11254 z1_hu = ssumask (ji,jj) / ( hu_n(ji,jj) + 1._wp - ssumask(ji,jj) )1255 z1_hv = ssvmask (ji,jj) / ( hv_n(ji,jj) + 1._wp - ssvmask(ji,jj) )1256 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu&1257 & * ( e1e2t (ji+1,jj)*ht_n(ji+1,jj)*ff_t(ji+1,jj) * ( vn_b(ji+1,jj) + vn_b(ji+1,jj-1) ) &1258 & + e1e2t (ji ,jj)*ht_n(ji ,jj)*ff_t(ji ,jj) * ( vn_b(ji ,jj) + vn_b(ji ,jj-1) ) )1463 DO jj = kdbj+1, kdej-1 1464 DO ji = kdbi+1, kdei-1 1465 z1_hu = ssumask_xtd(ji,jj) / ( phgtu(ji,jj) + 1._wp - ssumask_xtd(ji,jj) ) 1466 z1_hv = ssvmask_xtd(ji,jj) / ( phgtv(ji,jj) + 1._wp - ssvmask_xtd(ji,jj) ) 1467 pu_trd(ji,jj) = + r1_4 * r1_e1e2u_xtd(ji,jj) * z1_hu & 1468 & * ( e1e2t_xtd(ji+1,jj)*phgtt(ji+1,jj)*ff_t_xtd(ji+1,jj) * ( pvn(ji+1,jj) + pvn(ji+1,jj-1) ) & 1469 & + e1e2t_xtd(ji ,jj)*phgtt(ji ,jj)*ff_t_xtd(ji ,jj) * ( pvn(ji ,jj) + pvn(ji ,jj-1) ) ) 1259 1470 ! 1260 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv&1261 & * ( e1e2t (ji,jj+1)*ht_n(ji,jj+1)*ff_t(ji,jj+1) * ( un_b(ji,jj+1) + un_b(ji-1,jj+1) ) &1262 & + e1e2t (ji,jj )*ht_n(ji,jj )*ff_t(ji,jj ) * ( un_b(ji,jj ) + un_b(ji-1,jj ) ) )1471 pv_trd(ji,jj) = - r1_4 * r1_e1e2v_xtd(ji,jj) * z1_hv & 1472 & * ( e1e2t_xtd(ji,jj+1)*phgtt(ji,jj+1)*ff_t_xtd(ji,jj+1) * ( pun(ji,jj+1) + pun(ji-1,jj+1) ) & 1473 & + e1e2t_xtd(ji,jj )*phgtt(ji,jj )*ff_t_xtd(ji,jj ) * ( pun(ji,jj ) + pun(ji-1,jj ) ) ) 1263 1474 END DO 1264 END DO 1475 END DO 1265 1476 ! 1266 1477 CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX 1267 DO jj = 2, jpjm11268 DO ji = fs_2, fs_jpim1 ! vector opt.1269 zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj)1270 zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)1271 zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj)1272 zx2 = ( zhU(ji ,jj) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)1478 DO jj = kdbj+1, kdej-1 1479 DO ji = kdbi+1, kdei-1 1480 zy1 = ( phV(ji,jj-1) + phV(ji+1,jj-1) ) * r1_e1u_xtd(ji,jj) 1481 zy2 = ( phV(ji,jj ) + phV(ji+1,jj ) ) * r1_e1u_xtd(ji,jj) 1482 zx1 = ( phU(ji-1,jj) + phU(ji-1,jj+1) ) * r1_e2v_xtd(ji,jj) 1483 zx2 = ( phU(ji ,jj) + phU(ji ,jj+1) ) * r1_e2v_xtd(ji,jj) 1273 1484 ! energy conserving formulation for planetary vorticity term 1274 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )1275 zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )1485 pu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 1486 pv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 1276 1487 END DO 1277 1488 END DO 1278 1489 ! 1279 1490 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 1280 DO jj = 2, jpjm11281 DO ji = fs_2, fs_jpim1 ! vector opt.1282 zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) &1283 & + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj)1284 zx1 = - r1_8 * ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) &1285 & + zhU(ji ,jj ) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj)1286 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) )1287 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) )1288 END DO 1289 END DO 1290 ! 1291 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 1292 DO jj = 2, jpjm11293 DO ji = fs_2, fs_jpim1 ! vector opt.1294 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) &1295 & + ftnw(ji+1,jj) * zhV(ji+1,jj ) &1296 & + ftse(ji,jj ) * zhV(ji ,jj-1) &1297 & + ftsw(ji+1,jj) * zhV(ji+1,jj-1) )1298 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) &1299 & + ftse(ji,jj+1) * zhU(ji ,jj+1) &1300 & + ftnw(ji,jj ) * zhU(ji-1,jj ) &1301 & + ftne(ji,jj ) * zhU(ji ,jj ) )1491 DO jj = kdbj+1, kdej-1 1492 DO ji = kdbi+1, kdei-1 1493 zy1 = r1_8 * ( phV(ji ,jj-1) + phV(ji+1,jj-1) & 1494 & + phV(ji ,jj ) + phV(ji+1,jj ) ) * r1_e1u_xtd(ji,jj) 1495 zx1 = - r1_8 * ( phU(ji-1,jj ) + phU(ji-1,jj+1) & 1496 & + phU(ji ,jj ) + phU(ji ,jj+1) ) * r1_e2v_xtd(ji,jj) 1497 pu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 1498 pv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 1499 END DO 1500 END DO 1501 ! 1502 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 1503 DO jj = kdbj+1, kdej-1 1504 DO ji = kdbi+1, kdei-1 1505 pu_trd(ji,jj) = + r1_12 * r1_e1u_xtd(ji,jj) * ( ftne(ji,jj ) * phV(ji ,jj ) & 1506 & + ftnw(ji+1,jj) * phV(ji+1,jj ) & 1507 & + ftse(ji,jj ) * phV(ji ,jj-1) & 1508 & + ftsw(ji+1,jj) * phV(ji+1,jj-1) ) 1509 pv_trd(ji,jj) = - r1_12 * r1_e2v_xtd(ji,jj) * ( ftsw(ji,jj+1) * phU(ji-1,jj+1) & 1510 & + ftse(ji,jj+1) * phU(ji ,jj+1) & 1511 & + ftnw(ji,jj ) * phU(ji-1,jj ) & 1512 & + ftne(ji,jj ) * phU(ji ,jj ) ) 1302 1513 END DO 1303 1514 END DO … … 1305 1516 END SELECT 1306 1517 ! 1307 END SUBROUTINE dyn_cor_2 D1518 END SUBROUTINE dyn_cor_2d 1308 1519 1309 1520 … … 1448 1659 1449 1660 1450 SUBROUTINE d yn_drg_init( pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )1451 !!---------------------------------------------------------------------- 1452 !! *** ROUTINE d yn_drg_init ***1661 SUBROUTINE drg_init( kdbi, kdei, kdbj, kdej, pu_frc, pv_frc, pCdU_u, pCdU_v ) 1662 !!---------------------------------------------------------------------- 1663 !! *** ROUTINE drg_init *** 1453 1664 !! 1454 1665 !! ** Purpose : - add the baroclinic top/bottom drag contribution to … … 1458 1669 !! ** Method : computation done over the INNER domain only 1459 1670 !!---------------------------------------------------------------------- 1460 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pu_RHSi, pv_RHSi ! baroclinic part of the barotropic RHS 1461 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients 1671 INTEGER , INTENT(in ) :: kdbi, kdei, kdbj, kdej ! arrays size 1672 REAL(wp), DIMENSION(kdbi:kdei,kdbj:kdej), INTENT(inout) :: pu_frc, pv_frc ! baroclinic part of the barotropic RHS 1673 REAL(wp), DIMENSION(kdbi:kdei,kdbj:kdej), INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients 1462 1674 ! 1463 1675 INTEGER :: ji, jj ! dummy loop indices … … 1514 1726 DO jj = 2, jpjm1 1515 1727 DO ji = 2, jpim1 ! INNER domain 1516 pu_ RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( &1728 pu_frc(ji,jj) = pu_frc(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( & 1517 1729 & r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) 1518 pv_ RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( &1730 pv_frc(ji,jj) = pv_frc(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( & 1519 1731 & r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp ) 1520 1732 END DO … … 1524 1736 DO jj = 2, jpjm1 1525 1737 DO ji = 2, jpim1 ! INNER domain 1526 pu_ RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj)1527 pv_ RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj)1738 pu_frc(ji,jj) = pu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj) 1739 pv_frc(ji,jj) = pv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj) 1528 1740 END DO 1529 1741 END DO … … 1560 1772 DO jj = 2, jpjm1 1561 1773 DO ji = 2, jpim1 ! INNER domain 1562 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu_n(ji,jj) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj) 1563 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv_n(ji,jj) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj) 1564 END DO 1565 END DO 1566 ! 1567 ENDIF 1568 ! 1569 END SUBROUTINE dyn_drg_init 1570 1571 SUBROUTINE ts_bck_interp( jn, ll_init, & ! <== in 1774 pu_frc(ji,jj) = pu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj) 1775 pv_frc(ji,jj) = pv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj) 1776 END DO 1777 END DO 1778 ! 1779 ENDIF 1780 ! 1781 END SUBROUTINE drg_init 1782 1783 1784 SUBROUTINE ts_bck_interp( km, ld_init, & ! <== in 1572 1785 & za0, za1, za2, za3 ) ! ==> out 1573 1786 !!---------------------------------------------------------------------- 1574 INTEGER ,INTENT(in ) :: jn! index of sub time step1575 LOGICAL ,INTENT(in ) :: l l_init !1787 INTEGER ,INTENT(in ) :: km ! index of sub time step 1788 LOGICAL ,INTENT(in ) :: ld_init ! 1576 1789 REAL(wp),INTENT( out) :: za0, za1, za2, za3 ! Half-step back interpolation coefficient 1577 1790 ! … … 1579 1792 !!---------------------------------------------------------------------- 1580 1793 ! ! set Half-step back interpolation coefficient 1581 IF ( jn==1 .AND. ll_init ) THEN !* Forward-backward1794 IF ( km==1 .AND. ld_init ) THEN !* Forward-backward 1582 1795 za0 = 1._wp 1583 1796 za1 = 0._wp 1584 1797 za2 = 0._wp 1585 1798 za3 = 0._wp 1586 ELSEIF( jn==2 .AND. ll_init ) THEN !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/121799 ELSEIF( km==2 .AND. ld_init ) THEN !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12 1587 1800 za0 = 1.0833333333333_wp ! za0 = 1-gam-eps 1588 1801 za1 =-0.1666666666666_wp ! za1 = gam
Note: See TracChangeset
for help on using the changeset viewer.
