- Timestamp:
- 2020-12-02T16:32:24+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN
- Files:
-
- 3 edited
-
dynatf.F90 (modified) (16 diffs)
-
dynhpg.F90 (modified) (21 diffs)
-
dynspg_ts.F90 (modified) (70 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynatf.F90
r13981 r14017 13 13 !! - ! 2002-10 (C. Talandier, A-M. Treguier) Open boundary cond. 14 14 !! 2.0 ! 2005-11 (V. Garnier) Surface pressure gradient organization 15 !! 2.3 ! 2007-07 (D. Storkey) Calls to BDY routines. 15 !! 2.3 ! 2007-07 (D. Storkey) Calls to BDY routines. 16 16 !! 3.2 ! 2009-06 (G. Madec, R.Benshila) re-introduce the vvl option 17 17 !! 3.3 ! 2010-09 (D. Storkey, E.O'Dea) Bug fix for BDY module … … 22 22 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) Rename dynnxt.F90 -> dynatf.F90. Now just does time filtering. 23 23 !!------------------------------------------------------------------------- 24 24 25 25 !!---------------------------------------------------------------------------------------------- 26 26 !! dyn_atf : apply Asselin time filtering to "now" velocities and vertical scale factors … … 42 42 USE trdken ! trend manager: kinetic energy 43 43 USE isf_oce , ONLY: ln_isf ! ice shelf 44 USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine 44 USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine 45 45 ! 46 46 USE in_out_manager ! I/O manager … … 81 81 !!---------------------------------------------------------------------- 82 82 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 83 !! $Id$ 83 !! $Id$ 84 84 !! Software governed by the CeCILL license (see ./LICENSE) 85 85 !!---------------------------------------------------------------------- … … 89 89 !!---------------------------------------------------------------------- 90 90 !! *** ROUTINE dyn_atf *** 91 !! 92 !! ** Purpose : Finalize after horizontal velocity. Apply the boundary 91 !! 92 !! ** Purpose : Finalize after horizontal velocity. Apply the boundary 93 93 !! condition on the after velocity and apply the Asselin time 94 94 !! filter to the now fields. … … 97 97 !! estimate (ln_dynspg_ts=T) 98 98 !! 99 !! * Apply lateral boundary conditions on after velocity 99 !! * Apply lateral boundary conditions on after velocity 100 100 !! at the local domain boundaries through lbc_lnk call, 101 101 !! at the one-way open boundaries (ln_bdy=T), … … 104 104 !! * Apply the Asselin time filter to the now fields 105 105 !! arrays to start the next time step: 106 !! (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm)) 106 !! (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm)) 107 107 !! + rn_atfp [ (puu(Kbb),pvv(Kbb)) + (puu(Kaa),pvv(Kaa)) - 2 (puu(Kmm),pvv(Kmm)) ] 108 108 !! Note that with flux form advection and non linear free surface, … … 110 110 !! As a result, dyn_atf MUST be called after tra_atf. 111 111 !! 112 !! ** Action : puu(Kmm),pvv(Kmm) filtered now horizontal velocity 112 !! ** Action : puu(Kmm),pvv(Kmm) filtered now horizontal velocity 113 113 !!---------------------------------------------------------------------- 114 114 INTEGER , INTENT(in ) :: kt ! ocean time-step index … … 122 122 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zue, zve, zwfld 123 123 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zutau, zvtau 124 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze3t_f, ze3u_f, ze3v_f, zua, zva 124 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze3t_f, ze3u_f, ze3v_f, zua, zva 125 125 !!---------------------------------------------------------------------- 126 126 ! … … 150 150 ! 151 151 IF( .NOT.ln_bt_fw ) THEN 152 ! Remove advective velocity from "now velocities" 153 ! prior to asselin filtering 154 ! In the forward case, this is done below after asselin filtering 155 ! so that asselin contribution is removed at the same time 152 ! Remove advective velocity from "now velocities" 153 ! prior to asselin filtering 154 ! In the forward case, this is done below after asselin filtering 155 ! so that asselin contribution is removed at the same time 156 156 DO jk = 1, jpkm1 157 157 puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk) 158 158 pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk) 159 END DO 159 END DO 160 160 ENDIF 161 161 ENDIF 162 162 163 163 ! Update after velocity on domain lateral boundaries 164 ! -------------------------------------------------- 164 ! -------------------------------------------------- 165 165 # if defined key_agrif 166 166 CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries … … 194 194 ! Time filter and swap of dynamics arrays 195 195 ! ------------------------------------------ 196 197 IF( .NOT. l_1st_euler ) THEN !* Leap-Frog : Asselin time filter 196 197 IF( .NOT. l_1st_euler ) THEN !* Leap-Frog : Asselin time filter 198 198 ! ! =============! 199 199 IF( ln_linssh ) THEN ! Fixed volume ! … … 220 220 DO jk = 1, jpkm1 221 221 ze3t_f(:,:,jk) = ze3t_f(:,:,jk) - zcoef * zwfld(:,:) * tmask(:,:,jk) & 222 & * pe3t(:,:,jk,Kmm) / ( ht(:,:) + 1._wp - ssmask(:,:) ) 222 & * pe3t(:,:,jk,Kmm) / ( ht(:,:) + 1._wp - ssmask(:,:) ) 223 223 END DO 224 224 ! … … 257 257 pvv(ji,jj,jk,Kmm) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk) 258 258 END_3D 259 pe3u(:,:,1:jpkm1,Kmm) = ze3u_f(:,:,1:jpkm1) 259 pe3u(:,:,1:jpkm1,Kmm) = ze3u_f(:,:,1:jpkm1) 260 260 pe3v(:,:,1:jpkm1,Kmm) = ze3v_f(:,:,1:jpkm1) 261 261 ! … … 268 268 IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN 269 269 ! Revert filtered "now" velocities to time split estimate 270 ! Doing it here also means that asselin filter contribution is removed 270 ! Doing it here also means that asselin filter contribution is removed 271 271 zue(:,:) = pe3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1) 272 zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1) 272 zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1) 273 273 DO jk = 2, jpkm1 274 274 zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk) 275 zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk) 275 zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk) 276 276 END DO 277 277 DO jk = 1, jpkm1 … … 325 325 IF ( iom_use("utau") ) THEN 326 326 IF ( ln_drgice_imp.OR.ln_isfcav ) THEN 327 ALLOCATE(zutau(jpi,jpj)) 327 ALLOCATE(zutau(jpi,jpj)) 328 328 DO_2D( 0, 0, 0, 0 ) 329 jk = miku(ji,jj) 329 jk = miku(ji,jj) 330 330 zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * puu(ji,jj,jk,Kaa) 331 331 END_2D … … 353 353 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' nxt - puu(:,:,:,Kaa): ', mask1=umask, & 354 354 & tab3d_2=pvv(:,:,:,Kaa), clinfo2=' pvv(:,:,:,Kaa): ' , mask2=vmask ) 355 ! 355 ! 356 356 IF( ln_dynspg_ts ) DEALLOCATE( zue, zve ) 357 357 IF( l_trddyn ) DEALLOCATE( zua, zva ) -
NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynhpg.F90
r13295 r14017 43 43 USE in_out_manager ! I/O manager 44 44 USE prtctl ! Print control 45 USE lbclnk ! lateral boundary condition 45 USE lbclnk ! lateral boundary condition 46 46 USE lib_mpp ! MPP library 47 47 USE eosbn2 ! compute density … … 206 206 ! 207 207 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 208 ! 208 ! 209 209 IF(lwp) THEN 210 210 WRITE(numout,*) … … 219 219 WRITE(numout,*) 220 220 ENDIF 221 ! 221 ! 222 222 END SUBROUTINE dyn_hpg_init 223 223 … … 302 302 INTEGER :: iku, ikv ! temporary integers 303 303 REAL(wp) :: zcoef0, zcoef1, zcoef2, zcoef3 ! temporary scalars 304 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj 305 REAL(wp), DIMENSION(jpi,jpj) :: zgtsu, zgtsv, zgru, zgrv 304 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj 305 REAL(wp), DIMENSION(jpi,jpj,jpts) :: zgtsu, zgtsv 306 REAL(wp), DIMENSION(jpi,jpj) :: zgru, zgrv 306 307 !!---------------------------------------------------------------------- 307 308 ! … … 429 430 zcpx(ji,jj) = 0._wp 430 431 END IF 431 432 432 433 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 433 434 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & … … 470 471 IF( ln_wd_il ) THEN 471 472 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 472 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 473 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 473 474 zuap = zuap * zcpx(ji,jj) 474 475 zvap = zvap * zcpy(ji,jj) … … 497 498 IF( ln_wd_il ) THEN 498 499 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 499 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 500 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 500 501 zuap = zuap * zcpx(ji,jj) 501 502 zvap = zvap * zcpy(ji,jj) … … 528 529 !! pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj 529 530 !! iceload is added 530 !! 531 !! 531 532 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 532 533 !!---------------------------------------------------------------------- … … 546 547 znad=1._wp ! To use density and not density anomaly 547 548 ! 548 ! ! iniitialised to 0. zhpi zhpi 549 ! ! iniitialised to 0. zhpi zhpi 549 550 zhpi(:,:,:) = 0._wp ; zhpj(:,:,:) = 0._wp 550 551 … … 560 561 CALL eos( zts_top, risfdep, zrhdtop_oce ) 561 562 562 !================================================================================== 563 !===== Compute surface value ===================================================== 563 !================================================================================== 564 !===== Compute surface value ===================================================== 564 565 !================================================================================== 565 566 DO_2D( 0, 0, 0, 0 ) … … 573 574 & - 0.5_wp * e3w(ji,jj,ikt,Kmm) & 574 575 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 575 & + ( risfload(ji+1,jj) - risfload(ji,jj)) ) 576 & + ( risfload(ji+1,jj) - risfload(ji,jj)) ) 576 577 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w(ji,jj+1,iktp1j,Kmm) & 577 578 & * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) & 578 & - 0.5_wp * e3w(ji,jj,ikt,Kmm) & 579 & - 0.5_wp * e3w(ji,jj,ikt,Kmm) & 579 580 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 580 & + ( risfload(ji,jj+1) - risfload(ji,jj)) ) 581 & + ( risfload(ji,jj+1) - risfload(ji,jj)) ) 581 582 ! s-coordinate pressure gradient correction (=0 if z coordinate) 582 583 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & … … 588 589 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1) 589 590 END_2D 590 !================================================================================== 591 !===== Compute interior value ===================================================== 591 !================================================================================== 592 !===== Compute interior value ===================================================== 592 593 !================================================================================== 593 594 ! interior value (2=<jk=<jpkm1) … … 660 661 zcpx(ji,jj) = 0._wp 661 662 END IF 662 663 663 664 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 664 665 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & … … 835 836 IF( ln_wd_il ) THEN 836 837 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 837 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 838 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 838 839 ENDIF 839 840 ! add to the general momentum trend … … 855 856 IF( ln_wd_il ) THEN 856 857 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 857 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 858 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 858 859 ENDIF 859 860 ! add to the general momentum trend … … 926 927 zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 927 928 & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) ) 928 929 929 930 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 930 931 ELSE 931 932 zcpx(ji,jj) = 0._wp 932 933 END IF 933 934 934 935 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 935 936 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & … … 1012 1013 !!gm BUG ? if it is ssh at u- & v-point then it should be: 1013 1014 ! zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji+1,jj) * ssh(ji+1,jj,Kmm)) * & 1014 ! & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1015 ! & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1015 1016 ! zsshv_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji,jj+1) * ssh(ji,jj+1,Kmm)) * & 1016 ! & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1017 ! & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1017 1018 !!gm not this: 1018 1019 zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * & 1019 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1020 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1020 1021 zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * & 1021 & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1022 & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1022 1023 END_2D 1023 1024 … … 1025 1026 1026 1027 DO_2D( 0, 0, 0, 0 ) 1027 zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad) 1028 zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad) 1028 1029 zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) * znad) 1029 1030 END_2D … … 1108 1109 zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj) 1109 1110 ENDIF 1110 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 1111 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 1111 1112 ENDIF 1112 1113 … … 1164 1165 ENDIF 1165 1166 IF( ln_wd_il ) THEN 1166 zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj) 1167 zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj) 1167 zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj) 1168 zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj) 1168 1169 ENDIF 1169 1170 … … 1199 1200 !!---------------------------------------------------------------------- 1200 1201 ! 1201 !!gm WHAT !!!!! THIS IS VERY DANGEROUS !!!!! 1202 !!gm WHAT !!!!! THIS IS VERY DANGEROUS !!!!! 1202 1203 jpi = size(fsp,1) 1203 1204 jpj = size(fsp,2) -
NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/DYN/dynspg_ts.F90
r13546 r14017 1 1 MODULE dynspg_ts 2 2 3 !! Includes ROMS wd scheme with diagnostic outputs ; puu(:,:,:,Kmm) and puu(:,:,:,Krhs) updates are commented out ! 3 !! Includes ROMS wd scheme with diagnostic outputs ; puu(:,:,:,Kmm) and puu(:,:,:,Krhs) updates are commented out ! 4 4 5 5 !!====================================================================== … … 23 23 24 24 !!---------------------------------------------------------------------- 25 !! dyn_spg_ts : compute surface pressure gradient trend using a time-splitting scheme 25 !! dyn_spg_ts : compute surface pressure gradient trend using a time-splitting scheme 26 26 !! dyn_spg_ts_init: initialisation of the time-splitting scheme 27 27 !! ts_wgt : set time-splitting weights for temporal averaging (or not) … … 50 50 USE agrif_oce 51 51 #endif 52 #if defined key_asminc 52 #if defined key_asminc 53 53 USE asminc ! Assimilation increment 54 54 #endif … … 66 66 PRIVATE 67 67 68 PUBLIC dyn_spg_ts ! called by dyn_spg 68 PUBLIC dyn_spg_ts ! called by dyn_spg 69 69 PUBLIC dyn_spg_ts_init ! - - dyn_spg_init 70 70 … … 122 122 !! ** Purpose : - Compute the now trend due to the explicit time stepping 123 123 !! of the quasi-linear barotropic system, and add it to the 124 !! general momentum trend. 124 !! general momentum trend. 125 125 !! 126 126 !! ** Method : - split-explicit schem (time splitting) : 127 127 !! Barotropic variables are advanced from internal time steps 128 128 !! "n" to "n+1" if ln_bt_fw=T 129 !! or from 129 !! or from 130 130 !! "n-1" to "n+1" if ln_bt_fw=F 131 131 !! thanks to a generalized forward-backward time stepping (see ref. below). … … 136 136 !! -Compute barotropic advective fluxes at step "n" : un_adv, vn_adv 137 137 !! These are used to advect tracers and are compliant with discrete 138 !! continuity equation taken at the baroclinic time steps. This 138 !! continuity equation taken at the baroclinic time steps. This 139 139 !! ensures tracers conservation. 140 140 !! - (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) momentum trend updated with barotropic component. 141 141 !! 142 !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005. 142 !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005. 143 143 !!--------------------------------------------------------------------- 144 144 INTEGER , INTENT( in ) :: kt ! ocean time-step index … … 148 148 ! 149 149 INTEGER :: ji, jj, jk, jn ! dummy loop indices 150 LOGICAL :: ll_fw_start ! =T : forward integration 150 LOGICAL :: ll_fw_start ! =T : forward integration 151 151 LOGICAL :: ll_init ! =T : special startup of 2d equations 152 152 INTEGER :: noffset ! local integers : time offset for bdy update … … 164 164 REAL(wp), DIMENSION(jpi, jpj, jpk) :: ze3u, ze3v 165 165 ! 166 REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True. 166 REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True. 167 167 168 168 INTEGER :: iwdg, jwdg, kwdg ! short-hand values for the indices of the output point … … 179 179 IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj)) 180 180 ! 181 zwdramp = r_rn_wdmin1 ! simplest ramp 181 zwdramp = r_rn_wdmin1 ! simplest ramp 182 182 ! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp 183 ! ! inverse of baroclinic time step 184 r1_Dt_b = 1._wp / rDt 185 ! 186 ll_init = ln_bt_av ! if no time averaging, then no specific restart 183 ! ! inverse of baroclinic time step 184 r1_Dt_b = 1._wp / rDt 185 ! 186 ll_init = ln_bt_av ! if no time averaging, then no specific restart 187 187 ll_fw_start = .FALSE. 188 188 ! ! time offset in steps for bdy data update 189 189 IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_e 190 ELSE ; noffset = 0 190 ELSE ; noffset = 0 191 191 ENDIF 192 192 ! … … 215 215 ! and the averaging weights. We don't have an easy way of telling whether we did 216 216 ! an Euler timestep on the first timestep (because l_1st_euler is reset to .false. 217 ! at the end of the first timestep) so just do this in all cases. 217 ! at the end of the first timestep) so just do this in all cases. 218 218 ll_fw_start = .FALSE. 219 219 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) … … 225 225 ! Phase 1 : Coupling between general trend and barotropic estimates (1st step) 226 226 ! ----------------------------------------------------------------------------- 227 ! 227 ! 228 228 ! 229 229 ! != zu_frc = 1/H e3*d/dt(Ua) =! (Vertical mean of Ua, the 3D trends) … … 243 243 vv(:,:,jk,Krhs) = ( vv(:,:,jk,Krhs) - zv_frc(:,:) ) * vmask(:,:,jk) 244 244 END DO 245 245 246 246 !!gm Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum.... 247 247 !!gm Is it correct to do so ? I think so... 248 248 249 249 ! != remove 2D Coriolis and pressure gradient trends =! 250 250 ! ! ------------------------------------------------- ! … … 254 254 ! 255 255 ! !* 2D Coriolis trends 256 zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:) ! now fluxes 256 zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:) ! now fluxes 257 257 zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:) ! NB: FULL domain : put a value in last row and column 258 258 ! … … 273 273 DO_2D( 0, 0, 0, 0 ) 274 274 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj ,Kmm) - pssh(ji ,jj ,Kmm) ) * r1_e1u(ji,jj) 275 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( pssh(ji ,jj+1,Kmm) - pssh(ji ,jj ,Kmm) ) * r1_e2v(ji,jj) 275 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( pssh(ji ,jj+1,Kmm) - pssh(ji ,jj ,Kmm) ) * r1_e2v(ji,jj) 276 276 END_2D 277 277 ENDIF … … 319 319 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kmm) 320 320 END_2D 321 ENDIF 321 ENDIF 322 322 ! 323 323 ! !----------------! … … 369 369 ! 370 370 ! ----------------------------------------------------------------------- 371 ! Phase 2 : Integration of the barotropic equations 371 ! Phase 2 : Integration of the barotropic equations 372 372 ! ----------------------------------------------------------------------- 373 373 ! 374 374 ! ! ==================== ! 375 375 ! ! Initialisations ! 376 ! ! ==================== ! 377 ! Initialize barotropic variables: 376 ! ! ==================== ! 377 ! Initialize barotropic variables: 378 378 IF( ll_init )THEN 379 379 sshbb_e(:,:) = 0._wp … … 391 391 ENDIF 392 392 ! 393 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 394 sshn_e(:,:) = pssh(:,:,Kmm) 395 un_e (:,:) = puu_b(:,:,Kmm) 393 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 394 sshn_e(:,:) = pssh(:,:,Kmm) 395 un_e (:,:) = puu_b(:,:,Kmm) 396 396 vn_e (:,:) = pvv_b(:,:,Kmm) 397 397 ! 398 hu_e (:,:) = hu(:,:,Kmm) 399 hv_e (:,:) = hv(:,:,Kmm) 400 hur_e (:,:) = r1_hu(:,:,Kmm) 398 hu_e (:,:) = hu(:,:,Kmm) 399 hv_e (:,:) = hv(:,:,Kmm) 400 hur_e (:,:) = r1_hu(:,:,Kmm) 401 401 hvr_e (:,:) = r1_hv(:,:,Kmm) 402 402 ELSE ! CENTRED integration: start from BEFORE fields 403 403 sshn_e(:,:) = pssh(:,:,Kbb) 404 un_e (:,:) = puu_b(:,:,Kbb) 404 un_e (:,:) = puu_b(:,:,Kbb) 405 405 vn_e (:,:) = pvv_b(:,:,Kbb) 406 406 ! 407 hu_e (:,:) = hu(:,:,Kbb) 408 hv_e (:,:) = hv(:,:,Kbb) 409 hur_e (:,:) = r1_hu(:,:,Kbb) 407 hu_e (:,:) = hu(:,:,Kbb) 408 hv_e (:,:) = hv(:,:,Kbb) 409 hur_e (:,:) = r1_hu(:,:,Kbb) 410 410 hvr_e (:,:) = r1_hv(:,:,Kbb) 411 411 ENDIF 412 412 ! 413 413 ! Initialize sums: 414 puu_b (:,:,Kaa) = 0._wp ! After barotropic velocities (or transport if flux form) 414 puu_b (:,:,Kaa) = 0._wp ! After barotropic velocities (or transport if flux form) 415 415 pvv_b (:,:,Kaa) = 0._wp 416 416 pssh (:,:,Kaa) = 0._wp ! Sum for after averaged sea level … … 419 419 ! 420 420 IF( ln_wd_dl ) THEN 421 zuwdmask(:,:) = 0._wp ! set to zero for definiteness (not sure this is necessary) 422 zvwdmask(:,:) = 0._wp ! 423 zuwdav2 (:,:) = 0._wp 424 zvwdav2 (:,:) = 0._wp 425 END IF 421 zuwdmask(:,:) = 0._wp ! set to zero for definiteness (not sure this is necessary) 422 zvwdmask(:,:) = 0._wp ! 423 zuwdav2 (:,:) = 0._wp 424 zvwdav2 (:,:) = 0._wp 425 END IF 426 426 427 427 ! ! ==================== ! … … 443 443 ! 444 444 ! !* Set extrapolation coefficients for predictor step: 445 IF ((jn<3).AND.ll_init) THEN ! Forward 446 za1 = 1._wp 447 za2 = 0._wp 448 za3 = 0._wp 449 ELSE ! AB3-AM4 Coefficients: bet=0.281105 445 IF ((jn<3).AND.ll_init) THEN ! Forward 446 za1 = 1._wp 447 za2 = 0._wp 448 za3 = 0._wp 449 ELSE ! AB3-AM4 Coefficients: bet=0.281105 450 450 za1 = 1.781105_wp ! za1 = 3/2 + bet 451 451 za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet) … … 467 467 !--------------------------------------------------------------------------------! 468 468 zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 469 469 470 470 ! set wetting & drying mask at tracer points for this barotropic mid-step 471 471 IF( ln_wd_dl ) CALL wad_tmsk( zsshp2_e, ztwdmask ) … … 493 493 ! ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d 494 494 IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e ) 495 ! 495 ! 496 496 ! ! resulting flux at mid-step (not over the full domain) 497 497 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 … … 504 504 IF( ln_wd_il ) CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rDt_e) !!gm wad_lmt_bt use of lbc_lnk on zhU, zhV 505 505 506 IF( ln_wd_dl ) THEN ! un_e and vn_e are set to zero at faces where 506 IF( ln_wd_dl ) THEN ! un_e and vn_e are set to zero at faces where 507 507 ! ! the direction of the flow is from dry cells 508 508 CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask ) ! not jpi colomn for U, not jpj row for V 509 509 ! 510 ENDIF 510 ENDIF 511 511 ! 512 512 ! … … 532 532 un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:) 533 533 vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:) 534 ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True) 534 ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True) 535 535 IF ( ln_wd_dl_bc ) THEN 536 536 zuwdav2(1:jpim1,1:jpj ) = zuwdav2(1:jpim1,1:jpj ) + za2 * zuwdmask(1:jpim1,1:jpj ) ! not jpi-column … … 538 538 END IF 539 539 ! 540 ! 540 ! 541 541 ! Sea Surface Height at u-,v-points (vvl case only) 542 IF( .NOT.ln_linssh ) THEN 542 IF( .NOT.ln_linssh ) THEN 543 543 DO_2D( 0, 0, 0, 0 ) 544 544 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & … … 549 549 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) 550 550 END_2D 551 ENDIF 552 ! 551 ENDIF 552 ! 553 553 ! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2 554 554 !-- m+1/2 m+1 m m-1 m-2 --! … … 607 607 IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form 608 608 DO_2D( 0, 0, 0, 0 ) 609 ua_e(ji,jj) = ( un_e(ji,jj) & 609 ua_e(ji,jj) = ( un_e(ji,jj) & 610 610 & + rDt_e * ( zu_spg(ji,jj) & 611 611 & + zu_trd(ji,jj) & 612 & + zu_frc(ji,jj) ) & 612 & + zu_frc(ji,jj) ) & 613 613 & ) * ssumask(ji,jj) 614 614 … … 632 632 z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) ) 633 633 ! 634 ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) & 634 ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) & 635 635 & + rDt_e * ( zhu_bck * zu_spg (ji,jj) & ! 636 636 & + zhup2_e(ji,jj) * zu_trd (ji,jj) & ! … … 650 650 END_2D 651 651 ENDIF 652 652 653 653 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 654 654 hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1) … … 667 667 ! ! open boundaries 668 668 IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 669 #if defined key_agrif 669 #if defined key_agrif 670 670 IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif 671 671 #endif … … 686 686 ! !* Sum over whole bt loop 687 687 ! ! ---------------------- 688 za1 = wgtbtp1(jn) 688 za1 = wgtbtp1(jn) 689 689 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities 690 puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) 691 pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) 690 puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) 691 pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) 692 692 ELSE ! Sum transports 693 IF ( .NOT.ln_wd_dl ) THEN 693 IF ( .NOT.ln_wd_dl ) THEN 694 694 puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:) 695 695 pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:) 696 ELSE 696 ELSE 697 697 puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:) 698 698 pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:) 699 END IF 699 END IF 700 700 ENDIF 701 701 ! ! Sum sea level … … 715 715 zun_save = un_adv(ji,jj) 716 716 zvn_save = vn_adv(ji,jj) 717 ! ! apply the previously computed correction 717 ! ! apply the previously computed correction 718 718 un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - rn_atfp * un_bf(ji,jj) ) 719 719 vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - rn_atfp * vn_bf(ji,jj) ) … … 763 763 764 764 765 ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases) 765 ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases) 766 766 DO jk = 1, jpkm1 767 767 puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) + un_adv(:,:)*r1_hu(:,:,Kmm) - puu_b(:,:,Kmm) ) * umask(:,:,jk) … … 769 769 END DO 770 770 771 IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN 771 IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN 772 772 DO jk = 1, jpkm1 773 773 puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu(:,:,Kmm) & 774 & + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk) 775 pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) & 776 & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk) 774 & + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk) 775 pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) & 776 & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk) 777 777 END DO 778 END IF 779 780 778 END IF 779 780 781 781 CALL iom_put( "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) ) ! barotropic i-current 782 782 CALL iom_put( "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) ) ! barotropic i-current … … 796 796 vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:) 797 797 ENDIF 798 #endif 798 #endif 799 799 ! !* write time-spliting arrays in the restart 800 800 IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' ) … … 817 817 LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true. 818 818 LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true. 819 INTEGER, INTENT(inout) :: jpit ! cycle length 819 INTEGER, INTENT(inout) :: jpit ! cycle length 820 820 REAL(wp), DIMENSION(3*nn_e), INTENT(inout) :: zwgt1, & ! Primary weights 821 821 zwgt2 ! Secondary weights 822 822 823 823 INTEGER :: jic, jn, ji ! temporary integers 824 824 REAL(wp) :: za1, za2 … … 829 829 830 830 ! Set time index when averaged value is requested 831 IF (ll_fw) THEN 831 IF (ll_fw) THEN 832 832 jic = nn_e 833 833 ELSE … … 837 837 ! Set primary weights: 838 838 IF (ll_av) THEN 839 ! Define simple boxcar window for primary weights 840 ! (width = nn_e, centered around jic) 839 ! Define simple boxcar window for primary weights 840 ! (width = nn_e, centered around jic) 841 841 SELECT CASE ( nn_bt_flt ) 842 842 CASE( 0 ) ! No averaging … … 846 846 CASE( 1 ) ! Boxcar, width = nn_e 847 847 DO jn = 1, 3*nn_e 848 za1 = ABS(float(jn-jic))/float(nn_e) 848 za1 = ABS(float(jn-jic))/float(nn_e) 849 849 IF (za1 < 0.5_wp) THEN 850 850 zwgt1(jn) = 1._wp … … 855 855 CASE( 2 ) ! Boxcar, width = 2 * nn_e 856 856 DO jn = 1, 3*nn_e 857 za1 = ABS(float(jn-jic))/float(nn_e) 857 za1 = ABS(float(jn-jic))/float(nn_e) 858 858 IF (za1 < 1._wp) THEN 859 859 zwgt1(jn) = 1._wp … … 868 868 jpit = jic 869 869 ENDIF 870 870 871 871 ! Set secondary weights 872 872 DO jn = 1, jpit … … 897 897 !!---------------------------------------------------------------------- 898 898 ! 899 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 899 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 900 900 ! ! --------------- 901 901 IF( ln_rstart .AND. ln_bt_fw .AND. (.NOT.l_1st_euler) ) THEN !* Read the restart file 902 CALL iom_get( numror, jpdom_auto, 'ub2_b' , ub2_b (:,:), cd_type = 'U', psgn = -1._wp , ldxios = lrxios )903 CALL iom_get( numror, jpdom_auto, 'vb2_b' , vb2_b (:,:), cd_type = 'V', psgn = -1._wp , ldxios = lrxios )904 CALL iom_get( numror, jpdom_auto, 'un_bf' , un_bf (:,:), cd_type = 'U', psgn = -1._wp , ldxios = lrxios )905 CALL iom_get( numror, jpdom_auto, 'vn_bf' , vn_bf (:,:), cd_type = 'V', psgn = -1._wp , ldxios = lrxios )902 CALL iom_get( numror, jpdom_auto, 'ub2_b' , ub2_b (:,:), cd_type = 'U', psgn = -1._wp ) 903 CALL iom_get( numror, jpdom_auto, 'vb2_b' , vb2_b (:,:), cd_type = 'V', psgn = -1._wp ) 904 CALL iom_get( numror, jpdom_auto, 'un_bf' , un_bf (:,:), cd_type = 'U', psgn = -1._wp ) 905 CALL iom_get( numror, jpdom_auto, 'vn_bf' , vn_bf (:,:), cd_type = 'V', psgn = -1._wp ) 906 906 IF( .NOT.ln_bt_av ) THEN 907 CALL iom_get( numror, jpdom_auto, 'sshbb_e' , sshbb_e(:,:), cd_type = 'T', psgn = 1._wp , ldxios = lrxios )908 CALL iom_get( numror, jpdom_auto, 'ubb_e' , ubb_e(:,:), cd_type = 'U', psgn = -1._wp , ldxios = lrxios )909 CALL iom_get( numror, jpdom_auto, 'vbb_e' , vbb_e(:,:), cd_type = 'V', psgn = -1._wp , ldxios = lrxios)910 CALL iom_get( numror, jpdom_auto, 'sshb_e' , sshb_e(:,:), cd_type = 'T', psgn = 1._wp , ldxios = lrxios )911 CALL iom_get( numror, jpdom_auto, 'ub_e' , ub_e(:,:), cd_type = 'U', psgn = -1._wp , ldxios = lrxios )912 CALL iom_get( numror, jpdom_auto, 'vb_e' , vb_e(:,:), cd_type = 'V', psgn = -1._wp , ldxios = lrxios)907 CALL iom_get( numror, jpdom_auto, 'sshbb_e' , sshbb_e(:,:), cd_type = 'T', psgn = 1._wp ) 908 CALL iom_get( numror, jpdom_auto, 'ubb_e' , ubb_e(:,:), cd_type = 'U', psgn = -1._wp ) 909 CALL iom_get( numror, jpdom_auto, 'vbb_e' , vbb_e(:,:), cd_type = 'V', psgn = -1._wp ) 910 CALL iom_get( numror, jpdom_auto, 'sshb_e' , sshb_e(:,:), cd_type = 'T', psgn = 1._wp ) 911 CALL iom_get( numror, jpdom_auto, 'ub_e' , ub_e(:,:), cd_type = 'U', psgn = -1._wp ) 912 CALL iom_get( numror, jpdom_auto, 'vb_e' , vb_e(:,:), cd_type = 'V', psgn = -1._wp ) 913 913 ENDIF 914 914 #if defined key_agrif 915 915 ! Read time integrated fluxes 916 916 IF ( .NOT.Agrif_Root() ) THEN 917 CALL iom_get( numror, jpdom_auto, 'ub2_i_b' , ub2_i_b(:,:), cd_type = 'U', psgn = -1._wp , ldxios = lrxios )918 CALL iom_get( numror, jpdom_auto, 'vb2_i_b' , vb2_i_b(:,:), cd_type = 'V', psgn = -1._wp , ldxios = lrxios)917 CALL iom_get( numror, jpdom_auto, 'ub2_i_b' , ub2_i_b(:,:), cd_type = 'U', psgn = -1._wp ) 918 CALL iom_get( numror, jpdom_auto, 'vb2_i_b' , vb2_i_b(:,:), cd_type = 'V', psgn = -1._wp ) 919 919 ELSE 920 920 ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif … … 935 935 ! ! ------------------- 936 936 IF(lwp) WRITE(numout,*) '---- ts_rst ----' 937 IF( lwxios ) CALL iom_swap( cwxios_context ) 938 CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:), ldxios = lwxios ) 939 CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:), ldxios = lwxios ) 940 CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:), ldxios = lwxios ) 941 CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:), ldxios = lwxios ) 937 CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:) ) 938 CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:) ) 939 CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:) ) 940 CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:) ) 942 941 ! 943 942 IF (.NOT.ln_bt_av) THEN 944 CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) , ldxios = lwxios )945 CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:) , ldxios = lwxios)946 CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:) , ldxios = lwxios)947 CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:) , ldxios = lwxios)948 CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:) , ldxios = lwxios)949 CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:) , ldxios = lwxios)943 CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) ) 944 CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:) ) 945 CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:) ) 946 CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:) ) 947 CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:) ) 948 CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:) ) 950 949 ENDIF 951 950 #if defined key_agrif 952 951 ! Save time integrated fluxes 953 952 IF ( .NOT.Agrif_Root() ) THEN 954 CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:) , ldxios = lwxios)955 CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:) , ldxios = lwxios)953 CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:) ) 954 CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:) ) 956 955 ENDIF 957 956 #endif 958 IF( lwxios ) CALL iom_swap( cxios_context )959 957 ENDIF 960 958 ! … … 986 984 ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 987 985 IF( ln_bt_auto ) nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax) 988 986 989 987 rDt_e = rn_Dt / REAL( nn_e , wp ) 990 988 zcmax = zcmax * rDt_e … … 1022 1020 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac' 1023 1021 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_e' 1024 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_e' 1022 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_e' 1025 1023 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' ) 1026 1024 END SELECT … … 1040 1038 ENDIF 1041 1039 IF( zcmax>0.9_wp ) THEN 1042 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' ) 1040 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' ) 1043 1041 ENDIF 1044 1042 ! … … 1049 1047 CALL ts_rst( nit000, 'READ' ) 1050 1048 ! 1051 IF( lwxios ) THEN1052 ! define variables in restart file when writing with XIOS1053 CALL iom_set_rstw_var_active('ub2_b')1054 CALL iom_set_rstw_var_active('vb2_b')1055 CALL iom_set_rstw_var_active('un_bf')1056 CALL iom_set_rstw_var_active('vn_bf')1057 !1058 IF (.NOT.ln_bt_av) THEN1059 CALL iom_set_rstw_var_active('sshbb_e')1060 CALL iom_set_rstw_var_active('ubb_e')1061 CALL iom_set_rstw_var_active('vbb_e')1062 CALL iom_set_rstw_var_active('sshb_e')1063 CALL iom_set_rstw_var_active('ub_e')1064 CALL iom_set_rstw_var_active('vb_e')1065 ENDIF1066 #if defined key_agrif1067 ! Save time integrated fluxes1068 IF ( .NOT.Agrif_Root() ) THEN1069 CALL iom_set_rstw_var_active('ub2_i_b')1070 CALL iom_set_rstw_var_active('vb2_i_b')1071 ENDIF1072 #endif1073 ENDIF1074 !1075 1049 END SUBROUTINE dyn_spg_ts_init 1076 1050 1077 1051 1078 1052 SUBROUTINE dyn_cor_2D_init( Kmm ) 1079 1053 !!--------------------------------------------------------------------- … … 1102 1076 DO_2D( 1, 0, 1, 0 ) 1103 1077 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 1104 & ht(ji ,jj ) + ht(ji+1,jj ) ) * 0.25_wp 1078 & ht(ji ,jj ) + ht(ji+1,jj ) ) * 0.25_wp 1105 1079 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1106 1080 END_2D … … 1138 1112 zwz(:,:) = 0._wp 1139 1113 zhf(:,:) = 0._wp 1140 1141 !!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed 1114 1115 !!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed 1142 1116 !!gm A priori a better value should be something like : 1143 !!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1) 1144 !!gm divided by the sum of the corresponding mask 1145 !!gm 1146 !! 1117 !!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1) 1118 !!gm divided by the sum of the corresponding mask 1119 !!gm 1120 !! 1147 1121 IF( .NOT.ln_sco ) THEN 1148 1122 1149 1123 !!gm agree the JC comment : this should be done in a much clear way 1150 1124 1151 1125 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 1152 ! Set it to zero for the time being 1126 ! Set it to zero for the time being 1153 1127 ! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level 1154 1128 ! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth … … 1177 1151 END DO 1178 1152 CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp ) 1179 ! JC: TBC. hf should be greater than 0 1153 ! JC: TBC. hf should be greater than 0 1180 1154 DO_2D( 1, 1, 1, 1 ) 1181 1155 IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) … … 1183 1157 zwz(:,:) = ff_f(:,:) * zwz(:,:) 1184 1158 END SELECT 1185 1159 1186 1160 END SUBROUTINE dyn_cor_2d_init 1187 1161 … … 1209 1183 ! 1210 1184 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv & 1211 & * ( e1e2t(ji,jj+1)*pht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) ) & 1212 & + e1e2t(ji,jj )*pht(ji,jj )*ff_t(ji,jj ) * ( punb(ji,jj ) + punb(ji-1,jj ) ) ) 1213 END_2D 1214 ! 1185 & * ( e1e2t(ji,jj+1)*pht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) ) & 1186 & + e1e2t(ji,jj )*pht(ji,jj )*ff_t(ji,jj ) * ( punb(ji,jj ) + punb(ji-1,jj ) ) ) 1187 END_2D 1188 ! 1215 1189 CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX 1216 1190 DO_2D( 0, 0, 0, 0 ) … … 1234 1208 END_2D 1235 1209 ! 1236 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 1210 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 1237 1211 DO_2D( 0, 0, 0, 0 ) 1238 1212 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) & … … 1254 1228 !!---------------------------------------------------------------------- 1255 1229 !! *** ROUTINE wad_lmt *** 1256 !! 1257 !! ** Purpose : set wetting & drying mask at tracer points 1258 !! for the current barotropic sub-step 1259 !! 1260 !! ** Method : ??? 1230 !! 1231 !! ** Purpose : set wetting & drying mask at tracer points 1232 !! for the current barotropic sub-step 1233 !! 1234 !! ** Method : ??? 1261 1235 !! 1262 1236 !! ** Action : ptmsk : wetting & drying t-mask … … 1268 1242 !!---------------------------------------------------------------------- 1269 1243 ! 1270 IF( ln_wd_dl_rmp ) THEN 1244 IF( ln_wd_dl_rmp ) THEN 1271 1245 DO_2D( 1, 1, 1, 1 ) 1272 IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 1273 ! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 1246 IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 1247 ! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 1274 1248 ptmsk(ji,jj) = 1._wp 1275 1249 ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 1276 1250 ptmsk(ji,jj) = TANH( 50._wp*( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )*r_rn_wdmin1) ) 1277 ELSE 1251 ELSE 1278 1252 ptmsk(ji,jj) = 0._wp 1279 1253 ENDIF 1280 1254 END_2D 1281 ELSE 1255 ELSE 1282 1256 DO_2D( 1, 1, 1, 1 ) 1283 1257 IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp … … 1293 1267 !!---------------------------------------------------------------------- 1294 1268 !! *** ROUTINE wad_lmt *** 1295 !! 1296 !! ** Purpose : set wetting & drying mask at tracer points 1297 !! for the current barotropic sub-step 1298 !! 1299 !! ** Method : ??? 1269 !! 1270 !! ** Purpose : set wetting & drying mask at tracer points 1271 !! for the current barotropic sub-step 1272 !! 1273 !! ** Method : ??? 1300 1274 !! 1301 1275 !! ** Action : ptmsk : wetting & drying t-mask … … 1309 1283 ! 1310 1284 DO_2D( 1, 1, 1, 0 ) ! not jpi-column 1311 IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj) 1312 ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj) 1285 IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj) 1286 ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj) 1313 1287 ENDIF 1314 1288 phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj) … … 1318 1292 DO_2D( 1, 0, 1, 1 ) ! not jpj-row 1319 1293 IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj ) 1320 ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1) 1321 ENDIF 1322 phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 1294 ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1) 1295 ENDIF 1296 phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 1323 1297 pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj) 1324 1298 END_2D … … 1331 1305 !! *** ROUTINE wad_sp *** 1332 1306 !! 1333 !! ** Purpose : 1307 !! ** Purpose : 1334 1308 !!---------------------------------------------------------------------- 1335 1309 INTEGER :: ji ,jj ! dummy loop indices … … 1364 1338 & MAX( pshn(ji,jj) , pshn(ji,jj+1) ) > & 1365 1339 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1366 1340 1367 1341 IF(ll_tmp1) THEN 1368 1342 zcpy(ji,jj) = 1.0_wp … … 1376 1350 ENDIF 1377 1351 END_2D 1378 1352 1379 1353 END SUBROUTINE wad_spg 1380 1354 1381 1355 1382 1356 … … 1384 1358 !!---------------------------------------------------------------------- 1385 1359 !! *** ROUTINE dyn_drg_init *** 1386 !! 1387 !! ** Purpose : - add the baroclinic top/bottom drag contribution to 1360 !! 1361 !! ** Purpose : - add the baroclinic top/bottom drag contribution to 1388 1362 !! the baroclinic part of the barotropic RHS 1389 1363 !! - compute the barotropic drag coefficients 1390 1364 !! 1391 !! ** Method : computation done over the INNER domain only 1365 !! ** Method : computation done over the INNER domain only 1392 1366 !!---------------------------------------------------------------------- 1393 1367 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices … … 1406 1380 ! 1407 1381 IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! top+bottom friction (ocean cavities) 1408 1382 1409 1383 DO_2D( 0, 0, 0, 0 ) 1410 1384 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) … … 1421 1395 ! 1422 1396 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities 1423 1397 1424 1398 DO_2D( 0, 0, 0, 0 ) 1425 ikbu = mbku(ji,jj) 1426 ikbv = mbkv(ji,jj) 1399 ikbu = mbku(ji,jj) 1400 ikbv = mbkv(ji,jj) 1427 1401 zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm) 1428 1402 zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm) 1429 1403 END_2D 1430 1404 ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities 1431 1405 1432 1406 DO_2D( 0, 0, 0, 0 ) 1433 ikbu = mbku(ji,jj) 1434 ikbv = mbkv(ji,jj) 1407 ikbu = mbku(ji,jj) 1408 ikbv = mbkv(ji,jj) 1435 1409 zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb) 1436 1410 zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb) … … 1441 1415 zztmp = -1._wp / rDt_e 1442 1416 DO_2D( 0, 0, 0, 0 ) 1443 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( & 1417 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( & 1444 1418 & r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) 1445 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( & 1419 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( & 1446 1420 & r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp ) 1447 1421 END_2D 1448 1422 ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation) 1449 1423 1450 1424 DO_2D( 0, 0, 0, 0 ) 1451 1425 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) … … 1459 1433 ! 1460 1434 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity 1461 1435 1462 1436 DO_2D( 0, 0, 0, 0 ) 1463 1437 iktu = miku(ji,jj) … … 1467 1441 END_2D 1468 1442 ELSE ! CENTRED integration: use BEFORE top baroclinic velocity 1469 1443 1470 1444 DO_2D( 0, 0, 0, 0 ) 1471 1445 iktu = miku(ji,jj) … … 1477 1451 ! 1478 1452 ! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation) 1479 1453 1480 1454 DO_2D( 0, 0, 0, 0 ) 1481 1455 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) … … 1498 1472 ! ! set Half-step back interpolation coefficient 1499 1473 IF ( jn==1 .AND. ll_init ) THEN !* Forward-backward 1500 za0 = 1._wp 1501 za1 = 0._wp 1474 za0 = 1._wp 1475 za1 = 0._wp 1502 1476 za2 = 0._wp 1503 1477 za3 = 0._wp … … 1506 1480 za1 =-0.1666666666666_wp ! za1 = gam 1507 1481 za2 = 0.0833333333333_wp ! za2 = eps 1508 za3 = 0._wp 1509 ELSE !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 1510 IF( rn_bt_alpha == 0._wp ) THEN ! Time diffusion 1482 za3 = 0._wp 1483 ELSE !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 1484 IF( rn_bt_alpha == 0._wp ) THEN ! Time diffusion 1511 1485 za0 = 0.614_wp ! za0 = 1/2 + gam + 2*eps 1512 1486 za1 = 0.285_wp ! za1 = 1/2 - 2*gam - 3*eps … … 1520 1494 za2 = zgamma 1521 1495 za3 = zepsilon 1522 ENDIF 1496 ENDIF 1523 1497 ENDIF 1524 1498 END SUBROUTINE ts_bck_interp
Note: See TracChangeset
for help on using the changeset viewer.
