Changeset 12229
- Timestamp:
- 2019-12-12T17:41:04+01:00 (3 years ago)
- Location:
- NEMO/branches/2019/dev_r11943_MERGE_2019
- Files:
-
- 17 edited
- 1 copied
-
cfgs/SHARED/axis_def_nemo.xml (modified) (1 diff)
-
cfgs/SHARED/field_def_nemo-oce.xml (modified) (6 diffs)
-
cfgs/SHARED/namelist_ref (modified) (1 diff)
-
src/NST/agrif_all_update.F90 (modified) (3 diffs)
-
src/NST/agrif_ice_update.F90 (modified) (3 diffs)
-
src/NST/agrif_oce.F90 (modified) (5 diffs)
-
src/NST/agrif_oce_interp.F90 (modified) (37 diffs)
-
src/NST/agrif_oce_sponge.F90 (modified) (26 diffs)
-
src/NST/agrif_oce_update.F90 (modified) (53 diffs)
-
src/NST/agrif_top_interp.F90 (modified) (6 diffs)
-
src/NST/agrif_top_sponge.F90 (modified) (8 diffs)
-
src/NST/agrif_top_update.F90 (modified) (15 diffs)
-
src/NST/agrif_user.F90 (modified) (3 diffs)
-
src/NST/vremap.F90 (copied) (copied from NEMO/branches/2019/dev_AGRIF-01-05_merged/src/NST/vremap.F90)
-
src/OCE/DIA/diamlr.F90 (modified) (1 diff)
-
src/OCE/DYN/dynspg_ts.F90 (modified) (3 diffs)
-
tests/VORTEX/EXPREF/1_context_nemo.xml (modified) (2 diffs)
-
tests/VORTEX/EXPREF/context_nemo.xml (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11943_MERGE_2019/cfgs/SHARED/axis_def_nemo.xml
r12193 r12229 22 22 <axis id="basin" long_name="Sub-basin mask (1=Global 2=Atlantic 3=Indo-Pacific 4=Indian, 5=Pacific)" unit="1" /> 23 23 <axis id="nstrait" long_name="Number of straits" unit="1" /> 24 <!-- ABL vertical axis definition --> 25 <axis id="ght_abl" long_name="ABL Vertical T levels" unit="m" positive="up" /> 26 <axis id="ghw_abl" long_name="ABL Vertical W levels" unit="m" positive="up" /> 24 27 <axis id="section" n_glo="16" /> 25 28 <axis id="section_ice" n_glo="4" /> -
NEMO/branches/2019/dev_r11943_MERGE_2019/cfgs/SHARED/field_def_nemo-oce.xml
r12205 r12229 127 127 <field id="toce_vmean300" field_ref="toce_e3t_vsum300" unit="degree_C" grid_ref="grid_T_vsum" detect_missing_value="true" > toce_e3t_vsum300/e3t_vsum300 </field> 128 128 129 <!-- AGRIF sponge --> 130 <field id="agrif_spt" long_name=" AGRIF t-sponge coefficient" unit=" " /> 129 131 130 132 <!-- t-eddy viscosity coefficients (ldfdyn) --> … … 420 422 </field_group> 421 423 422 <!-- scalar variables -->423 <field_group id="SBC_0D" grid_ref="grid_1point" >424 </field_group>425 424 426 425 </field_group> <!-- SBC --> … … 492 491 <field id="uoces" long_name="ocean transport along i-axis times salinity (CRS)" unit="1e-3*m/s" grid_ref="grid_U_3D" /> 493 492 <field id="ssuww" long_name="ocean surface wind work along i-axis" standard_name="surface_x_wind_work" unit="N/m*s" > utau * ssu </field> 494 493 <!-- AGRIF sponge --> 494 <field id="agrif_spu" long_name=" AGRIF u-sponge coefficient" unit=" " /> 495 495 <!-- u-eddy diffusivity coefficients (available if ln_traldf_OFF=F) --> 496 496 <field id="ahtu_2d" long_name=" surface u-eddy diffusivity coefficient" unit="m2/s or m4/s" /> … … 550 550 <field id="voces" long_name="ocean transport along j-axis times salinity (CRS)" unit="1e-3*m/s" grid_ref="grid_V_3D" /> 551 551 <field id="ssvww" long_name="ocean surface wind work along j-axis" standard_name="surface_y_wind_work" unit="N/m*s" > vtau * ssv </field> 552 552 <!-- AGRIF sponge --> 553 <field id="agrif_spv" long_name=" AGRIF v-sponge coefficient" unit=" " /> 553 554 <!-- v-eddy diffusivity coefficients (available if ln_traldf_OFF=F) --> 554 555 <field id="ahtv_2d" long_name=" surface v-eddy diffusivity coefficient" unit="m2/s or (m4/s)^1/2" /> … … 639 640 640 641 <!-- F grid --> 642 <!-- AGRIF sponge --> 643 <field id="agrif_spf" long_name=" AGRIF f-sponge coefficient" unit=" " /> 641 644 <!-- f-eddy viscosity coefficients (ldfdyn) --> 642 645 <field id="ahmf_2d" long_name=" surface f-eddy viscosity coefficient" unit="m2/s or m4/s" /> … … 691 694 692 695 693 <!-- variables available with key_float-->696 <!-- variables available with ln_floats --> 694 697 695 698 <field_group id="floatvar" grid_ref="grid_T_nfloat" operation="instant" > -
NEMO/branches/2019/dev_r11943_MERGE_2019/cfgs/SHARED/namelist_ref
r12205 r12229 633 633 &namagrif ! AGRIF zoom ("key_agrif") 634 634 !----------------------------------------------------------------------- 635 ln_agrif_2way = .true. ! activate two way nesting 635 636 ln_spc_dyn = .true. ! use 0 as special value for dynamics 636 637 rn_sponge_tra = 2880. ! coefficient for tracer sponge layer [m2/s] 637 638 rn_sponge_dyn = 2880. ! coefficient for dynamics sponge layer [m2/s] 639 rn_trelax_tra = 0.01 ! inverse of relaxation time (in steps) for tracers [] 640 rn_trelax_dyn = 0.01 ! inverse of relaxation time (in steps) for dynamics [] 638 641 ln_chk_bathy = .false. ! =T check the parent bathymetry 639 642 / -
NEMO/branches/2019/dev_r11943_MERGE_2019/src/NST/agrif_all_update.F90
r10069 r12229 1 #define TWO_WAY 2 3 MODULE agrif_all_update 1 MODULE agrif_all_update 4 2 !!====================================================================== 5 3 !! *** MODULE agrif_all_update *** … … 41 39 !! Order of update matters here ! 42 40 !!---------------------------------------------------------------------- 43 # if defined TWO_WAY 44 IF (Agrif_Root()) RETURN 41 IF (( .NOT.ln_agrif_2way ).OR.(Agrif_Root())) RETURN 45 42 ! 46 43 IF (lwp.AND.lk_agrif_debug) Write(*,*) ' --> START AGRIF UPDATE from grid Number',Agrif_Fixed() … … 67 64 ! 68 65 Agrif_UseSpecialValueInUpdate = .FALSE. 69 #endif70 66 END SUBROUTINE agrif_Update_All 71 67 -
NEMO/branches/2019/dev_r11943_MERGE_2019/src/NST/agrif_ice_update.F90
r10069 r12229 1 #define TWO_WAY2 !!#undef TWO_WAY3 1 #undef DECAL_FEEDBACK /* SEPARATION of INTERFACES*/ 4 2 … … 63 61 Agrif_UseSpecialValueInUpdate = .TRUE. 64 62 65 # if defined TWO_WAY66 63 # if ! defined DECAL_FEEDBACK 67 64 CALL Agrif_Update_Variable( tra_ice_id , procname = update_tra_ice ) … … 79 76 ! CALL Agrif_Update_Variable( u_ice_id , locupdate=(/0,1/), procname = update_u_ice ) 80 77 ! CALL Agrif_Update_Variable( v_ice_id , locupdate=(/0,1/), procname = update_v_ice ) 81 # endif82 78 Agrif_SpecialValueFineGrid = 0. 83 79 Agrif_UseSpecialValueInUpdate = .FALSE. -
NEMO/branches/2019/dev_r11943_MERGE_2019/src/NST/agrif_oce.F90
r11949 r12229 17 17 18 18 PUBLIC agrif_oce_alloc ! routine called by nemo_init in nemogcm.F90 19 #if defined key_vertical 20 PUBLIC reconstructandremap ! remapping routine 21 #endif 19 22 20 ! !!* Namelist namagrif: AGRIF parameters 23 LOGICAL , PUBLIC :: ln_spc_dyn = .FALSE. !: 24 INTEGER , PUBLIC, PARAMETER :: nn_sponge_len = 2 !: Sponge width (in number of parent grid points) 21 LOGICAL , PUBLIC :: ln_agrif_2way = .TRUE. !: activate two way nesting 22 LOGICAL , PUBLIC :: ln_spc_dyn = .FALSE. !: use zeros (.false.) or not (.true.) in 23 !: bdys dynamical fields interpolation 25 24 REAL(wp), PUBLIC :: rn_sponge_tra = 2800. !: sponge coeff. for tracers 26 25 REAL(wp), PUBLIC :: rn_sponge_dyn = 2800. !: sponge coeff. for dynamics 26 REAL(wp), PUBLIC :: rn_trelax_tra = 0.01 !: time relaxation parameter for tracers 27 REAL(wp), PUBLIC :: rn_trelax_dyn = 0.01 !: time relaxation parameter for momentum 27 28 LOGICAL , PUBLIC :: ln_chk_bathy = .FALSE. !: check of parent bathymetry 28 LOGICAL , PUBLIC :: lk_agrif_clp = .FALSE. !: Force clamped bcs 29 ! !!! OLD namelist names 30 REAL(wp), PUBLIC :: visc_tra !: sponge coeff. for tracers 31 REAL(wp), PUBLIC :: visc_dyn !: sponge coeff. for dynamics 32 29 ! 30 INTEGER , PUBLIC, PARAMETER :: nn_sponge_len = 2 !: Sponge width (in number of parent grid points) 33 31 LOGICAL , PUBLIC :: spongedoneT = .FALSE. !: tracer sponge layer indicator 34 32 LOGICAL , PUBLIC :: spongedoneU = .FALSE. !: dynamics sponge layer indicator … … 42 40 LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_u 43 41 LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_v 44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fsaht_spu, fsaht_spv !: sponge diffusivities 45 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fsahm_spt, fsahm_spf !: sponge viscosities 42 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: utint_stage 43 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vtint_stage 44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fspu, fspv !: sponge arrays 45 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fspt, fspf !: " " 46 46 47 47 ! Barotropic arrays used to store open boundary data during time-splitting loop: 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ubdy_w, vbdy_w, hbdy_w 49 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ubdy_e, vbdy_e, hbdy_e 50 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ubdy_n, vbdy_n, hbdy_n 51 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ubdy_s, vbdy_s, hbdy_s 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ubdy, vbdy, hbdy 52 49 INTEGER , PUBLIC, SAVE :: Kbb_a, Kmm_a, Krhs_a !: AGRIF module-specific copies of time-level indices 53 50 51 # if defined key_vertical 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht0_parent, hu0_parent, hv0_parent 53 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt_parent, mbku_parent, mbkv_parent 54 # endif 54 55 55 56 INTEGER, PUBLIC :: tsn_id ! AGRIF profile for tracers interpolation and update … … 65 66 INTEGER, PUBLIC :: scales_t_id 66 67 INTEGER, PUBLIC :: avt_id, avm_id, en_id ! TKE related identificators 67 INTEGER, PUBLIC :: umsk_id, vmsk_id68 INTEGER, PUBLIC :: mbkt_id, ht0_id 68 69 INTEGER, PUBLIC :: kindic_agr 69 70 … … 83 84 ierr(:) = 0 84 85 ! 85 ALLOCATE( fsaht_spu(jpi,jpj), fsaht_spv(jpi,jpj), & 86 & fsahm_spt(jpi,jpj), fsahm_spf(jpi,jpj), & 87 & tabspongedone_tsn(jpi,jpj), & 86 ALLOCATE( fspu(jpi,jpj), fspv(jpi,jpj), & 87 & fspt(jpi,jpj), fspf(jpi,jpj), & 88 & tabspongedone_tsn(jpi,jpj), & 89 & utint_stage(jpi,jpj), vtint_stage(jpi,jpj), & 88 90 # if defined key_top 89 91 & tabspongedone_trn(jpi,jpj), & 90 # endif 92 # endif 93 # if defined key_vertical 94 & ht0_parent(jpi,jpj), mbkt_parent(jpi,jpj), & 95 & hu0_parent(jpi,jpj), mbku_parent(jpi,jpj), & 96 & hv0_parent(jpi,jpj), mbkv_parent(jpi,jpj), & 97 # endif 91 98 & tabspongedone_u (jpi,jpj), & 92 99 & tabspongedone_v (jpi,jpj), STAT = ierr(1) ) 93 100 94 ALLOCATE( ubdy_w(nbghostcells,jpj), vbdy_w(nbghostcells,jpj), hbdy_w(nbghostcells,jpj), & 95 & ubdy_e(nbghostcells,jpj), vbdy_e(nbghostcells,jpj), hbdy_e(nbghostcells,jpj), & 96 & ubdy_n(jpi,nbghostcells), vbdy_n(jpi,nbghostcells), hbdy_n(jpi,nbghostcells), & 97 & ubdy_s(jpi,nbghostcells), vbdy_s(jpi,nbghostcells), hbdy_s(jpi,nbghostcells), STAT = ierr(2) ) 101 ALLOCATE( ubdy(jpi,jpj), vbdy(jpi,jpj), hbdy(jpi,jpj), STAT = ierr(2) ) 98 102 99 103 agrif_oce_alloc = MAXVAL(ierr) … … 101 105 END FUNCTION agrif_oce_alloc 102 106 103 #if defined key_vertical104 SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout)105 !!----------------------------------------------------------------------106 !! *** FUNCTION reconstructandremap ***107 !!----------------------------------------------------------------------108 IMPLICIT NONE109 INTEGER N, Nout110 REAL(wp) tabin(N), tabout(Nout)111 REAL(wp) hin(N), hout(Nout)112 REAL(wp) coeffremap(N,3),zwork(N,3)113 REAL(wp) zwork2(N+1,3)114 INTEGER jk115 DOUBLE PRECISION, PARAMETER :: dsmll=1.0d-8116 REAL(wp) q,q01,q02,q001,q002,q0117 REAL(wp) z_win(1:N+1), z_wout(1:Nout+1)118 REAL(wp),PARAMETER :: dpthin = 1.D-3119 INTEGER :: k1, kbox, ktop, ka, kbot120 REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop121 122 z_win(1)=0.; z_wout(1)= 0.123 DO jk=1,N124 z_win(jk+1)=z_win(jk)+hin(jk)125 ENDDO126 127 DO jk=1,Nout128 z_wout(jk+1)=z_wout(jk)+hout(jk)129 ENDDO130 131 DO jk=2,N132 zwork(jk,1)=1./(hin(jk-1)+hin(jk))133 ENDDO134 135 DO jk=2,N-1136 q0 = 1./(hin(jk-1)+hin(jk)+hin(jk+1))137 zwork(jk,2)=hin(jk-1)+2.*hin(jk)+hin(jk+1)138 zwork(jk,3)=q0139 ENDDO140 141 DO jk= 2,N142 zwork2(jk,1)=zwork(jk,1)*(tabin(jk)-tabin(jk-1))143 ENDDO144 145 coeffremap(:,1) = tabin(:)146 147 DO jk=2,N-1148 q001 = hin(jk)*zwork2(jk+1,1)149 q002 = hin(jk)*zwork2(jk,1)150 IF (q001*q002 < 0) then151 q001 = 0.152 q002 = 0.153 ENDIF154 q=zwork(jk,2)155 q01=q*zwork2(jk+1,1)156 q02=q*zwork2(jk,1)157 IF (abs(q001) > abs(q02)) q001 = q02158 IF (abs(q002) > abs(q01)) q002 = q01159 160 q=(q001-q002)*zwork(jk,3)161 q001=q001-q*hin(jk+1)162 q002=q002+q*hin(jk-1)163 164 coeffremap(jk,3)=coeffremap(jk,1)+q001165 coeffremap(jk,2)=coeffremap(jk,1)-q002166 167 zwork2(jk,1)=(2.*q001-q002)**2168 zwork2(jk,2)=(2.*q002-q001)**2169 ENDDO170 171 DO jk=1,N172 IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN173 coeffremap(jk,3) = coeffremap(jk,1)174 coeffremap(jk,2) = coeffremap(jk,1)175 zwork2(jk,1) = 0.176 zwork2(jk,2) = 0.177 ENDIF178 ENDDO179 180 DO jk=2,N181 q002=max(zwork2(jk-1,2),dsmll)182 q001=max(zwork2(jk,1),dsmll)183 zwork2(jk,3)=(q001*coeffremap(jk-1,3)+q002*coeffremap(jk,2))/(q001+q002)184 ENDDO185 186 zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3)187 zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3)188 189 DO jk=1,N190 q01=zwork2(jk+1,3)-coeffremap(jk,1)191 q02=coeffremap(jk,1)-zwork2(jk,3)192 q001=2.*q01193 q002=2.*q02194 IF (q01*q02<0) then195 q01=0.196 q02=0.197 ELSEIF (abs(q01)>abs(q002)) then198 q01=q002199 ELSEIF (abs(q02)>abs(q001)) then200 q02=q001201 ENDIF202 coeffremap(jk,2)=coeffremap(jk,1)-q02203 coeffremap(jk,3)=coeffremap(jk,1)+q01204 ENDDO205 206 zbot=0.0207 kbot=1208 DO jk=1,Nout209 ztop=zbot !top is bottom of previous layer210 ktop=kbot211 IF (ztop.GE.z_win(ktop+1)) then212 ktop=ktop+1213 ENDIF214 215 zbot=z_wout(jk+1)216 zthk=zbot-ztop217 218 IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(Nout+1)) THEN219 220 kbot=ktop221 DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.N)222 kbot=kbot+1223 ENDDO224 zbox=zbot225 DO k1= jk+1,Nout226 IF (z_wout(k1+1)-z_wout(k1).GT.dpthin) THEN227 exit !thick layer228 ELSE229 zbox=z_wout(k1+1) !include thin adjacent layers230 IF(zbox.EQ.z_wout(Nout+1)) THEN231 exit !at bottom232 ENDIF233 ENDIF234 ENDDO235 zthk=zbox-ztop236 237 kbox=ktop238 DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.N)239 kbox=kbox+1240 ENDDO241 242 IF(ktop.EQ.kbox) THEN243 IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN244 IF(hin(kbox).GT.dpthin) THEN245 q001 = (zbox-z_win(kbox))/hin(kbox)246 q002 = (ztop-z_win(kbox))/hin(kbox)247 q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002)248 q02=q01-1.+(q001+q002)249 q0=1.-q01-q02250 ELSE251 q0 = 1.0252 q01 = 0.253 q02 = 0.254 ENDIF255 tabout(jk)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3)256 ELSE257 tabout(jk) = tabin(kbox)258 ENDIF259 ELSE260 IF(ktop.LE.jk .AND. kbox.GE.jk) THEN261 ka = jk262 ELSEIF (kbox-ktop.GE.3) THEN263 ka = (kbox+ktop)/2264 ELSEIF (hin(ktop).GE.hin(kbox)) THEN265 ka = ktop266 ELSE267 ka = kbox268 ENDIF !choose ka269 270 offset=coeffremap(ka,1)271 272 qtop = z_win(ktop+1)-ztop !partial layer thickness273 IF(hin(ktop).GT.dpthin) THEN274 q=(ztop-z_win(ktop))/hin(ktop)275 q01=q*(q-1.)276 q02=q01+q277 q0=1-q01-q02278 ELSE279 q0 = 1.280 q01 = 0.281 q02 = 0.282 ENDIF283 284 tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop285 286 DO k1= ktop+1,kbox-1287 tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1)288 ENDDO !k1289 290 qbot = zbox-z_win(kbox) !partial layer thickness291 IF(hin(kbox).GT.dpthin) THEN292 q=qbot/hin(kbox)293 q01=(q-1.)**2294 q02=q01-1.+q295 q0=1-q01-q02296 ELSE297 q0 = 1.0298 q01 = 0.299 q02 = 0.300 ENDIF301 302 tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot303 304 rpsum=1.0d0/zthk305 tabout(jk)=offset+tsum*rpsum306 307 ENDIF !single or multiple layers308 ELSE309 IF (jk==1) THEN310 write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1)311 ENDIF312 tabout(jk) = tabout(jk-1)313 314 ENDIF !normal:thin layer315 ENDDO !jk316 317 return318 end subroutine reconstructandremap319 #endif320 321 107 #endif 322 108 !!====================================================================== -
NEMO/branches/2019/dev_r11943_MERGE_2019/src/NST/agrif_oce_interp.F90
r11949 r12229 33 33 USE agrif_oce_sponge 34 34 USE lib_mpp 35 USE vremap 35 36 36 37 IMPLICIT NONE 37 38 PRIVATE 38 39 39 PUBLIC Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_ ssh_ts, Agrif_dta_ts40 PUBLIC Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_dyn_ts_flux, Agrif_ssh_ts, Agrif_dta_ts 40 41 PUBLIC Agrif_tra, Agrif_avm 41 42 PUBLIC interpun , interpvn 42 43 PUBLIC interptsn, interpsshn, interpavm 43 44 PUBLIC interpunb, interpvnb , interpub2b, interpvb2b 44 PUBLIC interpe3t, interpumsk, interpvmsk 45 45 PUBLIC interpe3t 46 #if defined key_vertical 47 PUBLIC interpht0, interpmbkt 48 # endif 46 49 INTEGER :: bdy_tinterp = 0 47 50 … … 78 81 ! 79 82 INTEGER :: ji, jj, jk ! dummy loop indices 80 INTEGER :: j1, j2, i1, i281 83 INTEGER :: ibdy1, jbdy1, ibdy2, jbdy2 82 84 REAL(wp), DIMENSION(jpi,jpj) :: zub, zvb … … 93 95 Agrif_UseSpecialValue = .FALSE. 94 96 ! 95 ! prevent smoothing in ghost cells96 i1 = 1 ; i2 = nlci97 j1 = 1 ; j2 = nlcj98 IF( nbondj == -1 .OR. nbondj == 2 ) j1 = 2 + nbghostcells99 IF( nbondj == +1 .OR. nbondj == 2 ) j2 = nlcj - nbghostcells - 1100 IF( nbondi == -1 .OR. nbondi == 2 ) i1 = 2 + nbghostcells101 IF( nbondi == +1 .OR. nbondi == 2 ) i2 = nlci - nbghostcells - 1102 103 97 ! --- West --- ! 104 IF( nbondi == -1 .OR. nbondi == 2 ) THEN 105 ibdy1 = 2 106 ibdy2 = 1+nbghostcells 107 ! 108 IF( .NOT.ln_dynspg_ts ) THEN ! Store transport 109 uu_b(ibdy1:ibdy2,:,Krhs_a) = 0._wp 98 ibdy1 = 2 99 ibdy2 = 1+nbghostcells 100 ! 101 IF( .NOT.ln_dynspg_ts ) THEN ! Store transport 102 DO ji = mi0(ibdy1), mi1(ibdy2) 103 uu_b(ji,:,Krhs_a) = 0._wp 104 110 105 DO jk = 1, jpkm1 111 106 DO jj = 1, jpj 112 uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj, Krhs_a) & 113 & + e3u(ibdy1:ibdy2,jj,jk,Krhs_a) & 114 & * uu(ibdy1:ibdy2,jj,jk,Krhs_a) * umask(ibdy1:ibdy2,jj,jk) 115 END DO 116 END DO 107 uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 108 END DO 109 END DO 110 117 111 DO jj = 1, jpj 118 uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj,Krhs_a) * r1_hu(ibdy1:ibdy2,jj,Krhs_a) 119 END DO 120 ENDIF 121 ! 122 IF( .NOT.lk_agrif_clp ) THEN 123 DO jk=1,jpkm1 ! Smooth 124 DO jj=j1,j2 125 uu(ibdy2,jj,jk,Krhs_a) = 0.25_wp*( uu(ibdy2-1,jj,jk,Krhs_a)+2._wp*uu(ibdy2,jj,jk,Krhs_a) & 126 & + uu(ibdy2+1,jj,jk,Krhs_a) ) 127 END DO 128 END DO 129 ENDIF 130 ! 131 zub(ibdy1:ibdy2,:) = 0._wp ! Correct transport 112 uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) * r1_hu(ji,jj,Krhs_a) 113 END DO 114 END DO 115 ENDIF 116 ! 117 DO ji = mi0(ibdy1), mi1(ibdy2) 118 zub(ji,:) = 0._wp ! Correct transport 132 119 DO jk = 1, jpkm1 133 120 DO jj = 1, jpj 134 zub( ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) + e3u(ibdy1:ibdy2,jj,jk,Krhs_a)&135 & * uu(ibdy1:ibdy2,jj,jk,Krhs_a) * umask(ibdy1:ibdy2,jj,jk)121 zub(ji,jj) = zub(ji,jj) & 122 & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a)*umask(ji,jj,jk) 136 123 END DO 137 124 END DO 138 125 DO jj=1,jpj 139 zub( ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu(ibdy1:ibdy2,jj,Krhs_a)126 zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 140 127 END DO 141 128 142 129 DO jk = 1, jpkm1 143 130 DO jj = 1, jpj 144 uu(ibdy1:ibdy2,jj,jk,Krhs_a) = ( uu(ibdy1:ibdy2,jj,jk,Krhs_a) & 145 & + uu_b(ibdy1:ibdy2,jj, Krhs_a) & 146 & - zub(ibdy1:ibdy2,jj) ) & 147 & * umask(ibdy1:ibdy2,jj,jk) 148 END DO 149 END DO 131 uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) + uu_b(ji,jj,Krhs_a)-zub(ji,jj)) * umask(ji,jj,jk) 132 END DO 133 END DO 134 END DO 150 135 151 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 152 zvb(ibdy1:ibdy2,:) = 0._wp 136 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 137 DO ji = mi0(ibdy1), mi1(ibdy2) 138 zvb(ji,:) = 0._wp 153 139 DO jk = 1, jpkm1 154 140 DO jj = 1, jpj 155 zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) + e3v(ibdy1:ibdy2,jj,jk,Krhs_a) & 156 & * vv(ibdy1:ibdy2,jj,jk,Krhs_a) * vmask(ibdy1:ibdy2,jj,jk) 141 zvb(ji,jj) = zvb(ji,jj) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 157 142 END DO 158 143 END DO 159 144 DO jj = 1, jpj 160 zvb( ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv(ibdy1:ibdy2,jj,Krhs_a)145 zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 161 146 END DO 162 147 DO jk = 1, jpkm1 163 148 DO jj = 1, jpj 164 vv(ibdy1:ibdy2,jj,jk,Krhs_a) = ( vv(ibdy1:ibdy2,jj,jk,Krhs_a) & 165 & + vv_b(ibdy1:ibdy2,jj, Krhs_a) & 166 & - zvb(ibdy1:ibdy2,jj) ) & 167 & * vmask(ibdy1:ibdy2,jj,jk) 168 END DO 169 END DO 170 ENDIF 171 ! 172 DO jk = 1, jpkm1 ! Mask domain edges 173 DO jj = 1, jpj 174 uu(1,jj,jk,Krhs_a) = 0._wp 175 vv(1,jj,jk,Krhs_a) = 0._wp 176 END DO 177 END DO 149 vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) + vv_b(ji,jj,Krhs_a)-zvb(ji,jj))*vmask(ji,jj,jk) 150 END DO 151 END DO 152 END DO 178 153 ENDIF 179 154 180 155 ! --- East --- ! 181 IF( nbondi == 1 .OR. nbondi == 2 ) THEN182 ibdy1 = nlci-1-nbghostcells183 ibdy2 = nlci-2184 !185 IF( .NOT.ln_dynspg_ts ) THEN ! Store transport186 uu_b( ibdy1:ibdy2,:,Krhs_a) = 0._wp156 ibdy1 = jpiglo-1-nbghostcells 157 ibdy2 = jpiglo-2 158 ! 159 IF( .NOT.ln_dynspg_ts ) THEN ! Store transport 160 DO ji = mi0(ibdy1), mi1(ibdy2) 161 uu_b(ji,:,Krhs_a) = 0._wp 187 162 DO jk = 1, jpkm1 188 163 DO jj = 1, jpj 189 uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj, Krhs_a) & 190 & + e3u(ibdy1:ibdy2,jj,jk,Krhs_a) & 191 & * uu(ibdy1:ibdy2,jj,jk,Krhs_a) & 192 & * umask(ibdy1:ibdy2,jj,jk) 164 uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) & 165 & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 193 166 END DO 194 167 END DO 195 168 DO jj = 1, jpj 196 uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj,Krhs_a) * r1_hu(ibdy1:ibdy2,jj,Krhs_a) 197 END DO 198 ENDIF 199 ! 200 IF( .NOT.lk_agrif_clp ) THEN 201 DO jk=1,jpkm1 ! Smooth 202 DO jj=j1,j2 203 uu(ibdy1,jj,jk,Krhs_a) = 0.25_wp*( uu(ibdy1-1,jj,jk,Krhs_a) & 204 & + 2._wp*uu(ibdy1 ,jj,jk,Krhs_a) & 205 & + uu(ibdy1+1,jj,jk,Krhs_a) ) 206 END DO 207 END DO 208 ENDIF 209 ! 210 zub(ibdy1:ibdy2,:) = 0._wp ! Correct transport 169 uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) * r1_hu(ji,jj,Krhs_a) 170 END DO 171 END DO 172 ENDIF 173 ! 174 DO ji = mi0(ibdy1), mi1(ibdy2) 175 zub(ji,:) = 0._wp ! Correct transport 211 176 DO jk = 1, jpkm1 212 177 DO jj = 1, jpj 213 zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) & 214 & + e3u(ibdy1:ibdy2,jj,jk,Krhs_a) & 215 & * uu(ibdy1:ibdy2,jj,jk,Krhs_a) * umask(ibdy1:ibdy2,jj,jk) 178 zub(ji,jj) = zub(ji,jj) & 179 & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 216 180 END DO 217 181 END DO 218 182 DO jj=1,jpj 219 zub( ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu(ibdy1:ibdy2,jj,Krhs_a)183 zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 220 184 END DO 221 185 222 186 DO jk = 1, jpkm1 223 187 DO jj = 1, jpj 224 uu(ibdy1:ibdy2,jj,jk,Krhs_a) = ( uu(ibdy1:ibdy2,jj,jk,Krhs_a) & 225 & + uu_b(ibdy1:ibdy2,jj, Krhs_a) & 226 & - zub(ibdy1:ibdy2,jj) & 227 & ) * umask(ibdy1:ibdy2,jj,jk) 228 END DO 229 END DO 188 uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) & 189 & + uu_b(ji,jj,Krhs_a)-zub(ji,jj))*umask(ji,jj,jk) 190 END DO 191 END DO 192 END DO 230 193 231 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 232 ibdy1 = ibdy1 + 1 233 ibdy2 = ibdy2 + 1 234 zvb(ibdy1:ibdy2,:) = 0._wp 194 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 195 ibdy1 = jpiglo-nbghostcells 196 ibdy2 = jpiglo-1 197 DO ji = mi0(ibdy1), mi1(ibdy2) 198 zvb(ji,:) = 0._wp 235 199 DO jk = 1, jpkm1 236 200 DO jj = 1, jpj 237 zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) & 238 & + e3v(ibdy1:ibdy2,jj,jk,Krhs_a) & 239 & * vv(ibdy1:ibdy2,jj,jk,Krhs_a) & 240 & * vmask(ibdy1:ibdy2,jj,jk) 201 zvb(ji,jj) = zvb(ji,jj) & 202 & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 241 203 END DO 242 204 END DO 243 205 DO jj = 1, jpj 244 zvb( ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv(ibdy1:ibdy2,jj,Krhs_a)206 zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 245 207 END DO 246 208 DO jk = 1, jpkm1 247 209 DO jj = 1, jpj 248 vv(ibdy1:ibdy2,jj,jk,Krhs_a) = ( vv(ibdy1:ibdy2,jj,jk,Krhs_a) & 249 & + vv_b(ibdy1:ibdy2,jj, Krhs_a) & 250 & - zvb(ibdy1:ibdy2,jj) & 251 & ) * vmask(ibdy1:ibdy2,jj,jk) 252 END DO 253 END DO 254 ENDIF 255 ! 256 DO jk = 1, jpkm1 ! Mask domain edges 257 DO jj = 1, jpj 258 uu(nlci-1,jj,jk,Krhs_a) = 0._wp 259 vv(nlci ,jj,jk,Krhs_a) = 0._wp 260 END DO 261 END DO 210 vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) & 211 & + vv_b(ji,jj,Krhs_a)-zvb(ji,jj)) * vmask(ji,jj,jk) 212 END DO 213 END DO 214 END DO 262 215 ENDIF 263 216 264 217 ! --- South --- ! 265 IF( nbondj == -1 .OR. nbondj == 2 ) THEN266 jbdy1 = 2267 jbdy2 = 1+nbghostcells268 !269 IF( .NOT.ln_dynspg_ts ) THEN ! Store transport270 vv_b(:,j bdy1:jbdy2,Krhs_a) = 0._wp218 jbdy1 = 2 219 jbdy2 = 1+nbghostcells 220 ! 221 IF( .NOT.ln_dynspg_ts ) THEN ! Store transport 222 DO jj = mj0(jbdy1), mj1(jbdy2) 223 vv_b(:,jj,Krhs_a) = 0._wp 271 224 DO jk = 1, jpkm1 272 225 DO ji = 1, jpi 273 vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2, Krhs_a) & 274 & + e3v(ji,jbdy1:jbdy2,jk,Krhs_a) & 275 & * vv(ji,jbdy1:jbdy2,jk,Krhs_a) & 276 & * vmask(ji,jbdy1:jbdy2,jk) 226 vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) & 227 & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 277 228 END DO 278 229 END DO 279 230 DO ji=1,jpi 280 vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2,Krhs_a) * r1_hv(ji,jbdy1:jbdy2,Krhs_a) 281 END DO 282 ENDIF 283 ! 284 IF ( .NOT.lk_agrif_clp ) THEN 285 DO jk = 1, jpkm1 ! Smooth 286 DO ji = i1, i2 287 vv(ji,jbdy2,jk,Krhs_a) = 0.25_wp*( vv(ji,jbdy2-1,jk,Krhs_a) & 288 & + 2._wp*vv(ji,jbdy2 ,jk,Krhs_a) & 289 & + vv(ji,jbdy2+1,jk,Krhs_a) ) 290 END DO 291 END DO 292 ENDIF 293 ! 294 zvb(:,jbdy1:jbdy2) = 0._wp ! Correct transport 231 vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) * r1_hv(ji,jj,Krhs_a) 232 END DO 233 END DO 234 ENDIF 235 ! 236 DO jj = mj0(jbdy1), mj1(jbdy2) 237 zvb(:,jj) = 0._wp ! Correct transport 295 238 DO jk=1,jpkm1 296 239 DO ji=1,jpi 297 zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) & 298 & + e3v(ji,jbdy1:jbdy2,jk,Krhs_a) & 299 & * vv(ji,jbdy1:jbdy2,jk,Krhs_a) & 300 & * vmask(ji,jbdy1:jbdy2,jk) 240 zvb(ji,jj) = zvb(ji,jj) & 241 & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 301 242 END DO 302 243 END DO 303 244 DO ji = 1, jpi 304 zvb(ji,j bdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv(ji,jbdy1:jbdy2,Krhs_a)245 zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 305 246 END DO 306 247 307 248 DO jk = 1, jpkm1 308 249 DO ji = 1, jpi 309 vv(ji,jbdy1:jbdy2,jk,Krhs_a) = ( vv(ji,jbdy1:jbdy2,jk,Krhs_a) & 310 & + vv_b(ji,jbdy1:jbdy2, Krhs_a) & 311 & - zvb(ji,jbdy1:jbdy2) & 312 & ) * vmask(ji,jbdy1:jbdy2,jk) 313 END DO 314 END DO 250 vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) & 251 & + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) ) * vmask(ji,jj,jk) 252 END DO 253 END DO 254 END DO 315 255 316 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 317 zub(:,jbdy1:jbdy2) = 0._wp 256 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 257 DO jj = mj0(jbdy1), mj1(jbdy2) 258 zub(:,jj) = 0._wp 318 259 DO jk = 1, jpkm1 319 260 DO ji = 1, jpi 320 zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) & 321 & + e3u(ji,jbdy1:jbdy2,jk,Krhs_a) & 322 & * uu(ji,jbdy1:jbdy2,jk,Krhs_a) & 323 & * umask(ji,jbdy1:jbdy2,jk) 261 zub(ji,jj) = zub(ji,jj) & 262 & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 324 263 END DO 325 264 END DO 326 265 DO ji = 1, jpi 327 zub(ji,j bdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu(ji,jbdy1:jbdy2,Krhs_a)266 zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 328 267 END DO 329 268 330 269 DO jk = 1, jpkm1 331 270 DO ji = 1, jpi 332 uu(ji,jbdy1:jbdy2,jk,Krhs_a) = ( uu(ji,jbdy1:jbdy2,jk,Krhs_a) & 333 & + uu_b(ji,jbdy1:jbdy2, Krhs_a) & 334 & - zub(ji,jbdy1:jbdy2) & 335 & ) * umask(ji,jbdy1:jbdy2,jk) 336 END DO 337 END DO 338 ENDIF 339 ! 340 DO jk = 1, jpkm1 ! Mask domain edges 341 DO ji = 1, jpi 342 uu(ji,1,jk,Krhs_a) = 0._wp 343 vv(ji,1,jk,Krhs_a) = 0._wp 344 END DO 345 END DO 271 uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) & 272 & + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk) 273 END DO 274 END DO 275 END DO 346 276 ENDIF 347 277 348 278 ! --- North --- ! 349 IF( nbondj == 1 .OR. nbondj == 2 ) THEN350 jbdy1 = nlcj-1-nbghostcells351 jbdy2 = nlcj-2352 !353 IF( .NOT.ln_dynspg_ts ) THEN ! Store transport354 vv_b(:,j bdy1:jbdy2,Krhs_a) = 0._wp279 jbdy1 = jpjglo-1-nbghostcells 280 jbdy2 = jpjglo-2 281 ! 282 IF( .NOT.ln_dynspg_ts ) THEN ! Store transport 283 DO jj = mj0(jbdy1), mj1(jbdy2) 284 vv_b(:,jj,Krhs_a) = 0._wp 355 285 DO jk = 1, jpkm1 356 286 DO ji = 1, jpi 357 vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2, Krhs_a) & 358 & + e3v(ji,jbdy1:jbdy2,jk,Krhs_a) & 359 & * vv(ji,jbdy1:jbdy2,jk,Krhs_a) & 360 & * vmask(ji,jbdy1:jbdy2,jk) 287 vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) & 288 & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 361 289 END DO 362 290 END DO 363 291 DO ji=1,jpi 364 vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2,Krhs_a) * r1_hv(ji,jbdy1:jbdy2,Krhs_a) 365 END DO 366 ENDIF 367 ! 368 IF ( .NOT.lk_agrif_clp ) THEN 369 DO jk = 1, jpkm1 ! Smooth 370 DO ji = i1, i2 371 vv(ji,jbdy1,jk,Krhs_a) = 0.25_wp*( vv(ji,jbdy1-1,jk,Krhs_a) & 372 & + 2._wp*vv(ji,jbdy1 ,jk,Krhs_a) & 373 & + vv(ji,jbdy1+1,jk,Krhs_a) ) 374 END DO 375 END DO 376 ENDIF 377 ! 378 zvb(:,jbdy1:jbdy2) = 0._wp ! Correct transport 292 vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) * r1_hv(ji,jj,Krhs_a) 293 END DO 294 END DO 295 ENDIF 296 ! 297 DO jj = mj0(jbdy1), mj1(jbdy2) 298 zvb(:,jj) = 0._wp ! Correct transport 379 299 DO jk=1,jpkm1 380 300 DO ji=1,jpi 381 zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) & 382 & + e3v(ji,jbdy1:jbdy2,jk,Krhs_a) & 383 & * vv(ji,jbdy1:jbdy2,jk,Krhs_a) & 384 & * vmask(ji,jbdy1:jbdy2,jk) 301 zvb(ji,jj) = zvb(ji,jj) & 302 & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 385 303 END DO 386 304 END DO 387 305 DO ji = 1, jpi 388 zvb(ji,j bdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv(ji,jbdy1:jbdy2,Krhs_a)306 zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 389 307 END DO 390 308 391 309 DO jk = 1, jpkm1 392 310 DO ji = 1, jpi 393 vv(ji,jbdy1:jbdy2,jk,Krhs_a) = ( vv(ji,jbdy1:jbdy2,jk,Krhs_a) & 394 & + vv_b(ji,jbdy1:jbdy2, Krhs_a) & 395 & - zvb(ji,jbdy1:jbdy2) & 396 & ) * vmask(ji,jbdy1:jbdy2,jk) 397 END DO 398 END DO 311 vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) & 312 & + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) ) * vmask(ji,jj,jk) 313 END DO 314 END DO 315 END DO 399 316 400 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 401 jbdy1 = jbdy1 + 1 402 jbdy2 = jbdy2 + 1 403 zub(:,jbdy1:jbdy2) = 0._wp 317 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 318 jbdy1 = jpjglo-nbghostcells 319 jbdy2 = jpjglo-1 320 DO jj = mj0(jbdy1), mj1(jbdy2) 321 zub(:,jj) = 0._wp 404 322 DO jk = 1, jpkm1 405 323 DO ji = 1, jpi 406 zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) & 407 & + e3u(ji,jbdy1:jbdy2,jk,Krhs_a) & 408 & * uu(ji,jbdy1:jbdy2,jk,Krhs_a) & 409 & * umask(ji,jbdy1:jbdy2,jk) 324 zub(ji,jj) = zub(ji,jj) & 325 & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 410 326 END DO 411 327 END DO 412 328 DO ji = 1, jpi 413 zub(ji,j bdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu(ji,jbdy1:jbdy2,Krhs_a)329 zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 414 330 END DO 415 331 416 332 DO jk = 1, jpkm1 417 333 DO ji = 1, jpi 418 uu(ji,jbdy1:jbdy2,jk,Krhs_a) = ( uu(ji,jbdy1:jbdy2,jk,Krhs_a) & 419 & + uu_b(ji,jbdy1:jbdy2, Krhs_a) & 420 & - zub(ji,jbdy1:jbdy2) & 421 & ) * umask(ji,jbdy1:jbdy2,jk) 422 END DO 423 END DO 424 ENDIF 425 ! 426 DO jk = 1, jpkm1 ! Mask domain edges 427 DO ji = 1, jpi 428 uu(ji,nlcj ,jk,Krhs_a) = 0._wp 429 vv(ji,nlcj-1,jk,Krhs_a) = 0._wp 430 END DO 431 END DO 334 uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) & 335 & + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk) 336 END DO 337 END DO 338 END DO 432 339 ENDIF 433 340 ! … … 442 349 !! 443 350 INTEGER :: ji, jj 351 INTEGER :: istart, iend, jstart, jend 444 352 !!---------------------------------------------------------------------- 445 353 ! 446 354 IF( Agrif_Root() ) RETURN 447 355 ! 448 IF((nbondi == -1).OR.(nbondi == 2)) THEN 356 !--- West ---! 357 istart = 2 358 iend = nbghostcells+1 359 DO ji = mi0(istart), mi1(iend) 449 360 DO jj=1,jpj 450 va_e(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * hvr_e(2:nbghostcells+1,jj) 451 ! Specified fluxes: 452 ua_e(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * hur_e(2:nbghostcells+1,jj) 453 ! Characteristics method (only if ghostcells=1): 454 !alt ua_e(2,jj) = 0.5_wp * ( ubdy_w(jj) * hur_e(2,jj) + ua_e(3,jj) & 455 !alt & - sqrt(grav * hur_e(2,jj)) * (sshn_e(3,jj) - hbdy_w(jj)) ) 456 END DO 457 ENDIF 458 ! 459 IF((nbondi == 1).OR.(nbondi == 2)) THEN 361 va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj) 362 ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj) 363 END DO 364 END DO 365 ! 366 !--- East ---! 367 istart = jpiglo-nbghostcells 368 iend = jpiglo-1 369 DO ji = mi0(istart), mi1(iend) 460 370 DO jj=1,jpj 461 va_e(nlci-nbghostcells:nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * hvr_e(nlci-nbghostcells:nlci-1,jj) 462 ! Specified fluxes: 463 ua_e(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * hur_e(nlci-nbghostcells-1:nlci-2,jj) 464 ! Characteristics method (only if ghostcells=1): 465 !alt ua_e(nlci-2,jj) = 0.5_wp * ( ubdy_e(jj) * hur_e(nlci-2,jj) + ua_e(nlci-3,jj) & 466 !alt & + sqrt(grav * hur_e(nlci-2,jj)) * (sshn_e(nlci-2,jj) - hbdy_e(jj)) ) 467 END DO 468 ENDIF 469 ! 470 IF((nbondj == -1).OR.(nbondj == 2)) THEN 371 va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj) 372 END DO 373 END DO 374 istart = jpiglo-nbghostcells-1 375 iend = jpiglo-2 376 DO ji = mi0(istart), mi1(iend) 377 DO jj=1,jpj 378 ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj) 379 END DO 380 END DO 381 ! 382 !--- South ---! 383 jstart = 2 384 jend = nbghostcells+1 385 DO jj = mj0(jstart), mj1(jend) 471 386 DO ji=1,jpi 472 ua_e(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * hur_e(ji,2:nbghostcells+1) 473 ! Specified fluxes: 474 va_e(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * hvr_e(ji,2:nbghostcells+1) 475 ! Characteristics method (only if ghostcells=1): 476 !alt va_e(ji,2) = 0.5_wp * ( vbdy_s(ji) * hvr_e(ji,2) + va_e(ji,3) & 477 !alt & - sqrt(grav * hvr_e(ji,2)) * (sshn_e(ji,3) - hbdy_s(ji)) ) 478 END DO 479 ENDIF 480 ! 481 IF((nbondj == 1).OR.(nbondj == 2)) THEN 387 ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj) 388 va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj) 389 END DO 390 END DO 391 ! 392 !--- North ---! 393 jstart = jpjglo-nbghostcells 394 jend = jpjglo-1 395 DO jj = mj0(jstart), mj1(jend) 482 396 DO ji=1,jpi 483 ua_e(ji,nlcj-nbghostcells:nlcj-1) = ubdy_n(ji,1:nbghostcells) * hur_e(ji,nlcj-nbghostcells:nlcj-1) 484 ! Specified fluxes: 485 va_e(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * hvr_e(ji,nlcj-nbghostcells-1:nlcj-2) 486 ! Characteristics method (only if ghostcells=1): 487 !alt va_e(ji,nlcj-2) = 0.5_wp * ( vbdy_n(ji) * hvr_e(ji,nlcj-2) + va_e(ji,nlcj-3) & 488 !alt & + sqrt(grav * hvr_e(ji,nlcj-2)) * (sshn_e(ji,nlcj-2) - hbdy_n(ji)) ) 489 END DO 490 ENDIF 397 ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj) 398 END DO 399 END DO 400 jstart = jpjglo-nbghostcells-1 401 jend = jpjglo-2 402 DO jj = mj0(jstart), mj1(jend) 403 DO ji=1,jpi 404 va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj) 405 END DO 406 END DO 491 407 ! 492 408 END SUBROUTINE Agrif_dyn_ts 493 409 410 SUBROUTINE Agrif_dyn_ts_flux( jn, zu, zv ) 411 !!---------------------------------------------------------------------- 412 !! *** ROUTINE Agrif_dyn_ts_flux *** 413 !!---------------------------------------------------------------------- 414 INTEGER, INTENT(in) :: jn 415 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zu, zv 416 !! 417 INTEGER :: ji, jj 418 INTEGER :: istart, iend, jstart, jend 419 !!---------------------------------------------------------------------- 420 ! 421 IF( Agrif_Root() ) RETURN 422 ! 423 !--- West ---! 424 istart = 2 425 iend = nbghostcells+1 426 DO ji = mi0(istart), mi1(iend) 427 DO jj=1,jpj 428 zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj) 429 zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj) 430 END DO 431 END DO 432 ! 433 !--- East ---! 434 istart = jpiglo-nbghostcells 435 iend = jpiglo-1 436 DO ji = mi0(istart), mi1(iend) 437 DO jj=1,jpj 438 zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj) 439 END DO 440 END DO 441 istart = jpiglo-nbghostcells-1 442 iend = jpiglo-2 443 DO ji = mi0(istart), mi1(iend) 444 DO jj=1,jpj 445 zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj) 446 END DO 447 END DO 448 ! 449 !--- South ---! 450 jstart = 2 451 jend = nbghostcells+1 452 DO jj = mj0(jstart), mj1(jend) 453 DO ji=1,jpi 454 zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj) 455 zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj) 456 END DO 457 END DO 458 ! 459 !--- North ---! 460 jstart = jpjglo-nbghostcells 461 jend = jpjglo-1 462 DO jj = mj0(jstart), mj1(jend) 463 DO ji=1,jpi 464 zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj) 465 END DO 466 END DO 467 jstart = jpjglo-nbghostcells-1 468 jend = jpjglo-2 469 DO jj = mj0(jstart), mj1(jend) 470 DO ji=1,jpi 471 zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj) 472 END DO 473 END DO 474 ! 475 END SUBROUTINE Agrif_dyn_ts_flux 494 476 495 477 SUBROUTINE Agrif_dta_ts( kt ) … … 511 493 ! 512 494 ! Interpolate barotropic fluxes 513 Agrif_SpecialValue =0._wp495 Agrif_SpecialValue = 0._wp 514 496 Agrif_UseSpecialValue = ln_spc_dyn 497 ! 498 ! Set bdy time interpolation stage to 0 (latter incremented locally do deal with corners) 499 utint_stage(:,:) = 0 500 vtint_stage(:,:) = 0 515 501 ! 516 502 IF( ll_int_cons ) THEN ! Conservative interpolation 517 503 ! order matters here !!!!!! 518 504 CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b ) ! Time integrated 519 CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b ) 505 CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b ) 506 ! 520 507 bdy_tinterp = 1 521 508 CALL Agrif_Bc_variable( unb_id , calledweight=1._wp, procname=interpunb ) ! After 522 CALL Agrif_Bc_variable( vnb_id , calledweight=1._wp, procname=interpvnb ) 509 CALL Agrif_Bc_variable( vnb_id , calledweight=1._wp, procname=interpvnb ) 510 ! 523 511 bdy_tinterp = 2 524 512 CALL Agrif_Bc_variable( unb_id , calledweight=0._wp, procname=interpunb ) ! Before 525 CALL Agrif_Bc_variable( vnb_id , calledweight=0._wp, procname=interpvnb ) 513 CALL Agrif_Bc_variable( vnb_id , calledweight=0._wp, procname=interpvnb ) 526 514 ELSE ! Linear interpolation 527 bdy_tinterp = 0 528 ubdy_w(:,:) = 0._wp ; vbdy_w(:,:) = 0._wp 529 ubdy_e(:,:) = 0._wp ; vbdy_e(:,:) = 0._wp 530 ubdy_n(:,:) = 0._wp ; vbdy_n(:,:) = 0._wp 531 ubdy_s(:,:) = 0._wp ; vbdy_s(:,:) = 0._wp 515 ! 516 ubdy(:,:) = 0._wp ; vbdy(:,:) = 0._wp 532 517 CALL Agrif_Bc_variable( unb_id, procname=interpunb ) 533 518 CALL Agrif_Bc_variable( vnb_id, procname=interpvnb ) … … 544 529 INTEGER, INTENT(in) :: kt 545 530 ! 546 INTEGER :: ji, jj, indx, indy 531 INTEGER :: ji, jj 532 INTEGER :: istart, iend, jstart, jend 547 533 !!---------------------------------------------------------------------- 548 534 ! … … 557 543 ! 558 544 ! --- West --- ! 559 IF((nbondi == -1).OR.(nbondi == 2)) THEN 560 indx = 1+nbghostcells 545 istart = 2 546 iend = 1 + nbghostcells 547 DO ji = mi0(istart), mi1(iend) 561 548 DO jj = 1, jpj 562 DO ji = 2, indx 563 ssh(ji,jj,Krhs_a) = hbdy_w(ji-1,jj) 564 ENDDO 549 ssh(ji,jj,Krhs_a) = hbdy(ji,jj) 565 550 ENDDO 566 END IF551 ENDDO 567 552 ! 568 553 ! --- East --- ! 569 IF((nbondi == 1).OR.(nbondi == 2)) THEN 570 indx = nlci-nbghostcells 554 istart = jpiglo - nbghostcells 555 iend = jpiglo - 1 556 DO ji = mi0(istart), mi1(iend) 571 557 DO jj = 1, jpj 572 DO ji = indx, nlci-1 573 ssh(ji,jj,Krhs_a) = hbdy_e(ji-indx+1,jj) 574 ENDDO 558 ssh(ji,jj,Krhs_a) = hbdy(ji,jj) 575 559 ENDDO 576 END IF560 ENDDO 577 561 ! 578 562 ! --- South --- ! 579 IF((nbondj == -1).OR.(nbondj == 2)) THEN 580 indy = 1+nbghostcells 581 DO jj = 2, indy 582 DO ji = 1, jpi 583 ssh(ji,jj,Krhs_a) = hbdy_s(ji,jj-1) 584 ENDDO 563 jstart = 2 564 jend = 1 + nbghostcells 565 DO jj = mj0(jstart), mj1(jend) 566 DO ji = 1, jpi 567 ssh(ji,jj,Krhs_a) = hbdy(ji,jj) 585 568 ENDDO 586 END IF569 ENDDO 587 570 ! 588 571 ! --- North --- ! 589 IF((nbondj == 1).OR.(nbondj == 2)) THEN 590 indy = nlcj-nbghostcells 591 DO jj = indy, nlcj-1 592 DO ji = 1, jpi 593 ssh(ji,jj,Krhs_a) = hbdy_n(ji,jj-indy+1) 594 ENDDO 572 jstart = jpjglo - nbghostcells 573 jend = jpjglo - 1 574 DO jj = mj0(jstart), mj1(jend) 575 DO ji = 1, jpi 576 ssh(ji,jj,Krhs_a) = hbdy(ji,jj) 595 577 ENDDO 596 END IF578 ENDDO 597 579 ! 598 580 END SUBROUTINE Agrif_ssh … … 605 587 INTEGER, INTENT(in) :: jn 606 588 !! 607 INTEGER :: ji, jj , indx, indy608 !!----------------------------------------------------------------------609 !! clem ghost (starting at i,j=1 is important I think otherwise you introduce a grad(ssh)/=0 at point 2)589 INTEGER :: ji, jj 590 INTEGER :: istart, iend, jstart, jend 591 !!---------------------------------------------------------------------- 610 592 ! 611 593 IF( Agrif_Root() ) RETURN 612 594 ! 613 595 ! --- West --- ! 614 IF((nbondi == -1).OR.(nbondi == 2)) THEN 615 indx = 1+nbghostcells 596 istart = 2 597 iend = 1+nbghostcells 598 DO ji = mi0(istart), mi1(iend) 616 599 DO jj = 1, jpj 617 DO ji = 2, indx 618 ssha_e(ji,jj) = hbdy_w(ji-1,jj) 619 ENDDO 600 ssha_e(ji,jj) = hbdy(ji,jj) 620 601 ENDDO 621 END IF602 ENDDO 622 603 ! 623 604 ! --- East --- ! 624 IF((nbondi == 1).OR.(nbondi == 2)) THEN 625 indx = nlci-nbghostcells 605 istart = jpiglo - nbghostcells 606 iend = jpiglo - 1 607 DO ji = mi0(istart), mi1(iend) 626 608 DO jj = 1, jpj 627 DO ji = indx, nlci-1 628 ssha_e(ji,jj) = hbdy_e(ji-indx+1,jj) 629 ENDDO 609 ssha_e(ji,jj) = hbdy(ji,jj) 630 610 ENDDO 631 END IF611 ENDDO 632 612 ! 633 613 ! --- South --- ! 634 IF((nbondj == -1).OR.(nbondj == 2)) THEN 635 indy = 1+nbghostcells 636 DO jj = 2, indy 637 DO ji = 1, jpi 638 ssha_e(ji,jj) = hbdy_s(ji,jj-1) 639 ENDDO 614 jstart = 2 615 jend = 1+nbghostcells 616 DO jj = mj0(jstart), mj1(jend) 617 DO ji = 1, jpi 618 ssha_e(ji,jj) = hbdy(ji,jj) 640 619 ENDDO 641 END IF620 ENDDO 642 621 ! 643 622 ! --- North --- ! 644 IF((nbondj == 1).OR.(nbondj == 2)) THEN 645 indy = nlcj-nbghostcells 646 DO jj = indy, nlcj-1 647 DO ji = 1, jpi 648 ssha_e(ji,jj) = hbdy_n(ji,jj-indy+1) 649 ENDDO 623 jstart = jpjglo - nbghostcells 624 jend = jpjglo - 1 625 DO jj = mj0(jstart), mj1(jend) 626 DO ji = 1, jpi 627 ssha_e(ji,jj) = hbdy(ji,jj) 650 628 ENDDO 651 END IF629 ENDDO 652 630 ! 653 631 END SUBROUTINE Agrif_ssh_ts … … 675 653 676 654 677 SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before , nb, ndir)655 SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before ) 678 656 !!---------------------------------------------------------------------- 679 657 !! *** ROUTINE interptsn *** … … 682 660 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2 683 661 LOGICAL , INTENT(in ) :: before 684 INTEGER , INTENT(in ) :: nb , ndir 685 ! 686 INTEGER :: ji, jj, jk, jn, iref, jref, ibdy, jbdy ! dummy loop indices 687 INTEGER :: imin, imax, jmin, jmax, N_in, N_out 688 REAL(wp) :: zrho, z1, z2, z3, z4, z5, z6, z7 689 LOGICAL :: western_side, eastern_side,northern_side,southern_side 662 ! 663 INTEGER :: ji, jj, jk, jn ! dummy loop indices 664 INTEGER :: N_in, N_out 690 665 ! vertical interpolation: 691 REAL(wp) , DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child692 REAL(wp), DIMENSION(k1:k2, n1:n2-1) :: tabin666 REAL(wp) :: zhtot 667 REAL(wp), DIMENSION(k1:k2,1:jpts) :: tabin 693 668 REAL(wp), DIMENSION(k1:k2) :: h_in 694 669 REAL(wp), DIMENSION(1:jpk) :: h_out 695 REAL(wp) :: h_diff670 !!---------------------------------------------------------------------- 696 671 697 672 IF( before ) THEN … … 707 682 708 683 # if defined key_vertical 684 ! Interpolate thicknesses 685 ! Warning: these are masked, hence extrapolated prior interpolation. 709 686 DO jk=k1,k2 710 687 DO jj=j1,j2 711 688 DO ji=i1,i2 712 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a)689 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 713 690 END DO 714 691 END DO 715 692 END DO 693 694 ! Extrapolate thicknesses in partial bottom cells: 695 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 696 IF (ln_zps) THEN 697 DO jj=j1,j2 698 DO ji=i1,i2 699 jk = mbkt(ji,jj) 700 ptab(ji,jj,jk,jpts+1) = 0._wp 701 END DO 702 END DO 703 END IF 704 705 ! Save ssh at last level: 706 IF (.NOT.ln_linssh) THEN 707 ptab(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 708 ELSE 709 ptab(i1:i2,j1:j2,k2,jpts+1) = 0._wp 710 END IF 716 711 # endif 717 712 ELSE 718 713 719 western_side = (nb == 1).AND.(ndir == 1) ; eastern_side = (nb == 1).AND.(ndir == 2) 720 southern_side = (nb == 2).AND.(ndir == 1) ; northern_side = (nb == 2).AND.(ndir == 2) 721 722 # if defined key_vertical 714 # if defined key_vertical 715 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp 716 723 717 DO jj=j1,j2 724 718 DO ji=i1,i2 725 iref = ji726 jref = jj727 if(western_side) iref=MAX(2 ,ji)728 if(eastern_side) iref=MIN(nlci-1,ji)729 if(southern_side) jref=MAX(2 ,jj)730 if(northern_side) jref=MIN(nlcj-1,jj)731 N_in = 0732 DO jk=k1,k2 !k2 = jpk of parent grid733 IF (ptab(ji,jj,jk,n2) == 0) EXIT734 N_in = N_in + 1719 ts(ji,jj,:,:,Krhs_a) = 0._wp 720 N_in = mbkt_parent(ji,jj) 721 zhtot = 0._wp 722 DO jk=1,N_in !k2 = jpk of parent grid 723 IF (jk==N_in) THEN 724 h_in(jk) = ht0_parent(ji,jj) + ptab(ji,jj,k2,n2) - zhtot 725 ELSE 726 h_in(jk) = ptab(ji,jj,jk,n2) 727 ENDIF 728 zhtot = zhtot + h_in(jk) 735 729 tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1) 736 h_in(N_in) = ptab(ji,jj,jk,n2)737 730 END DO 738 731 N_out = 0 739 732 DO jk=1,jpk ! jpk of child grid 740 IF (tmask( iref,jref,jk) == 0) EXIT733 IF (tmask(ji,jj,jk) == 0._wp) EXIT 741 734 N_out = N_out + 1 742 h_out(jk) = e3t( iref,jref,jk,Kmm_a)735 h_out(jk) = e3t(ji,jj,jk,Krhs_a) 743 736 ENDDO 744 IF (N_in > 0) THEN 745 DO jn=1,jpts 746 call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 747 ENDDO 737 IF (N_in*N_out > 0) THEN 738 CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),h_out(1:N_out),N_in,N_out,jpts) 748 739 ENDIF 749 740 ENDDO 750 741 ENDDO 751 742 # else 752 ptab_child(i1:i2,j1:j2,1:jpk,1:jpts) = ptab(i1:i2,j1:j2,1:jpk,1:jpts)753 # endif754 743 ! 755 744 DO jn=1, jpts 756 ts(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab_child(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 757 END DO 758 759 IF ( .NOT.lk_agrif_clp ) THEN 760 ! 761 imin = i1 ; imax = i2 762 jmin = j1 ; jmax = j2 763 ! 764 ! Remove CORNERS 765 IF((nbondj == -1).OR.(nbondj == 2)) jmin = 2 + nbghostcells 766 IF((nbondj == +1).OR.(nbondj == 2)) jmax = nlcj - nbghostcells - 1 767 IF((nbondi == -1).OR.(nbondi == 2)) imin = 2 + nbghostcells 768 IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci - nbghostcells - 1 769 ! 770 IF( eastern_side ) THEN 771 zrho = Agrif_Rhox() 772 z1 = ( zrho - 1._wp ) * 0.5_wp 773 z3 = ( zrho - 1._wp ) / ( zrho + 1._wp ) 774 z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp ) 775 z7 = - ( zrho - 1._wp ) / ( zrho + 3._wp ) 776 z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7 777 ! 778 ibdy = nlci-nbghostcells 779 DO jn = 1, jpts 780 ts(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs_a) = z1 * ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) & 781 & + z2 * ptab_child(ibdy ,jmin:jmax,1:jpkm1,jn) 782 DO jk = 1, jpkm1 783 DO jj = jmin,jmax 784 IF( umask(ibdy-1,jj,jk) == 0._wp ) THEN 785 ts(ibdy,jj,jk,jn,Krhs_a) = ts(ibdy+1,jj,jk,jn,Krhs_a) * tmask(ibdy,jj,jk) 786 ELSE 787 ts(ibdy,jj,jk,jn,Krhs_a) = ( z4 * ts(ibdy+1,jj,jk,jn,Krhs_a) & 788 & + z3 * ts(ibdy-1,jj,jk,jn,Krhs_a) & 789 & ) * tmask(ibdy ,jj,jk) 790 IF( uu(ibdy-1,jj,jk,Kmm_a) > 0._wp ) THEN 791 ts(ibdy,jj,jk,jn,Krhs_a) = ( z6 * ts(ibdy-1,jj,jk,jn,Krhs_a) & 792 & + z5 * ts(ibdy+1,jj,jk,jn,Krhs_a) & 793 & + z7 * ts(ibdy-2,jj,jk,jn,Krhs_a) & 794 & ) * tmask(ibdy ,jj,jk) 795 ENDIF 796 ENDIF 797 END DO 798 END DO 799 ! Restore ghost points: 800 ts(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs_a) = ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) & 801 & * tmask(ibdy+1,jmin:jmax,1:jpkm1) 802 END DO 803 ENDIF 804 ! 805 IF( northern_side ) THEN 806 zrho = Agrif_Rhoy() 807 z1 = ( zrho - 1._wp ) * 0.5_wp 808 z3 = ( zrho - 1._wp ) / ( zrho + 1._wp ) 809 z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp ) 810 z7 = - ( zrho - 1._wp ) / ( zrho + 3._wp ) 811 z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7 812 ! 813 jbdy = nlcj-nbghostcells 814 DO jn = 1, jpts 815 ts(imin:imax,jbdy+1,1:jpkm1,jn,Krhs_a) = z1 * ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) & 816 & + z2 * ptab_child(imin:imax,jbdy ,1:jpkm1,jn) 817 DO jk = 1, jpkm1 818 DO ji = imin,imax 819 IF( vmask(ji,jbdy-1,jk) == 0._wp ) THEN 820 ts(ji,jbdy,jk,jn,Krhs_a) = ts(ji,jbdy+1,jk,jn,Krhs_a) * tmask(ji,jbdy,jk) 821 ELSE 822 ts(ji,jbdy,jk,jn,Krhs_a)=( z4 * ts(ji,jbdy+1,jk,jn,Krhs_a) & 823 & + z3 * ts(ji,jbdy-1,jk,jn,Krhs_a) & 824 & ) * tmask(ji,jbdy ,jk) 825 IF (vv(ji,jbdy-1,jk,Kmm_a) > 0._wp ) THEN 826 ts(ji,jbdy,jk,jn,Krhs_a)=( z6 * ts(ji,jbdy-1,jk,jn,Krhs_a) & 827 & + z5 * ts(ji,jbdy+1,jk,jn,Krhs_a) & 828 & + z7 * ts(ji,jbdy-2,jk,jn,Krhs_a) & 829 & ) * tmask(ji,jbdy ,jk) 830 ENDIF 831 ENDIF 832 END DO 833 END DO 834 ! Restore ghost points: 835 ts(imin:imax,jbdy+1,1:jpkm1,jn,Krhs_a) = ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) & 836 & * tmask(imin:imax,jbdy+1,1:jpkm1) 837 END DO 838 ENDIF 839 ! 840 IF( western_side ) THEN 841 zrho = Agrif_Rhox() 842 z1 = ( zrho - 1._wp ) * 0.5_wp 843 z3 = ( zrho - 1._wp ) / ( zrho + 1._wp ) 844 z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp ) 845 z7 = - ( zrho - 1._wp ) / ( zrho + 3._wp ) 846 z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7 847 ! 848 ibdy = 1+nbghostcells 849 DO jn = 1, jpts 850 ts(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs_a) = z1 * ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) & 851 & + z2 * ptab_child(ibdy ,jmin:jmax,1:jpkm1,jn) 852 DO jk = 1, jpkm1 853 DO jj = jmin,jmax 854 IF( umask(ibdy,jj,jk) == 0._wp ) THEN 855 ts(ibdy,jj,jk,jn,Krhs_a) = ts(ibdy-1,jj,jk,jn,Krhs_a) * tmask(ibdy,jj,jk) 856 ELSE 857 ts(ibdy,jj,jk,jn,Krhs_a) = ( z4 * ts(ibdy-1,jj,jk,jn,Krhs_a) & 858 & + z3 * ts(ibdy+1,jj,jk,jn,Krhs_a) & 859 & ) * tmask(ibdy ,jj,jk) 860 IF( uu(ibdy,jj,jk,Kmm_a) < 0._wp ) THEN 861 ts(ibdy,jj,jk,jn,Krhs_a) = ( z6 * ts(ibdy+1,jj,jk,jn,Krhs_a) & 862 & + z5 * ts(ibdy-1,jj,jk,jn,Krhs_a) & 863 & + z7 * ts(ibdy+2,jj,jk,jn,Krhs_a) & 864 & ) * tmask(ibdy,jj,jk) 865 ENDIF 866 ENDIF 867 END DO 868 END DO 869 ! Restore ghost points: 870 ts(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs_a) = ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) & 871 & * tmask(ibdy-1,jmin:jmax,1:jpkm1) 872 END DO 873 ENDIF 874 ! 875 IF( southern_side ) THEN 876 zrho = Agrif_Rhoy() 877 z1 = ( zrho - 1._wp ) * 0.5_wp 878 z3 = ( zrho - 1._wp ) / ( zrho + 1._wp ) 879 z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp ) 880 z7 = - ( zrho - 1._wp ) / ( zrho + 3._wp ) 881 z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7 882 ! 883 jbdy=1+nbghostcells 884 DO jn = 1, jpts 885 ts(imin:imax,jbdy-1,1:jpkm1,jn,Krhs_a) = z1 * ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) & 886 & + z2 * ptab_child(imin:imax,jbdy ,1:jpkm1,jn) 887 DO jk = 1, jpkm1 888 DO ji = imin,imax 889 IF( vmask(ji,jbdy,jk) == 0._wp ) THEN 890 ts(ji,jbdy,jk,jn,Krhs_a) = ts(ji,jbdy-1,jk,jn,Krhs_a) * tmask(ji,jbdy,jk) 891 ELSE 892 ts(ji,jbdy,jk,jn,Krhs_a) = ( z4 * ts(ji,jbdy-1,jk,jn,Krhs_a) & 893 & + z3 * ts(ji,jbdy+1,jk,jn,Krhs_a) & 894 & ) * tmask(ji,jbdy ,jk) 895 IF( vv(ji,jbdy,jk,Kmm_a) < 0._wp ) THEN 896 ts(ji,jbdy,jk,jn,Krhs_a) = ( z6 * ts(ji,jbdy+1,jk,jn,Krhs_a) & 897 & + z5 * ts(ji,jbdy-1,jk,jn,Krhs_a) & 898 & + z7 * ts(ji,jbdy+2,jk,jn,Krhs_a) & 899 & ) * tmask(ji,jbdy ,jk) 900 ENDIF 901 ENDIF 902 END DO 903 END DO 904 ! Restore ghost points: 905 ts(imin:imax,jbdy-1,1:jpkm1,jn,Krhs_a) = ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) & 906 & * tmask(imin:imax,jbdy-1,1:jpkm1) 907 END DO 908 ENDIF 909 ! 910 ENDIF 745 ts(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 746 END DO 747 # endif 748 911 749 ENDIF 912 750 ! 913 751 END SUBROUTINE interptsn 914 752 915 SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before , nb, ndir)753 SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before ) 916 754 !!---------------------------------------------------------------------- 917 755 !! *** ROUTINE interpsshn *** … … 920 758 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 921 759 LOGICAL , INTENT(in ) :: before 922 INTEGER , INTENT(in ) :: nb , ndir 923 ! 924 LOGICAL :: western_side, eastern_side,northern_side,southern_side 760 ! 925 761 !!---------------------------------------------------------------------- 926 762 ! … … 928 764 ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kmm_a) 929 765 ELSE 930 western_side = (nb == 1).AND.(ndir == 1) 931 eastern_side = (nb == 1).AND.(ndir == 2) 932 southern_side = (nb == 2).AND.(ndir == 1) 933 northern_side = (nb == 2).AND.(ndir == 2) 934 !! clem ghost 935 IF(western_side) hbdy_w(1:nbghostcells,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1) 936 IF(eastern_side) hbdy_e(1:nbghostcells,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1) 937 IF(southern_side) hbdy_s(i1:i2,1:nbghostcells) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1) 938 IF(northern_side) hbdy_n(i1:i2,1:nbghostcells) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1) 766 hbdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1) 939 767 ENDIF 940 768 ! 941 769 END SUBROUTINE interpsshn 942 770 943 SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before , nb, ndir)771 SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before ) 944 772 !!---------------------------------------------------------------------- 945 773 !! *** ROUTINE interpun *** … … 949 777 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 950 778 LOGICAL, INTENT(in) :: before 951 INTEGER, INTENT(in) :: nb , ndir952 779 !! 953 780 INTEGER :: ji,jj,jk 954 REAL(wp) :: zrhoy 781 REAL(wp) :: zrhoy, zhtot 955 782 ! vertical interpolation: 956 783 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 957 784 REAL(wp), DIMENSION(1:jpk) :: h_out 958 INTEGER :: N_in, N_out , iref785 INTEGER :: N_in, N_out 959 786 REAL(wp) :: h_diff 960 LOGICAL :: western_side, eastern_side961 787 !!--------------------------------------------- 962 788 ! … … 965 791 DO jj=j1,j2 966 792 DO ji=i1,i2 967 ptab(ji,jj,jk,1) = ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) & 968 & * uu(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) ) 793 ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)*umask(ji,jj,jk)) 969 794 # if defined key_vertical 970 ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)) 795 ! Interpolate thicknesses (masked for subsequent extrapolation) 796 ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 971 797 # endif 972 798 END DO 973 799 END DO 974 800 END DO 801 # if defined key_vertical 802 ! Extrapolate thicknesses in partial bottom cells: 803 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 804 IF (ln_zps) THEN 805 DO jj=j1,j2 806 DO ji=i1,i2 807 jk = mbku(ji,jj) 808 ptab(ji,jj,jk,2) = 0._wp 809 END DO 810 END DO 811 END IF 812 ! Save ssh at last level: 813 ptab(i1:i2,j1:j2,k2,2) = 0._wp 814 IF (.NOT.ln_linssh) THEN 815 ! This vertical sum below should be replaced by the sea-level at U-points (optimization): 816 DO jk=1,jpk 817 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3u(i1:i2,j1:j2,jk,Kmm_a) * umask(i1:i2,j1:j2,jk) 818 END DO 819 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hu_0(i1:i2,j1:j2) 820 END IF 821 # endif 822 ! 975 823 ELSE 976 824 zrhoy = Agrif_rhoy() 977 825 # if defined key_vertical 978 826 ! VERTICAL REFINEMENT BEGIN 979 western_side = (nb == 1).AND.(ndir == 1) 980 eastern_side = (nb == 1).AND.(ndir == 2)827 828 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 981 829 982 830 DO ji=i1,i2 983 iref = ji984 IF (western_side) iref = MAX(2,ji)985 IF (eastern_side) iref = MIN(nlci-2,ji)986 831 DO jj=j1,j2 987 N_in = 0 988 DO jk=k1,k2 989 IF (ptab(ji,jj,jk,2) == 0) EXIT 990 N_in = N_in + 1 991 tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2) 992 h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 832 uu(ji,jj,:,Krhs_a) = 0._wp 833 N_in = mbku_parent(ji,jj) 834 zhtot = 0._wp 835 DO jk=1,N_in 836 IF (jk==N_in) THEN 837 h_in(jk) = hu0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 838 ELSE 839 h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 840 ENDIF 841 zhtot = zhtot + h_in(jk) 842 tabin(jk) = ptab(ji,jj,jk,1)/(e2u(ji,jj)*zrhoy*h_in(jk)) 993 843 ENDDO 994 995 IF (N_in == 0) THEN 996 uu(ji,jj,:,Krhs_a) = 0._wp 997 CYCLE 998 ENDIF 999 844 1000 845 N_out = 0 1001 846 DO jk=1,jpk 1002 if (umask( iref,jj,jk) == 0) EXIT847 if (umask(ji,jj,jk) == 0) EXIT 1003 848 N_out = N_out + 1 1004 h_out(N_out) = e3u( iref,jj,jk,Krhs_a)849 h_out(N_out) = e3u(ji,jj,jk,Krhs_a) 1005 850 ENDDO 1006 1007 IF (N_out == 0) THEN 1008 uu(ji,jj,:,Krhs_a) = 0._wp 1009 CYCLE 851 IF (N_in*N_out > 0) THEN 852 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),uu(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1) 1010 853 ENDIF 1011 1012 IF (N_in * N_out > 0) THEN1013 h_diff = SUM( h_out(1:N_out) ) - SUM( h_in(1:N_in) )1014 ! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly1015 if (h_diff < -1.e4) then1016 print *,'CHECK YOUR BATHY ...', h_diff, SUM( h_out(1:N_out) ), SUM( h_in(1:N_in) )1017 ! stop1018 endif1019 ENDIF1020 call reconstructandremap( tabin(1:N_in) , h_in(1:N_in), uu(ji,jj,1:N_out,Krhs_a), &1021 & h_out(1:N_out), N_in , N_out )1022 854 ENDDO 1023 855 ENDDO … … 1035 867 END SUBROUTINE interpun 1036 868 1037 SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before , nb, ndir)869 SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before ) 1038 870 !!---------------------------------------------------------------------- 1039 871 !! *** ROUTINE interpvn *** … … 1043 875 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 1044 876 LOGICAL, INTENT(in) :: before 1045 INTEGER, INTENT(in) :: nb , ndir1046 877 ! 1047 878 INTEGER :: ji,jj,jk … … 1050 881 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 1051 882 REAL(wp), DIMENSION(1:jpk) :: h_out 1052 INTEGER :: N_in, N_out, jref 1053 REAL(wp) :: h_diff 1054 LOGICAL :: northern_side,southern_side 883 INTEGER :: N_in, N_out 884 REAL(wp) :: h_diff, zhtot 1055 885 !!--------------------------------------------- 1056 886 ! … … 1059 889 DO jj=j1,j2 1060 890 DO ji=i1,i2 1061 ptab(ji,jj,jk,1) = ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) & 1062 & * vv(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) ) 891 ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)*vmask(ji,jj,jk)) 1063 892 # if defined key_vertical 893 ! Interpolate thicknesses (masked for subsequent extrapolation) 1064 894 ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) 1065 895 # endif … … 1067 897 END DO 1068 898 END DO 899 # if defined key_vertical 900 ! Extrapolate thicknesses in partial bottom cells: 901 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 902 IF (ln_zps) THEN 903 DO jj=j1,j2 904 DO ji=i1,i2 905 jk = mbkv(ji,jj) 906 ptab(ji,jj,jk,2) = 0._wp 907 END DO 908 END DO 909 END IF 910 ! Save ssh at last level: 911 ptab(i1:i2,j1:j2,k2,2) = 0._wp 912 IF (.NOT.ln_linssh) THEN 913 ! This vertical sum below should be replaced by the sea-level at V-points (optimization): 914 DO jk=1,jpk 915 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3v(i1:i2,j1:j2,jk,Kmm_a) * vmask(i1:i2,j1:j2,jk) 916 END DO 917 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hv_0(i1:i2,j1:j2) 918 END IF 919 # endif 1069 920 ELSE 1070 921 zrhox = Agrif_rhox() 1071 922 # if defined key_vertical 1072 923 1073 southern_side = (nb == 2).AND.(ndir == 1) 1074 northern_side = (nb == 2).AND.(ndir == 2) 924 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 1075 925 1076 926 DO jj=j1,j2 1077 jref = jj1078 IF (southern_side) jref = MAX(2,jj)1079 IF (northern_side) jref = MIN(nlcj-2,jj)1080 927 DO ji=i1,i2 1081 N_in = 0 1082 DO jk=k1,k2 1083 if (ptab(ji,jj,jk,2) == 0) EXIT 1084 N_in = N_in + 1 1085 tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2) 1086 h_in(N_in) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 1087 END DO 1088 IF (N_in == 0) THEN 1089 vv(ji,jj,:,Krhs_a) = 0._wp 1090 CYCLE 1091 ENDIF 928 vv(ji,jj,:,Krhs_a) = 0._wp 929 N_in = mbkv_parent(ji,jj) 930 zhtot = 0._wp 931 DO jk=1,N_in 932 IF (jk==N_in) THEN 933 h_in(jk) = hv0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 934 ELSE 935 h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 936 ENDIF 937 zhtot = zhtot + h_in(jk) 938 tabin(jk) = ptab(ji,jj,jk,1)/(e1v(ji,jj)*zrhox*h_in(jk)) 939 ENDDO 1092 940 1093 941 N_out = 0 1094 942 DO jk=1,jpk 1095 if (vmask(ji,j ref,jk) == 0) EXIT943 if (vmask(ji,jj,jk) == 0) EXIT 1096 944 N_out = N_out + 1 1097 h_out(N_out) = e3v(ji,jref,jk,Krhs_a) 1098 END DO 1099 IF (N_out == 0) THEN 1100 vv(ji,jj,:,Krhs_a) = 0._wp 1101 CYCLE 945 h_out(N_out) = e3v(ji,jj,jk,Krhs_a) 946 END DO 947 IF (N_in*N_out > 0) THEN 948 call reconstructandremap(tabin(1:N_in),h_in(1:N_in),vv(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1) 1102 949 ENDIF 1103 call reconstructandremap( tabin(1:N_in) , h_in(1:N_in), vv(ji,jj,1:N_out,Krhs_a), &1104 & h_out(1:N_out), N_in , N_out )1105 950 END DO 1106 951 END DO … … 1114 959 END SUBROUTINE interpvn 1115 960 1116 SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before , nb, ndir)961 SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before) 1117 962 !!---------------------------------------------------------------------- 1118 963 !! *** ROUTINE interpunb *** … … 1121 966 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 1122 967 LOGICAL , INTENT(in ) :: before 1123 INTEGER , INTENT(in ) :: nb , ndir1124 968 ! 1125 969 INTEGER :: ji, jj 1126 970 REAL(wp) :: zrhoy, zrhot, zt0, zt1, ztcoeff 1127 LOGICAL :: western_side, eastern_side,northern_side,southern_side1128 971 !!---------------------------------------------------------------------- 1129 972 ! … … 1131 974 ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu(i1:i2,j1:j2,Kmm_a) * uu_b(i1:i2,j1:j2,Kmm_a) 1132 975 ELSE 1133 western_side = (nb == 1).AND.(ndir == 1)1134 eastern_side = (nb == 1).AND.(ndir == 2)1135 southern_side = (nb == 2).AND.(ndir == 1)1136 northern_side = (nb == 2).AND.(ndir == 2)1137 976 zrhoy = Agrif_Rhoy() 1138 977 zrhot = Agrif_rhot() … … 1140 979 zt0 = REAL(Agrif_NbStepint() , wp) / zrhot 1141 980 zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot 1142 ! Polynomial interpolation coefficients: 1143 IF( bdy_tinterp == 1 ) THEN 1144 ztcoeff = zrhot * ( zt1**2._wp * ( zt1 - 1._wp) & 1145 & - zt0**2._wp * ( zt0 - 1._wp) ) 1146 ELSEIF( bdy_tinterp == 2 ) THEN 1147 ztcoeff = zrhot * ( zt1 * ( zt1 - 1._wp)**2._wp & 1148 & - zt0 * ( zt0 - 1._wp)**2._wp ) 1149 ELSE 1150 ztcoeff = 1 1151 ENDIF 1152 ! 1153 IF(western_side) ubdy_w(1:nbghostcells,j1:j2) = ubdy_w(1:nbghostcells,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2) 1154 IF(eastern_side) ubdy_e(1:nbghostcells,j1:j2) = ubdy_e(1:nbghostcells,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2) 1155 IF(southern_side) ubdy_s(i1:i2,1:nbghostcells) = ubdy_s(i1:i2,1:nbghostcells) + ztcoeff * ptab(i1:i2,j1:j2) 1156 IF(northern_side) ubdy_n(i1:i2,1:nbghostcells) = ubdy_n(i1:i2,1:nbghostcells) + ztcoeff * ptab(i1:i2,j1:j2) 1157 ! 1158 IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN 1159 IF(western_side) ubdy_w(1:nbghostcells,j1:j2) = ubdy_w(1:nbghostcells,j1:j2) / (zrhoy*e2u(i1:i2,j1:j2)) * umask(i1:i2,j1:j2,1) 1160 IF(eastern_side) ubdy_e(1:nbghostcells,j1:j2) = ubdy_e(1:nbghostcells,j1:j2) / (zrhoy*e2u(i1:i2,j1:j2)) * umask(i1:i2,j1:j2,1) 1161 IF(southern_side) ubdy_s(i1:i2,1:nbghostcells) = ubdy_s(i1:i2,1:nbghostcells) / (zrhoy*e2u(i1:i2,j1:j2)) * umask(i1:i2,j1:j2,1) 1162 IF(northern_side) ubdy_n(i1:i2,1:nbghostcells) = ubdy_n(i1:i2,1:nbghostcells) / (zrhoy*e2u(i1:i2,j1:j2)) * umask(i1:i2,j1:j2,1) 1163 ENDIF 1164 ENDIF 981 ! 982 DO ji = i1, i2 983 DO jj = j1, j2 984 IF ( utint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN 985 IF ( utint_stage(ji,jj) == 1 ) THEN 986 ztcoeff = zrhot * ( zt1**2._wp * ( zt1 - 1._wp) & 987 & - zt0**2._wp * ( zt0 - 1._wp) ) 988 ELSEIF( utint_stage(ji,jj) == 2 ) THEN 989 ztcoeff = zrhot * ( zt1 * ( zt1 - 1._wp)**2._wp & 990 & - zt0 * ( zt0 - 1._wp)**2._wp ) 991 ELSEIF( utint_stage(ji,jj) == 0 ) THEN 992 ztcoeff = 1._wp 993 ELSE 994 ztcoeff = 0._wp 995 ENDIF 996 ! 997 ubdy(ji,jj) = ubdy(ji,jj) + ztcoeff * ptab(ji,jj) 998 ! 999 IF (( utint_stage(ji,jj) == 2 ).OR.( utint_stage(ji,jj) == 0 )) THEN 1000 ubdy(ji,jj) = ubdy(ji,jj) / (zrhoy*e2u(ji,jj)) * umask(ji,jj,1) 1001 ENDIF 1002 ! 1003 utint_stage(ji,jj) = utint_stage(ji,jj) + 1 1004 ENDIF 1005 END DO 1006 END DO 1007 END IF 1165 1008 ! 1166 1009 END SUBROUTINE interpunb 1167 1010 1168 1011 1169 SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before , nb, ndir)1012 SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before ) 1170 1013 !!---------------------------------------------------------------------- 1171 1014 !! *** ROUTINE interpvnb *** … … 1174 1017 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 1175 1018 LOGICAL , INTENT(in ) :: before 1176 INTEGER , INTENT(in ) :: nb , ndir 1177 ! 1178 INTEGER :: ji,jj 1019 ! 1020 INTEGER :: ji, jj 1179 1021 REAL(wp) :: zrhox, zrhot, zt0, zt1, ztcoeff 1180 LOGICAL :: western_side, eastern_side,northern_side,southern_side1181 1022 !!---------------------------------------------------------------------- 1182 1023 ! … … 1184 1025 ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv(i1:i2,j1:j2,Kmm_a) * vv_b(i1:i2,j1:j2,Kmm_a) 1185 1026 ELSE 1186 western_side = (nb == 1).AND.(ndir == 1)1187 eastern_side = (nb == 1).AND.(ndir == 2)1188 southern_side = (nb == 2).AND.(ndir == 1)1189 northern_side = (nb == 2).AND.(ndir == 2)1190 1027 zrhox = Agrif_Rhox() 1191 1028 zrhot = Agrif_rhot() 1192 1029 ! Time indexes bounds for integration 1193 1030 zt0 = REAL(Agrif_NbStepint() , wp) / zrhot 1194 zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot 1195 IF( bdy_tinterp == 1 ) THEN 1196 ztcoeff = zrhot * ( zt1**2._wp * ( zt1 - 1._wp) & 1197 & - zt0**2._wp * ( zt0 - 1._wp) ) 1198 ELSEIF( bdy_tinterp == 2 ) THEN 1199 ztcoeff = zrhot * ( zt1 * ( zt1 - 1._wp)**2._wp & 1200 & - zt0 * ( zt0 - 1._wp)**2._wp ) 1201 ELSE 1202 ztcoeff = 1 1203 ENDIF 1204 !! clem ghost 1205 IF(western_side) vbdy_w(1:nbghostcells,j1:j2) = vbdy_w(1:nbghostcells,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2) 1206 IF(eastern_side) vbdy_e(1:nbghostcells,j1:j2) = vbdy_e(1:nbghostcells,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2) 1207 IF(southern_side) vbdy_s(i1:i2,1:nbghostcells) = vbdy_s(i1:i2,1:nbghostcells) + ztcoeff * ptab(i1:i2,j1:j2) 1208 IF(northern_side) vbdy_n(i1:i2,1:nbghostcells) = vbdy_n(i1:i2,1:nbghostcells) + ztcoeff * ptab(i1:i2,j1:j2) 1209 ! 1210 IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN 1211 IF(western_side) vbdy_w(1:nbghostcells,j1:j2) = vbdy_w(1:nbghostcells,j1:j2) / (zrhox*e1v(i1:i2,j1:j2)) * vmask(i1:i2,j1:j2,1) 1212 IF(eastern_side) vbdy_e(1:nbghostcells,j1:j2) = vbdy_e(1:nbghostcells,j1:j2) / (zrhox*e1v(i1:i2,j1:j2)) * vmask(i1:i2,j1:j2,1) 1213 IF(southern_side) vbdy_s(i1:i2,1:nbghostcells) = vbdy_s(i1:i2,1:nbghostcells) / (zrhox*e1v(i1:i2,j1:j2)) * vmask(i1:i2,j1:j2,1) 1214 IF(northern_side) vbdy_n(i1:i2,1:nbghostcells) = vbdy_n(i1:i2,1:nbghostcells) / (zrhox*e1v(i1:i2,j1:j2)) * vmask(i1:i2,j1:j2,1) 1215 ENDIF 1031 zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot 1032 ! 1033 DO ji = i1, i2 1034 DO jj = j1, j2 1035 IF ( vtint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN 1036 IF ( vtint_stage(ji,jj) == 1 ) THEN 1037 ztcoeff = zrhot * ( zt1**2._wp * ( zt1 - 1._wp) & 1038 & - zt0**2._wp * ( zt0 - 1._wp) ) 1039 ELSEIF( vtint_stage(ji,jj) == 2 ) THEN 1040 ztcoeff = zrhot * ( zt1 * ( zt1 - 1._wp)**2._wp & 1041 & - zt0 * ( zt0 - 1._wp)**2._wp ) 1042 ELSEIF( vtint_stage(ji,jj) == 0 ) THEN 1043 ztcoeff = 1._wp 1044 ELSE 1045 ztcoeff = 0._wp 1046 ENDIF 1047 ! 1048 vbdy(ji,jj) = vbdy(ji,jj) + ztcoeff * ptab(ji,jj) 1049 ! 1050 IF (( vtint_stage(ji,jj) == 2 ).OR.( vtint_stage(ji,jj) == 0 )) THEN 1051 vbdy(ji,jj) = vbdy(ji,jj) / (zrhox*e1v(ji,jj)) * vmask(ji,jj,1) 1052 ENDIF 1053 ! 1054 vtint_stage(ji,jj) = vtint_stage(ji,jj) + 1 1055 ENDIF 1056 END DO 1057 END DO 1216 1058 ENDIF 1217 1059 ! … … 1219 1061 1220 1062 1221 SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before , nb, ndir)1063 SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before ) 1222 1064 !!---------------------------------------------------------------------- 1223 1065 !! *** ROUTINE interpub2b *** … … 1226 1068 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 1227 1069 LOGICAL , INTENT(in ) :: before 1228 INTEGER , INTENT(in ) :: nb , ndir1229 1070 ! 1230 1071 INTEGER :: ji,jj 1231 REAL(wp) :: zrhot, zt0, zt1,zat 1232 LOGICAL :: western_side, eastern_side,northern_side,southern_side 1072 REAL(wp) :: zrhot, zt0, zt1, zat 1233 1073 !!---------------------------------------------------------------------- 1234 1074 IF( before ) THEN … … 1239 1079 ENDIF 1240 1080 ELSE 1241 western_side = (nb == 1).AND.(ndir == 1)1242 eastern_side = (nb == 1).AND.(ndir == 2)1243 southern_side = (nb == 2).AND.(ndir == 1)1244 northern_side = (nb == 2).AND.(ndir == 2)1245 zrhot = Agrif_rhot()1246 ! Time indexes bounds for integration1247 zt0 = REAL(Agrif_NbStepint() , wp) / zrhot1248 zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot1249 ! Polynomial interpolation coefficients:1250 zat = zrhot * ( zt1**2._wp * (-2._wp*zt1 + 3._wp) &1251 & - zt0**2._wp * (-2._wp*zt0 + 3._wp) )1252 !! clem ghost1253 IF(western_side ) ubdy_w(1:nbghostcells,j1:j2) = zat * ptab(i1:i2,j1:j2)1254 IF(eastern_side ) ubdy_e(1:nbghostcells,j1:j2) = zat * ptab(i1:i2,j1:j2)1255 IF(southern_side) ubdy_s(i1:i2,1:nbghostcells) = zat * ptab(i1:i2,j1:j2)1256 IF(northern_side) ubdy_n(i1:i2,1:nbghostcells) = zat * ptab(i1:i2,j1:j2)1257 ENDIF1258 !1259 END SUBROUTINE interpub2b1260 1261 1262 SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before, nb, ndir )1263 !!----------------------------------------------------------------------1264 !! *** ROUTINE interpvb2b ***1265 !!----------------------------------------------------------------------1266 INTEGER , INTENT(in ) :: i1, i2, j1, j21267 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab1268 LOGICAL , INTENT(in ) :: before1269 INTEGER , INTENT(in ) :: nb , ndir1270 !1271 INTEGER :: ji,jj1272 REAL(wp) :: zrhot, zt0, zt1,zat1273 LOGICAL :: western_side, eastern_side,northern_side,southern_side1274 !!----------------------------------------------------------------------1275 !1276 IF( before ) THEN1277 IF ( ln_bt_fw ) THEN1278 ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)1279 ELSE1280 ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2)1281 ENDIF1282 ELSE1283 western_side = (nb == 1).AND.(ndir == 1)1284 eastern_side = (nb == 1).AND.(ndir == 2)1285 southern_side = (nb == 2).AND.(ndir == 1)1286 northern_side = (nb == 2).AND.(ndir == 2)1287 1081 zrhot = Agrif_rhot() 1288 1082 ! Time indexes bounds for integration … … 1293 1087 & - zt0**2._wp * (-2._wp*zt0 + 3._wp) ) 1294 1088 ! 1295 IF(western_side ) vbdy_w(1:nbghostcells,j1:j2) = zat * ptab(i1:i2,j1:j2) 1296 IF(eastern_side ) vbdy_e(1:nbghostcells,j1:j2) = zat * ptab(i1:i2,j1:j2) 1297 IF(southern_side) vbdy_s(i1:i2,1:nbghostcells) = zat * ptab(i1:i2,j1:j2) 1298 IF(northern_side) vbdy_n(i1:i2,1:nbghostcells) = zat * ptab(i1:i2,j1:j2) 1089 ubdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2) 1090 ! 1091 ! Update interpolation stage: 1092 utint_stage(i1:i2,j1:j2) = 1 1093 ENDIF 1094 ! 1095 END SUBROUTINE interpub2b 1096 1097 1098 SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before ) 1099 !!---------------------------------------------------------------------- 1100 !! *** ROUTINE interpvb2b *** 1101 !!---------------------------------------------------------------------- 1102 INTEGER , INTENT(in ) :: i1, i2, j1, j2 1103 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 1104 LOGICAL , INTENT(in ) :: before 1105 ! 1106 INTEGER :: ji,jj 1107 REAL(wp) :: zrhot, zt0, zt1, zat 1108 !!---------------------------------------------------------------------- 1109 ! 1110 IF( before ) THEN 1111 IF ( ln_bt_fw ) THEN 1112 ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2) 1113 ELSE 1114 ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2) 1115 ENDIF 1116 ELSE 1117 zrhot = Agrif_rhot() 1118 ! Time indexes bounds for integration 1119 zt0 = REAL(Agrif_NbStepint() , wp) / zrhot 1120 zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot 1121 ! Polynomial interpolation coefficients: 1122 zat = zrhot * ( zt1**2._wp * (-2._wp*zt1 + 3._wp) & 1123 & - zt0**2._wp * (-2._wp*zt0 + 3._wp) ) 1124 ! 1125 vbdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2) 1126 ! 1127 ! update interpolation stage: 1128 vtint_stage(i1:i2,j1:j2) = 1 1299 1129 ENDIF 1300 1130 ! … … 1302 1132 1303 1133 1304 SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before , nb, ndir)1134 SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before ) 1305 1135 !!---------------------------------------------------------------------- 1306 1136 !! *** ROUTINE interpe3t *** … … 1309 1139 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab 1310 1140 LOGICAL , INTENT(in ) :: before 1311 INTEGER , INTENT(in ) :: nb , ndir1312 1141 ! 1313 1142 INTEGER :: ji, jj, jk 1314 LOGICAL :: western_side, eastern_side, northern_side, southern_side1315 1143 !!---------------------------------------------------------------------- 1316 1144 ! … … 1318 1146 ptab(i1:i2,j1:j2,k1:k2) = tmask(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2) 1319 1147 ELSE 1320 western_side = (nb == 1).AND.(ndir == 1)1321 eastern_side = (nb == 1).AND.(ndir == 2)1322 southern_side = (nb == 2).AND.(ndir == 1)1323 northern_side = (nb == 2).AND.(ndir == 2)1324 1148 ! 1325 1149 DO jk = k1, k2 1326 1150 DO jj = j1, j2 1327 1151 DO ji = i1, i2 1328 !1329 1152 IF( ABS( ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk) ) > 1.D-2) THEN 1330 IF (western_side.AND.(ptab(i1+nbghostcells-1,jj,jk)>0._wp)) THEN 1331 WRITE(numout,*) 'ERROR bathymetry merge at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk 1332 WRITE(numout,*) ptab(ji,jj,jk), e3t_0(ji,jj,jk) 1333 kindic_agr = kindic_agr + 1 1334 ELSEIF (eastern_side.AND.(ptab(i2-nbghostcells+1,jj,jk)>0._wp)) THEN 1335 WRITE(numout,*) 'ERROR bathymetry merge at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk 1336 WRITE(numout,*) ptab(ji,jj,jk), e3t_0(ji,jj,jk) 1337 kindic_agr = kindic_agr + 1 1338 ELSEIF (southern_side.AND.(ptab(ji,j1+nbghostcells-1,jk)>0._wp)) THEN 1339 WRITE(numout,*) 'ERROR bathymetry merge at the southern border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk 1340 WRITE(numout,*) ptab(ji,jj,jk), e3t_0(ji,jj,jk) 1341 kindic_agr = kindic_agr + 1 1342 ELSEIF (northern_side.AND.(ptab(ji,j2-nbghostcells+1,jk)>0._wp)) THEN 1343 WRITE(numout,*) 'ERROR bathymetry merge at the northen border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk 1344 WRITE(numout,*) ptab(ji,jj,jk), e3t_0(ji,jj,jk) 1345 kindic_agr = kindic_agr + 1 1346 ENDIF 1153 WRITE(numout,*) ' Agrif error for e3t_0: parent , child, i, j, k ', & 1154 & ptab(ji,jj,jk), tmask(ji,jj,jk) * e3t_0(ji,jj,jk), & 1155 & ji+nimpp-1, jj+njmpp-1, jk 1156 kindic_agr = kindic_agr + 1 1347 1157 ENDIF 1348 1158 END DO … … 1353 1163 ! 1354 1164 END SUBROUTINE interpe3t 1355 1356 1357 SUBROUTINE interpumsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )1358 !!----------------------------------------------------------------------1359 !! *** ROUTINE interpumsk ***1360 !!----------------------------------------------------------------------1361 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k21362 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab1363 LOGICAL , INTENT(in ) :: before1364 INTEGER , INTENT(in ) :: nb , ndir1365 !1366 INTEGER :: ji, jj, jk1367 LOGICAL :: western_side, eastern_side1368 !!----------------------------------------------------------------------1369 !1370 IF( before ) THEN1371 ptab(i1:i2,j1:j2,k1:k2) = umask(i1:i2,j1:j2,k1:k2)1372 ELSE1373 western_side = (nb == 1).AND.(ndir == 1)1374 eastern_side = (nb == 1).AND.(ndir == 2)1375 DO jk = k1, k21376 DO jj = j1, j21377 DO ji = i1, i21378 ! Velocity mask at boundary edge points:1379 IF (ABS(ptab(ji,jj,jk) - umask(ji,jj,jk)) > 1.D-2) THEN1380 IF (western_side) THEN1381 WRITE(numout,*) 'ERROR with umask at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk1382 WRITE(numout,*) ' masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)1383 kindic_agr = kindic_agr + 11384 ELSEIF (eastern_side) THEN1385 WRITE(numout,*) 'ERROR with umask at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk1386 WRITE(numout,*) ' masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)1387 kindic_agr = kindic_agr + 11388 ENDIF1389 ENDIF1390 END DO1391 END DO1392 END DO1393 !1394 ENDIF1395 !1396 END SUBROUTINE interpumsk1397 1398 1399 SUBROUTINE interpvmsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )1400 !!----------------------------------------------------------------------1401 !! *** ROUTINE interpvmsk ***1402 !!----------------------------------------------------------------------1403 INTEGER , INTENT(in ) :: i1,i2,j1,j2,k1,k21404 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab1405 LOGICAL , INTENT(in ) :: before1406 INTEGER , INTENT(in ) :: nb , ndir1407 !1408 INTEGER :: ji, jj, jk1409 LOGICAL :: northern_side, southern_side1410 !!----------------------------------------------------------------------1411 !1412 IF( before ) THEN1413 ptab(i1:i2,j1:j2,k1:k2) = vmask(i1:i2,j1:j2,k1:k2)1414 ELSE1415 southern_side = (nb == 2).AND.(ndir == 1)1416 northern_side = (nb == 2).AND.(ndir == 2)1417 DO jk = k1, k21418 DO jj = j1, j21419 DO ji = i1, i21420 ! Velocity mask at boundary edge points:1421 IF (ABS(ptab(ji,jj,jk) - vmask(ji,jj,jk)) > 1.D-2) THEN1422 IF (southern_side) THEN1423 WRITE(numout,*) 'ERROR with vmask at the southern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk1424 WRITE(numout,*) ' masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)1425 kindic_agr = kindic_agr + 11426 ELSEIF (northern_side) THEN1427 WRITE(numout,*) 'ERROR with vmask at the northern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk1428 WRITE(numout,*) ' masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)1429 kindic_agr = kindic_agr + 11430 ENDIF1431 ENDIF1432 END DO1433 END DO1434 END DO1435 !1436 ENDIF1437 !1438 END SUBROUTINE interpvmsk1439 1165 1440 1166 … … 1446 1172 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 1447 1173 LOGICAL , INTENT(in ) :: before 1448 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 1449 REAL(wp), DIMENSION(1:jpk) :: h_out 1450 INTEGER :: N_in, N_out, ji, jj, jk 1174 ! 1175 INTEGER :: ji, jj, jk 1176 INTEGER :: N_in, N_out 1177 REAL(wp), DIMENSION(k1:k2) :: tabin, z_in 1178 REAL(wp), DIMENSION(1:jpk) :: z_out 1451 1179 !!---------------------------------------------------------------------- 1452 1180 ! … … 1459 1187 END DO 1460 1188 END DO 1461 #ifdef key_vertical 1189 1190 # if defined key_vertical 1191 ! Interpolate thicknesses 1192 ! Warning: these are masked, hence extrapolated prior interpolation. 1462 1193 DO jk=k1,k2 1463 1194 DO jj=j1,j2 1464 1195 DO ji=i1,i2 1465 ptab(ji,jj,jk,2) = wmask(ji,jj,jk) * e3w(ji,jj,jk,Kmm_a)1196 ptab(ji,jj,jk,2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 1466 1197 END DO 1467 1198 END DO 1468 1199 END DO 1469 #endif 1200 1201 ! Extrapolate thicknesses in partial bottom cells: 1202 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 1203 IF (ln_zps) THEN 1204 DO jj=j1,j2 1205 DO ji=i1,i2 1206 jk = mbkt(ji,jj) 1207 ptab(ji,jj,jk,2) = 0._wp 1208 END DO 1209 END DO 1210 END IF 1211 1212 ! Save ssh at last level: 1213 IF (.NOT.ln_linssh) THEN 1214 ptab(i1:i2,j1:j2,k2,2) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 1215 ELSE 1216 ptab(i1:i2,j1:j2,k2,2) = 0._wp 1217 END IF 1218 # endif 1470 1219 ELSE 1471 1220 #ifdef key_vertical 1472 avm_k(i1:i2,j1:j2,1:jpk) = 0.1473 DO jj=j1,j21474 DO ji=i1,i21475 N_in = 01476 DO jk=k1,k2 !k2 = jpk of parent grid1477 IF (ptab(ji,jj,jk,2) == 0) EXIT1478 N_in = N_in + 11479 tabin(jk) = ptab(ji,jj,jk,1)1480 h_in(N_in) = ptab(ji,jj,jk,2)1481 END DO1482 N_out = 01483 DO jk=1,jpk ! jpk of child grid1484 IF (wmask(ji,jj,jk) == 0) EXIT1485 N_out = N_out + 11486 h_out(jk) = e3t(ji,jj,jk,Kmm_a)1221 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 1222 avm_k(i1:i2,j1:j2,k1:k2) = 0._wp 1223 1224 DO jj = j1, j2 1225 DO ji =i1, i2 1226 N_in = mbkt_parent(ji,jj) 1227 IF ( tmask(ji,jj,1) == 0._wp) N_in = 0 1228 z_in(N_in+1) = ht0_parent(ji,jj) + ptab(ji,jj,k2,2) 1229 DO jk = N_in, 1, -1 ! Parent vertical grid 1230 z_in(jk) = z_in(jk+1) - ptab(ji,jj,jk,2) 1231 tabin(jk) = ptab(ji,jj,jk,1) 1232 END DO 1233 N_out = mbkt(ji,jj) 1234 DO jk = 1, N_out ! Child vertical grid 1235 z_out(jk) = gdepw(ji,jj,jk,Kmm_a) 1487 1236 ENDDO 1488 IF (N_in > 0) THEN1489 CALL re constructandremap(tabin(1:N_in),h_in,avm_k(ji,jj,1:N_out),h_out,N_in,N_out)1237 IF (N_in*N_out > 0) THEN 1238 CALL remap_linear(tabin(1:N_in),z_in(1:N_in),avm_k(ji,jj,1:N_out),z_out(1:N_out),N_in,N_out,1) 1490 1239 ENDIF 1491 1240 ENDDO … … 1497 1246 ! 1498 1247 END SUBROUTINE interpavm 1248 1249 # if defined key_vertical 1250 SUBROUTINE interpmbkt( ptab, i1, i2, j1, j2, before ) 1251 !!---------------------------------------------------------------------- 1252 !! *** ROUTINE interpsshn *** 1253 !!---------------------------------------------------------------------- 1254 INTEGER , INTENT(in ) :: i1, i2, j1, j2 1255 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 1256 LOGICAL , INTENT(in ) :: before 1257 ! 1258 !!---------------------------------------------------------------------- 1259 ! 1260 IF( before) THEN 1261 ptab(i1:i2,j1:j2) = REAL(mbkt(i1:i2,j1:j2),wp) 1262 ELSE 1263 mbkt_parent(i1:i2,j1:j2) = NINT(ptab(i1:i2,j1:j2)) 1264 ENDIF 1265 ! 1266 END SUBROUTINE interpmbkt 1267 1268 SUBROUTINE interpht0( ptab, i1, i2, j1, j2, before ) 1269 !!---------------------------------------------------------------------- 1270 !! *** ROUTINE interpsshn *** 1271 !!---------------------------------------------------------------------- 1272 INTEGER , INTENT(in ) :: i1, i2, j1, j2 1273 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 1274 LOGICAL , INTENT(in ) :: before 1275 ! 1276 !!---------------------------------------------------------------------- 1277 ! 1278 IF( before) THEN 1279 ptab(i1:i2,j1:j2) = ht_0(i1:i2,j1:j2) 1280 ELSE 1281 ht0_parent(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) 1282 ENDIF 1283 ! 1284 END SUBROUTINE interpht0 1285 #endif 1499 1286 1500 1287 #else -
NEMO/branches/2019/dev_r11943_MERGE_2019/src/NST/agrif_oce_sponge.F90
r11949 r12229 22 22 USE agrif_oce 23 23 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 24 USE iom 25 USE vremap 24 26 25 27 IMPLICIT NONE … … 58 60 #endif 59 61 ! 62 CALL iom_put( 'agrif_spu', fspu(:,:)) 63 CALL iom_put( 'agrif_spv', fspv(:,:)) 64 ! 60 65 END SUBROUTINE Agrif_Sponge_Tra 61 66 … … 85 90 #endif 86 91 ! 92 CALL iom_put( 'agrif_spt', fspt(:,:)) 93 CALL iom_put( 'agrif_spf', fspf(:,:)) 94 ! 87 95 END SUBROUTINE Agrif_Sponge_dyn 88 96 … … 93 101 !!---------------------------------------------------------------------- 94 102 INTEGER :: ji, jj, ind1, ind2 95 INTEGER :: ispongearea 96 REAL(wp) :: z1_ spongearea103 INTEGER :: ispongearea, jspongearea 104 REAL(wp) :: z1_ispongearea, z1_jspongearea 97 105 REAL(wp), DIMENSION(jpi,jpj) :: ztabramp 98 !!---------------------------------------------------------------------- 99 ! 106 REAL(wp), DIMENSION(jpjmax) :: zmskwest, zmskeast 107 REAL(wp), DIMENSION(jpimax) :: zmsknorth, zmsksouth 108 !!---------------------------------------------------------------------- 109 ! 110 ! Sponge 1d example with: 111 ! iraf = 3 ; nbghost = 3 ; nn_sponge_len = 2 112 ! 113 !coarse : U T U T U T U 114 !| | | | | 115 !fine : t u t u t u t u t u t u t u t u t u t u t 116 !sponge val:0 0 0 1 5/6 4/6 3/6 2/6 1/6 0 0 117 ! | ghost | <-- sponge area -- > | 118 ! | points | | 119 ! |--> dynamical interface 120 100 121 #if defined SPONGE || defined SPONGE_TOP 101 122 IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN 123 ! 124 ! Retrieve masks at open boundaries: 125 126 ! --- West --- ! 127 ztabramp(:,:) = 0._wp 128 ind1 = 1+nbghostcells 129 DO ji = mi0(ind1), mi1(ind1) 130 ztabramp(ji,:) = ssumask(ji,:) 131 END DO 132 ! 133 zmskwest(:) = 0._wp 134 zmskwest(1:jpj) = MAXVAL(ztabramp(:,:), dim=1) 135 136 ! --- East --- ! 137 ztabramp(:,:) = 0._wp 138 ind1 = jpiglo - nbghostcells - 1 139 DO ji = mi0(ind1), mi1(ind1) 140 ztabramp(ji,:) = ssumask(ji,:) 141 END DO 142 ! 143 zmskeast(:) = 0._wp 144 zmskeast(1:jpj) = MAXVAL(ztabramp(:,:), dim=1) 145 146 ! --- South --- ! 147 ztabramp(:,:) = 0._wp 148 ind1 = 1+nbghostcells 149 DO jj = mj0(ind1), mj1(ind1) 150 ztabramp(:,jj) = ssvmask(:,jj) 151 END DO 152 ! 153 zmsksouth(:) = 0._wp 154 zmsksouth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2) 155 156 ! --- North --- ! 157 ztabramp(:,:) = 0._wp 158 ind1 = jpjglo - nbghostcells - 1 159 DO jj = mj0(ind1), mj1(ind1) 160 ztabramp(:,jj) = ssvmask(:,jj) 161 END DO 162 ! 163 zmsknorth(:) = 0._wp 164 zmsknorth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2) 165 ! JC: SPONGE MASKING TO BE SORTED OUT: 166 zmskwest(:) = 1._wp 167 zmskeast(:) = 1._wp 168 zmsknorth(:) = 1._wp 169 zmsksouth(:) = 1._wp 170 #if defined key_mpp_mpi 171 ! CALL mpp_max( 'AGRIF_sponge', zmskwest(:) , jpjmax ) 172 ! CALL mpp_max( 'AGRIF_Sponge', zmskeast(:) , jpjmax ) 173 ! CALL mpp_max( 'AGRIF_Sponge', zmsksouth(:), jpimax ) 174 ! CALL mpp_max( 'AGRIF_Sponge', zmsknorth(:), jpimax ) 175 #endif 176 102 177 ! Define ramp from boundaries towards domain interior at T-points 103 178 ! Store it in ztabramp 104 179 105 ispongearea = 1 + nn_sponge_len * Agrif_irhox() 106 z1_spongearea = 1._wp / REAL( ispongearea ) 180 ispongearea = nn_sponge_len * Agrif_irhox() 181 z1_ispongearea = 1._wp / REAL( ispongearea ) 182 jspongearea = nn_sponge_len * Agrif_irhoy() 183 z1_jspongearea = 1._wp / REAL( jspongearea ) 107 184 108 185 ztabramp(:,:) = 0._wp 109 186 187 ! Trick to remove sponge in 2DV domains: 188 IF ( nbcellsx <= 3 ) ispongearea = -1 189 IF ( nbcellsy <= 3 ) jspongearea = -1 190 110 191 ! --- West --- ! 111 IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN 112 ind1 = 1+nbghostcells 113 ind2 = 1+nbghostcells + ispongearea 192 ind1 = 1+nbghostcells 193 ind2 = 1+nbghostcells + ispongearea 194 DO ji = mi0(ind1), mi1(ind2) 195 DO jj = 1, jpj 196 ztabramp(ji,jj) = REAL( ind2 - mig(ji) ) * z1_ispongearea * zmskwest(jj) 197 END DO 198 END DO 199 200 ! ghost cells: 201 ind1 = 1 202 ind2 = nbghostcells + 1 203 DO ji = mi0(ind1), mi1(ind2) 204 DO jj = 1, jpj 205 ztabramp(ji,jj) = zmskwest(jj) 206 END DO 207 END DO 208 209 ! --- East --- ! 210 ind1 = jpiglo - nbghostcells - ispongearea 211 ind2 = jpiglo - nbghostcells 212 DO ji = mi0(ind1), mi1(ind2) 114 213 DO jj = 1, jpj 115 DO ji = ind1, ind2 116 ztabramp(ji,jj) = REAL( ind2 - ji ) * z1_spongearea * umask(ind1,jj,1) 117 END DO 214 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mig(ji) - ind1 ) * z1_ispongearea) * zmskeast(jj) 118 215 ENDDO 119 END IF120 121 ! --- East --- !122 IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN123 ind1 = nlci - nbghostcells - ispongearea124 ind2 = nlci - nbghostcells216 END DO 217 218 ! ghost cells: 219 ind1 = jpiglo - nbghostcells 220 ind2 = jpiglo 221 DO ji = mi0(ind1), mi1(ind2) 125 222 DO jj = 1, jpj 126 DO ji = ind1, ind2 127 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ji - ind1 ) * z1_spongearea * umask(ind2-1,jj,1) ) 128 ENDDO 223 ztabramp(ji,jj) = zmskeast(jj) 129 224 ENDDO 130 END IF225 END DO 131 226 132 227 ! --- South --- ! 133 IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN 134 ind1 = 1+nbghostcells 135 ind2 = 1+nbghostcells + ispongearea 136 DO jj = ind1, ind2 137 DO ji = 1, jpi 138 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ind2 - jj ) * z1_spongearea * vmask(ji,ind1,1) ) 139 END DO 140 ENDDO 141 ENDIF 228 ind1 = 1+nbghostcells 229 ind2 = 1+nbghostcells + jspongearea 230 DO jj = mj0(ind1), mj1(ind2) 231 DO ji = 1, jpi 232 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ind2 - mjg(jj) ) * z1_jspongearea) * zmsksouth(ji) 233 END DO 234 END DO 235 236 ! ghost cells: 237 ind1 = 1 238 ind2 = nbghostcells + 1 239 DO jj = mj0(ind1), mj1(ind2) 240 DO ji = 1, jpi 241 ztabramp(ji,jj) = zmsksouth(ji) 242 END DO 243 END DO 142 244 143 245 ! --- North --- ! 144 IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN 145 ind1 = nlcj - nbghostcells - ispongearea 146 ind2 = nlcj - nbghostcells 147 DO jj = ind1, ind2 148 DO ji = 1, jpi 149 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( jj - ind1 ) * z1_spongearea * vmask(ji,ind2-1,1) ) 150 END DO 151 ENDDO 152 ENDIF 246 ind1 = jpjglo - nbghostcells - jspongearea 247 ind2 = jpjglo - nbghostcells 248 DO jj = mj0(ind1), mj1(ind2) 249 DO ji = 1, jpi 250 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mjg(jj) - ind1 ) * z1_jspongearea) * zmsknorth(ji) 251 END DO 252 END DO 253 254 ! ghost cells: 255 ind1 = jpjglo - nbghostcells 256 ind2 = jpjglo 257 DO jj = mj0(ind1), mj1(ind2) 258 DO ji = 1, jpi 259 ztabramp(ji,jj) = zmsknorth(ji) 260 END DO 261 END DO 153 262 154 263 ENDIF … … 156 265 ! Tracers 157 266 IF( .NOT. spongedoneT ) THEN 158 fs aht_spu(:,:) = 0._wp159 fs aht_spv(:,:) = 0._wp267 fspu(:,:) = 0._wp 268 fspv(:,:) = 0._wp 160 269 DO jj = 2, jpjm1 161 270 DO ji = 2, jpim1 ! vector opt. 162 fs aht_spu(ji,jj) = 0.5_wp * visc_tra * ( ztabramp(ji,jj) + ztabramp(ji+1,jj ))163 fs aht_spv(ji,jj) = 0.5_wp * visc_tra * ( ztabramp(ji,jj) + ztabramp(ji ,jj+1))164 END DO 165 END DO 166 CALL lbc_lnk( 'agrif_ oce_sponge', fsaht_spu, 'U', 1. ) ! Lateral boundary conditions167 CALL lbc_lnk( 'agrif_ oce_sponge', fsaht_spv, 'V', 1. )168 271 fspu(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji+1,jj ) ) * ssumask(ji,jj) 272 fspv(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji ,jj+1) ) * ssvmask(ji,jj) 273 END DO 274 END DO 275 CALL lbc_lnk( 'agrif_Sponge', fspu, 'U', 1. ) ! Lateral boundary conditions 276 CALL lbc_lnk( 'agrif_Sponge', fspv, 'V', 1. ) 277 169 278 spongedoneT = .TRUE. 170 279 ENDIF … … 172 281 ! Dynamics 173 282 IF( .NOT. spongedoneU ) THEN 174 fs ahm_spt(:,:) = 0._wp175 fs ahm_spf(:,:) = 0._wp283 fspt(:,:) = 0._wp 284 fspf(:,:) = 0._wp 176 285 DO jj = 2, jpjm1 177 286 DO ji = 2, jpim1 ! vector opt. 178 fsahm_spt(ji,jj) = visc_dyn * ztabramp(ji,jj) 179 fsahm_spf(ji,jj) = 0.25_wp * visc_dyn * ( ztabramp(ji ,jj ) + ztabramp(ji ,jj+1) & 180 & +ztabramp(ji+1,jj+1) + ztabramp(ji+1,jj ) ) 181 END DO 182 END DO 183 CALL lbc_lnk( 'agrif_oce_sponge', fsahm_spt, 'T', 1. ) ! Lateral boundary conditions 184 CALL lbc_lnk( 'agrif_oce_sponge', fsahm_spf, 'F', 1. ) 287 fspt(ji,jj) = ztabramp(ji,jj) * ssmask(ji,jj) 288 fspf(ji,jj) = 0.25_wp * ( ztabramp(ji ,jj ) + ztabramp(ji ,jj+1) & 289 & +ztabramp(ji+1,jj+1) + ztabramp(ji+1,jj ) ) & 290 & * ssvmask(ji,jj) * ssvmask(ji,jj+1) 291 END DO 292 END DO 293 CALL lbc_lnk( 'agrif_Sponge', fspt, 'T', 1. ) ! Lateral boundary conditions 294 CALL lbc_lnk( 'agrif_Sponge', fspf, 'F', 1. ) 185 295 186 296 spongedoneU = .TRUE. 187 297 ENDIF 188 ! 298 299 #if defined key_vertical 300 ! Remove vertical interpolation where not needed: 301 DO jj = 2, jpjm1 302 DO ji = 2, jpim1 303 IF ((fspu(ji-1,jj)==0._wp).AND.(fspu(ji,jj)==0._wp).AND. & 304 & (fspv(ji,jj-1)==0._wp).AND.(fspv(ji,jj)==0._wp)) mbkt_parent(ji,jj) = 0 305 ! 306 IF ((fspt(ji+1,jj)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. & 307 & (fspf(ji,jj-1)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbku_parent(ji,jj) = 0 308 ! 309 IF ((fspt(ji,jj+1)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. & 310 & (fspf(ji-1,jj)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbkv_parent(ji,jj) = 0 311 ! 312 IF ( ssmask(ji,jj) == 0._wp) mbkt_parent(ji,jj) = 0 313 IF (ssumask(ji,jj) == 0._wp) mbku_parent(ji,jj) = 0 314 IF (ssvmask(ji,jj) == 0._wp) mbkv_parent(ji,jj) = 0 315 END DO 316 END DO 317 ! 318 ztabramp(:,:) = REAL( mbkt_parent(:,:), wp ) ; CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'T', 1. ) 319 mbkt_parent(:,:) = NINT( ztabramp(:,:) ) 320 ztabramp(:,:) = REAL( mbku_parent(:,:), wp ) ; CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'U', 1. ) 321 mbku_parent(:,:) = NINT( ztabramp(:,:) ) 322 ztabramp(:,:) = REAL( mbkv_parent(:,:), wp ) ; CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'V', 1. ) 323 mbkv_parent(:,:) = NINT( ztabramp(:,:) ) 189 324 #endif 190 325 ! 326 #endif 327 ! 191 328 END SUBROUTINE Agrif_Sponge 192 329 193 SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )330 SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 194 331 !!---------------------------------------------------------------------- 195 332 !! *** ROUTINE interptsn_sponge *** … … 201 338 INTEGER :: ji, jj, jk, jn ! dummy loop indices 202 339 INTEGER :: iku, ikv 203 REAL(wp) :: ztsa, zabe1, zabe2, zbtr 340 REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot, ztrelax 204 341 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv 205 342 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff … … 210 347 REAL(wp), DIMENSION(1:jpk) :: h_out 211 348 INTEGER :: N_in, N_out 212 REAL(wp) :: h_diff213 349 !!---------------------------------------------------------------------- 214 350 ! … … 225 361 226 362 # if defined key_vertical 227 DO jk=k1,k2 228 DO jj=j1,j2 229 DO ji=i1,i2 230 tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 231 END DO 232 END DO 233 END DO 363 ! Interpolate thicknesses 364 ! Warning: these are masked, hence extrapolated prior interpolation. 365 DO jk=k1,k2 366 DO jj=j1,j2 367 DO ji=i1,i2 368 tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kbb_a) 369 END DO 370 END DO 371 END DO 372 373 ! Extrapolate thicknesses in partial bottom cells: 374 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 375 IF (ln_zps) THEN 376 DO jj=j1,j2 377 DO ji=i1,i2 378 jk = mbkt(ji,jj) 379 tabres(ji,jj,jk,jpts+1) = 0._wp 380 END DO 381 END DO 382 END IF 383 384 ! Save ssh at last level: 385 IF (.NOT.ln_linssh) THEN 386 tabres(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1) 387 ELSE 388 tabres(i1:i2,j1:j2,k2,jpts+1) = 0._wp 389 END IF 234 390 # endif 235 391 … … 237 393 ! 238 394 # if defined key_vertical 239 tabres_child(:,:,:,:) = 0. 395 396 IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp 397 240 398 DO jj=j1,j2 241 399 DO ji=i1,i2 242 N_in = 0 243 DO jk=k1,k2 !k2 = jpk of parent grid 244 IF (tabres(ji,jj,jk,n2) == 0) EXIT 245 N_in = N_in + 1 400 tabres_child(ji,jj,:,:) = 0._wp 401 N_in = mbkt_parent(ji,jj) 402 zhtot = 0._wp 403 DO jk=1,N_in !k2 = jpk of parent grid 404 IF (jk==N_in) THEN 405 h_in(jk) = ht0_parent(ji,jj) + tabres(ji,jj,k2,n2) - zhtot 406 ELSE 407 h_in(jk) = tabres(ji,jj,jk,n2) 408 ENDIF 409 zhtot = zhtot + h_in(jk) 246 410 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1) 247 h_in(N_in) = tabres(ji,jj,jk,n2)248 411 END DO 249 412 N_out = 0 … … 251 414 IF (tmask(ji,jj,jk) == 0) EXIT 252 415 N_out = N_out + 1 253 h_out(jk) = e3t(ji,jj,jk,K mm_a) !Child grid scale factors. Could multiply by e1e2t here instead of division above416 h_out(jk) = e3t(ji,jj,jk,Kbb_a) !Child grid scale factors. Could multiply by e1e2t here instead of division above 254 417 ENDDO 255 IF (N_in > 0) THEN 256 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 257 tabres(ji,jj,k2,:) = tabres(ji,jj,k2-1,:) !what is this line for????? 258 DO jn=1,jpts 259 call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 260 ENDDO 418 419 ! Account for small differences in free-surface 420 IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 421 h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 422 ELSE 423 h_in(1) = h_in(1) - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 424 ENDIF 425 IF (N_in*N_out > 0) THEN 426 CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts) 261 427 ENDIF 262 428 ENDDO … … 268 434 DO jk=1,jpkm1 269 435 # if defined key_vertical 270 tsbdiff(ji,jj,jk,1:jpts) = ts(ji,jj,jk,1:jpts,Kbb_a) - tabres_child(ji,jj,jk,1:jpts)436 tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres_child(ji,jj,jk,1:jpts)) * tmask(ji,jj,jk) 271 437 # else 272 tsbdiff(ji,jj,jk,1:jpts) = ts(ji,jj,jk,1:jpts,Kbb_a) - tabres(ji,jj,jk,1:jpts)438 tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk) 273 439 # endif 274 440 ENDDO 275 441 ENDDO 276 442 ENDDO 443 444 !* set relaxation time scale 445 IF( neuler == 0 .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_tra / ( rdt ) 446 ELSE ; ztrelax = rn_trelax_tra / (2._wp * rdt ) 447 ENDIF 277 448 278 449 DO jn = 1, jpts … … 281 452 DO jj = j1,j2 282 453 DO ji = i1,i2-1 283 zabe1 = fsaht_spu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm_a)454 zabe1 = rn_sponge_tra * fspu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 284 455 ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj ,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 285 456 END DO … … 288 459 DO ji = i1,i2 289 460 DO jj = j1,j2-1 290 zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a)461 zabe2 = rn_sponge_tra * fspv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a) 291 462 ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 292 463 END DO … … 312 483 zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a) 313 484 ! horizontal diffusive trends 314 ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) 485 ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) & 486 & - ztrelax * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn) 315 487 ! add it to the general tracer trends 316 488 ts(ji,jj,jk,jn,Krhs_a) = ts(ji,jj,jk,jn,Krhs_a) + ztsa … … 339 511 340 512 ! sponge parameters 341 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff513 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot, ztrelax 342 514 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff 343 515 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff … … 346 518 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 347 519 REAL(wp), DIMENSION(1:jpk) :: h_out 348 INTEGER ::N_in, N_out520 INTEGER ::N_in, N_out 349 521 !!--------------------------------------------- 350 522 ! 351 523 IF( before ) THEN 352 DO jk= 1,jpkm1524 DO jk=k1,k2 353 525 DO jj=j1,j2 354 526 DO ji=i1,i2 355 527 tabres(ji,jj,jk,m1) = uu(ji,jj,jk,Kbb_a) 356 528 # if defined key_vertical 357 tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,K mm_a)*umask(ji,jj,jk)529 tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kbb_a)*umask(ji,jj,jk) 358 530 # endif 359 531 END DO … … 361 533 END DO 362 534 535 # if defined key_vertical 536 ! Extrapolate thicknesses in partial bottom cells: 537 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 538 IF (ln_zps) THEN 539 DO jj=j1,j2 540 DO ji=i1,i2 541 jk = mbku(ji,jj) 542 tabres(ji,jj,jk,m2) = 0._wp 543 END DO 544 END DO 545 END IF 546 ! Save ssh at last level: 547 tabres(i1:i2,j1:j2,k2,m2) = 0._wp 548 IF (.NOT.ln_linssh) THEN 549 ! This vertical sum below should be replaced by the sea-level at U-points (optimization): 550 DO jk=1,jpk 551 tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3u(i1:i2,j1:j2,jk,Kbb_a) * umask(i1:i2,j1:j2,jk) 552 END DO 553 tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hu_0(i1:i2,j1:j2) 554 END IF 555 # endif 556 363 557 ELSE 364 558 365 559 # if defined key_vertical 366 tabres_child(:,:,:) = 0._wp 560 IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp 561 367 562 DO jj=j1,j2 368 563 DO ji=i1,i2 369 N_in = 0 370 DO jk=k1,k2 371 IF (tabres(ji,jj,jk,m2) == 0) EXIT 372 N_in = N_in + 1 564 tabres_child(ji,jj,:) = 0._wp 565 N_in = mbku_parent(ji,jj) 566 zhtot = 0._wp 567 DO jk=1,N_in 568 IF (jk==N_in) THEN 569 h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 570 ELSE 571 h_in(jk) = tabres(ji,jj,jk,m2) 572 ENDIF 573 zhtot = zhtot + h_in(jk) 373 574 tabin(jk) = tabres(ji,jj,jk,m1) 374 h_in(N_in) = tabres(ji,jj,jk,m2) 375 ENDDO 376 ! 377 IF (N_in == 0) THEN 378 tabres_child(ji,jj,:) = 0. 379 CYCLE 380 ENDIF 381 382 N_out = 0 383 DO jk=1,jpk 384 if (umask(ji,jj,jk) == 0) EXIT 385 N_out = N_out + 1 386 h_out(N_out) = e3u(ji,jj,jk,Kmm_a) 387 ENDDO 388 389 IF (N_out == 0) THEN 390 tabres_child(ji,jj,:) = 0. 391 CYCLE 392 ENDIF 393 394 IF (N_in * N_out > 0) THEN 395 h_diff = SUM( h_out(1:N_out) ) - SUM( h_in(1:N_in) ) 396 if (h_diff < -1.e4) then 397 print *,'CHECK YOUR BATHY ...', h_diff, SUM( h_out(1:N_out) ), SUM( h_in(1:N_in) ) 398 endif 399 ENDIF 400 call reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 401 575 ENDDO 576 ! 577 N_out = 0 578 DO jk=1,jpk 579 IF (umask(ji,jj,jk) == 0) EXIT 580 N_out = N_out + 1 581 h_out(N_out) = e3u(ji,jj,jk,Kbb_a) 582 ENDDO 583 584 ! Account for small differences in free-surface 585 IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 586 h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 587 ELSE 588 h_in(1) = h_in(1) - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 589 ENDIF 590 591 IF (N_in * N_out > 0) THEN 592 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 593 ENDIF 402 594 ENDDO 403 595 ENDDO 404 596 405 ubdiff(i1:i2,j1:j2,:) = ( uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,: ))*umask(i1:i2,j1:j2,:)597 ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:) 406 598 #else 407 ubdiff(i1:i2,j1:j2,:) = ( uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)599 ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:) 408 600 #endif 601 !* set relaxation time scale 602 IF( neuler == 0 .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_dyn / ( rdt ) 603 ELSE ; ztrelax = rn_trelax_dyn / (2._wp * rdt ) 604 ENDIF 409 605 ! 410 606 DO jk = 1, jpkm1 ! Horizontal slab … … 416 612 DO jj = j1,j2 417 613 DO ji = i1+1,i2 ! vector opt. 418 zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,K mm_a) * fsahm_spt(ji,jj)419 hdivdiff(ji,jj,jk) = ( e2u(ji ,jj)*e3u(ji ,jj,jk,K mm_a) * ubdiff(ji ,jj,jk) &420 & -e2u(ji-1,jj)*e3u(ji-1,jj,jk,K mm_a) * ubdiff(ji-1,jj,jk) ) * zbtr614 zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb_a) * rn_sponge_dyn * fspt(ji,jj) 615 hdivdiff(ji,jj,jk) = ( e2u(ji ,jj)*e3u(ji ,jj,jk,Kbb_a) * ubdiff(ji ,jj,jk) & 616 & -e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb_a) * ubdiff(ji-1,jj,jk) ) * zbtr 421 617 END DO 422 618 END DO … … 424 620 DO jj = j1,j2-1 425 621 DO ji = i1,i2 ! vector opt. 426 zbtr = r1_e1e2f(ji,jj) * e3f(ji,jj,jk) * fsahm_spf(ji,jj)622 zbtr = r1_e1e2f(ji,jj) * e3f(ji,jj,jk) * rn_sponge_dyn * fspf(ji,jj) 427 623 rotdiff(ji,jj,jk) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1,jk) & 428 624 & +e1u(ji,jj ) * ubdiff(ji,jj ,jk) ) * fmask(ji,jj,jk) * zbtr … … 440 636 ! horizontal diffusive trends 441 637 zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) ) & 442 + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj) 638 & + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj) & 639 & - ztrelax * fspu(ji,jj) * ubdiff(ji,jj,jk) 443 640 444 641 ! add it to the general momentum trends 445 uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) + zua 446 642 uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) + zua 447 643 END DO 448 644 ENDIF … … 492 688 ! 493 689 INTEGER :: ji, jj, jk, imax 494 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff690 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot, ztrelax 495 691 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff 496 692 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff … … 503 699 504 700 IF( before ) THEN 505 DO jk= 1,jpkm1701 DO jk=k1,k2 506 702 DO jj=j1,j2 507 703 DO ji=i1,i2 508 704 tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb_a) 509 705 # if defined key_vertical 510 tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v(ji,jj,jk,K mm_a)706 tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v(ji,jj,jk,Kbb_a) 511 707 # endif 512 708 END DO 513 709 END DO 514 710 END DO 711 712 # if defined key_vertical 713 ! Extrapolate thicknesses in partial bottom cells: 714 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 715 IF (ln_zps) THEN 716 DO jj=j1,j2 717 DO ji=i1,i2 718 jk = mbkv(ji,jj) 719 tabres(ji,jj,jk,m2) = 0._wp 720 END DO 721 END DO 722 END IF 723 ! Save ssh at last level: 724 tabres(i1:i2,j1:j2,k2,m2) = 0._wp 725 IF (.NOT.ln_linssh) THEN 726 ! This vertical sum below should be replaced by the sea-level at V-points (optimization): 727 DO jk=1,jpk 728 tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3v(i1:i2,j1:j2,jk,Kbb_a) * vmask(i1:i2,j1:j2,jk) 729 END DO 730 tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hv_0(i1:i2,j1:j2) 731 END IF 732 # endif 733 515 734 ELSE 516 735 517 736 # if defined key_vertical 518 tabres_child(:,:,:) = 0._wp737 IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp 519 738 DO jj=j1,j2 520 739 DO ji=i1,i2 521 N_in = 0 522 DO jk=k1,k2 523 IF (tabres(ji,jj,jk,m2) == 0) EXIT 524 N_in = N_in + 1 740 tabres_child(ji,jj,:) = 0._wp 741 N_in = mbkv_parent(ji,jj) 742 zhtot = 0._wp 743 DO jk=1,N_in 744 IF (jk==N_in) THEN 745 h_in(jk) = hv0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 746 ELSE 747 h_in(jk) = tabres(ji,jj,jk,m2) 748 ENDIF 749 zhtot = zhtot + h_in(jk) 525 750 tabin(jk) = tabres(ji,jj,jk,m1) 526 h_in(N_in) = tabres(ji,jj,jk,m2) 527 ENDDO 751 ENDDO 752 ! 753 N_out = 0 754 DO jk=1,jpk 755 IF (vmask(ji,jj,jk) == 0) EXIT 756 N_out = N_out + 1 757 h_out(N_out) = e3v(ji,jj,jk,Kbb_a) 758 ENDDO 759 760 ! Account for small differences in free-surface 761 IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 762 h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 763 ELSE 764 h_in(1) = h_in(1) - ( sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 765 ENDIF 528 766 529 IF (N_in == 0) THEN 530 tabres_child(ji,jj,:) = 0. 531 CYCLE 532 ENDIF 533 534 N_out = 0 535 DO jk=1,jpk 536 if (vmask(ji,jj,jk) == 0) EXIT 537 N_out = N_out + 1 538 h_out(N_out) = e3v(ji,jj,jk,Kmm_a) 539 ENDDO 540 541 IF (N_in * N_out > 0) THEN 542 h_diff = SUM( h_out(1:N_out) ) - SUM( h_in(1:N_in) ) 543 if (h_diff < -1.e4) then 544 print *,'CHECK YOUR BATHY ...', h_diff, SUM( h_out(1:N_out) ), SUM( h_in(1:N_in) ) 545 endif 546 ENDIF 547 call reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 767 IF (N_in * N_out > 0) THEN 768 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 769 ENDIF 548 770 ENDDO 549 771 ENDDO 550 772 551 vbdiff(i1:i2,j1:j2,:) = ( vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,: ))*vmask(i1:i2,j1:j2,:)773 vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:) 552 774 # else 553 vbdiff(i1:i2,j1:j2,:) = ( vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)775 vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:) 554 776 # endif 777 !* set relaxation time scale 778 IF( neuler == 0 .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_dyn / ( rdt ) 779 ELSE ; ztrelax = rn_trelax_dyn / (2._wp * rdt ) 780 ENDIF 555 781 ! 556 782 DO jk = 1, jpkm1 ! Horizontal slab … … 562 788 DO jj = j1+1,j2 563 789 DO ji = i1,i2 ! vector opt. 564 zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,K mm_a) * fsahm_spt(ji,jj)565 hdivdiff(ji,jj,jk) = ( e1v(ji,jj ) * e3v(ji,jj ,jk,K mm_a) * vbdiff(ji,jj ,jk) &566 & -e1v(ji,jj-1) * e3v(ji,jj-1,jk,K mm_a) * vbdiff(ji,jj-1,jk) ) * zbtr790 zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb_a) * rn_sponge_dyn * fspt(ji,jj) 791 hdivdiff(ji,jj,jk) = ( e1v(ji,jj ) * e3v(ji,jj ,jk,Kbb_a) * vbdiff(ji,jj ,jk) & 792 & -e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kbb_a) * vbdiff(ji,jj-1,jk) ) * zbtr 567 793 END DO 568 794 END DO 569 795 DO jj = j1,j2 570 796 DO ji = i1,i2-1 ! vector opt. 571 zbtr = r1_e1e2f(ji,jj) * e3f(ji,jj,jk) * fsahm_spf(ji,jj)572 rotdiff(ji,jj,jk) = ( 573 & - 797 zbtr = r1_e1e2f(ji,jj) * e3f(ji,jj,jk) * rn_sponge_dyn * fspf(ji,jj) 798 rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) & 799 & -e2v(ji ,jj) * vbdiff(ji ,jj,jk) ) * fmask(ji,jj,jk) * zbtr 574 800 END DO 575 801 END DO … … 586 812 IF( .NOT. tabspongedone_u(ji,jj) ) THEN 587 813 DO jk = 1, jpkm1 588 uu(ji,jj,jk,Krhs_a) = uu(ji ,jj,jk,Krhs_a) & 589 & - ( rotdiff(ji ,jj,jk) - rotdiff(ji,jj-1,jk) ) & 590 & / ( e2u(ji ,jj) * e3u(ji,jj ,jk,Kmm_a) ) & 591 & + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj ,jk) ) * r1_e1u(ji,jj) 814 uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) & 815 & - ( rotdiff (ji ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) ) & 816 & + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj ,jk)) * r1_e1u(ji,jj) 592 817 END DO 593 818 ENDIF … … 601 826 IF( .NOT. tabspongedone_v(ji,jj) ) THEN 602 827 DO jk = 1, jpkm1 603 vv(ji,jj,jk,Krhs_a) = vv(ji,jj ,jk,Krhs_a)&604 & + ( rotdiff(ji,jj ,jk) - rotdiff(ji-1,jj,jk) )&605 & / ( e1v(ji,jj ) * e3v(ji ,jj,jk,Kmm_a) )&606 & + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji ,jj,jk) ) * r1_e2v(ji,jj)828 vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) & 829 & + ( rotdiff (ji,jj ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) ) & 830 & + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji ,jj,jk) ) * r1_e2v(ji,jj) & 831 & - ztrelax * fspv(ji,jj) * vbdiff(ji,jj,jk) 607 832 END DO 608 833 ENDIF -
NEMO/branches/2019/dev_r11943_MERGE_2019/src/NST/agrif_oce_update.F90
r11949 r12229 1 # define TWO_WAY /* TWO WAY NESTING*/2 #undef DECAL_FEEDBACK /* SEPARATION of INTERFACES*/3 #undef VOL_REFLUX /* VOLUME REFLUXING*/1 #undef DECAL_FEEDBACK /* SEPARATION of INTERFACES */ 2 #undef DECAL_FEEDBACK_2D /* SEPARATION of INTERFACES (Barotropic mode) */ 3 #undef VOL_REFLUX /* VOLUME REFLUXING*/ 4 4 5 5 MODULE agrif_oce_update … … 25 25 USE lib_mpp ! MPP library 26 26 USE domvvl ! Need interpolation routines 27 USE vremap ! Vertical remapping 27 28 28 29 IMPLICIT NONE … … 46 47 IF (Agrif_Root()) RETURN 47 48 ! 48 #if defined TWO_WAY49 49 IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers from grid Number',Agrif_Fixed() 50 50 51 #if defined key_vertical 52 ! Effect of this has to be carrefully checked 53 ! depending on what the nesting tools ensure for 54 ! volume conservation: 55 Agrif_UseSpecialValueInUpdate = .FALSE. 56 #else 51 57 Agrif_UseSpecialValueInUpdate = .TRUE. 58 #endif 52 59 Agrif_SpecialValueFineGrid = 0._wp 53 60 ! … … 64 71 Agrif_UseSpecialValueInUpdate = .FALSE. 65 72 ! 66 #endif67 73 ! 68 74 END SUBROUTINE Agrif_Update_Tra … … 75 81 IF (Agrif_Root()) RETURN 76 82 ! 77 #if defined TWO_WAY78 83 IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update momentum from grid Number',Agrif_Fixed() 79 84 … … 95 100 # endif 96 101 97 # if ! defined DECAL_FEEDBACK 102 # if ! defined DECAL_FEEDBACK_2D 98 103 CALL Agrif_Update_Variable(e1u_id,procname = updateU2d) 99 104 CALL Agrif_Update_Variable(e2v_id,procname = updateV2d) … … 103 108 # endif 104 109 ! 105 # if ! defined DECAL_FEEDBACK 110 # if ! defined DECAL_FEEDBACK_2D 106 111 ! Account for updated thicknesses at boundary edges 107 112 IF (.NOT.ln_linssh) THEN … … 113 118 IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN 114 119 ! Update time integrated transports 115 # if ! defined DECAL_FEEDBACK 120 # if ! defined DECAL_FEEDBACK_2D 116 121 CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b) 117 122 CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b) … … 121 126 # endif 122 127 END IF 123 #endif124 128 ! 125 129 END SUBROUTINE Agrif_Update_Dyn … … 131 135 ! 132 136 IF (Agrif_Root()) RETURN 133 !134 #if defined TWO_WAY135 137 ! 136 138 Agrif_UseSpecialValueInUpdate = .TRUE. 137 139 Agrif_SpecialValueFineGrid = 0. 138 # if ! defined DECAL_FEEDBACK 140 # if ! defined DECAL_FEEDBACK_2D 139 141 CALL Agrif_Update_Variable(sshn_id,procname = updateSSH) 140 142 # else … … 147 149 IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN 148 150 ! Refluxing on ssh: 149 # if defined DECAL_FEEDBACK 151 # if defined DECAL_FEEDBACK_2D 150 152 CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0, 0/),locupdate2=(/1, 1/),procname = reflux_sshu) 151 153 CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1, 1/),locupdate2=(/0, 0/),procname = reflux_sshv) … … 157 159 # endif 158 160 ! 159 #endif160 !161 161 END SUBROUTINE Agrif_Update_ssh 162 162 … … 170 170 IF (Agrif_Root()) RETURN 171 171 ! 172 # if defined TWO_WAY173 174 172 Agrif_UseSpecialValueInUpdate = .TRUE. 175 173 Agrif_SpecialValueFineGrid = 0. … … 180 178 181 179 Agrif_UseSpecialValueInUpdate = .FALSE. 182 183 # endif184 180 185 181 END SUBROUTINE Agrif_Update_Tke … … 192 188 ! 193 189 IF (Agrif_Root()) RETURN 194 !195 #if defined TWO_WAY196 190 ! 197 191 IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update e3 from grid Number',Agrif_Fixed(), 'Step', Agrif_Nb_Step() … … 209 203 CALL dom_vvl_update_UVF 210 204 CALL Agrif_ParentGrid_To_ChildGrid() 211 !212 #endif213 205 ! 214 206 END SUBROUTINE Agrif_Update_vvl … … 300 292 !! 301 293 INTEGER :: ji,jj,jk,jn 302 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child 294 INTEGER :: N_in, N_out 295 REAL(wp) :: ztb, ztnu, ztno 303 296 REAL(wp) :: h_in(k1:k2) 304 297 REAL(wp) :: h_out(1:jpk) 305 INTEGER :: N_in, N_out 306 REAL(wp) :: zrho_xy, h_diff 307 REAL(wp) :: tabin(k1:k2,n1:n2) 298 REAL(wp) :: tabin(k1:k2,1:jpts) 299 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,1:jpts) :: tabres_child 308 300 !!--------------------------------------------- 309 301 ! 310 302 IF (before) THEN 311 AGRIF_SpecialValue = -999._wp 312 zrho_xy = Agrif_rhox() * Agrif_rhoy() 303 !jc_alt 304 ! AGRIF_SpecialValue = -999._wp 313 305 DO jn = n1,n2-1 314 306 DO jk=k1,k2 315 307 DO jj=j1,j2 316 308 DO ji=i1,i2 317 tabres(ji,jj,jk,jn) = ( ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) ) & 318 & * tmask(ji,jj,jk) + (tmask(ji,jj,jk) - 1._wp)*999._wp 309 !jc_alt 310 ! tabres(ji,jj,jk,jn) = (ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) ) & 311 ! & * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1._wp) * 999._wp 312 tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) 319 313 END DO 320 314 END DO … … 324 318 DO jj=j1,j2 325 319 DO ji=i1,i2 326 tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) & 327 & + (tmask(ji,jj,jk) - 1._wp)*999._wp 320 !jc_alt 321 ! tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) & 322 ! & + (tmask(ji,jj,jk) - 1._wp) * 999._wp 323 tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 328 324 END DO 329 325 END DO … … 336 332 N_in = 0 337 333 DO jk=k1,k2 !k2 = jpk of child grid 338 IF (tabres(ji,jj,jk,n2) == 0 ) EXIT 334 ! jc_alt 335 ! IF (tabres(ji,jj,jk,n2) < -900._wp ) EXIT 336 IF (tabres(ji,jj,jk,n2) == 0._wp ) EXIT 339 337 N_in = N_in + 1 340 338 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) … … 343 341 N_out = 0 344 342 DO jk=1,jpk ! jpk of parent grid 345 IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF343 IF (tmask(ji,jj,jk) == 0 ) EXIT ! TODO: Will not work with ISF 346 344 N_out = N_out + 1 347 345 h_out(N_out) = e3t(ji,jj,jk,Kmm_a) 348 346 ENDDO 349 IF (N_in > 0) THEN !Remove this? 350 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 351 IF (h_diff < -1.e-4) THEN 352 print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out)) 353 print *,h_in(1:N_in) 354 print *,h_out(1:N_out) 355 STOP 356 ENDIF 357 DO jn=n1,n2-1 358 CALL reconstructandremap( tabin(1:N_in,jn), h_in(1:N_in), tabres_child(ji,jj,1:N_out,jn), & 359 & h_out(1:N_out) , N_in , N_out ) 360 ENDDO 347 IF (N_in*N_out > 0) THEN !Remove this? 348 CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts) 361 349 ENDIF 362 350 ENDDO … … 365 353 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 366 354 ! Add asselin part 367 DO jn = n1,n2-1 368 DO jk=1,jpk 369 DO jj=j1,j2 370 DO ji=i1,i2 371 IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 372 ts(ji,jj,jk,jn,Kbb_a) = ts(ji,jj,jk,jn,Kbb_a) & 373 & + atfp * ( tabres_child(ji,jj,jk,jn) & 374 & - ts(ji,jj,jk,jn,Kmm_a) & 375 & ) * tmask(ji,jj,jk) 355 DO jn = 1,jpts 356 DO jk = 1, jpkm1 357 DO jj = j1, j2 358 DO ji = i1, i2 359 IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN 360 ztb = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 361 ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a) 362 ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 363 ts(ji,jj,jk,jn,Kbb_a) = ( ztb + atfp * ( ztnu - ztno) ) & 364 & * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 376 365 ENDIF 377 END DO378 END DO379 END DO380 END DO381 ENDIF 382 DO jn = n1,n2-1383 DO jk =1,jpk384 DO jj =j1,j2385 DO ji =i1,i2386 IF( tabres_child(ji,jj,jk,jn) .NE. 0.) THEN387 ts(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk)366 END DO 367 END DO 368 END DO 369 END DO 370 ENDIF 371 DO jn = 1,jpts 372 DO jk = 1, jpkm1 373 DO jj = j1, j2 374 DO ji = i1, i2 375 IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN 376 ts(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn) 388 377 END IF 389 378 END DO … … 391 380 END DO 392 381 END DO 382 ! 383 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 384 ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kbb_a) = ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kmm_a) 385 ENDIF 393 386 ENDIF 394 387 ! … … 416 409 !> jc tmp 417 410 tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk) 418 ! tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a)411 ! tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) 419 412 !< jc tmp 420 413 END DO … … 440 433 ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 441 434 ts(ji,jj,jk,jn,Kbb_a) = ( ztb + atfp * ( ztnu - ztno) ) & 442 & 435 & * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 443 436 ENDIF 444 437 END DO … … 480 473 ! 481 474 INTEGER :: ji, jj, jk 482 REAL(wp):: zrhoy 475 REAL(wp):: zrhoy, zub, zunu, zuno 483 476 ! VERTICAL REFINEMENT BEGIN 484 477 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child … … 493 486 IF( before ) THEN 494 487 zrhoy = Agrif_Rhoy() 495 AGRIF_SpecialValue = -999._wp 488 !jc_alt 489 ! AGRIF_SpecialValue = -999._wp 496 490 DO jk=k1,k2 497 491 DO jj=j1,j2 498 492 DO ji=i1,i2 499 tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm_a) & 500 & + (umask(ji,jj,jk)-1._wp)*999._wp 501 tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) & 502 & + (umask(ji,jj,jk)-1._wp)*999._wp 493 !jc_alt 494 ! tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm_a) & 495 ! & + (umask(ji,jj,jk)-1._wp)*999._wp 496 tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm_a) 497 !jc_alt 498 ! tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) & 499 ! & + (umask(ji,jj,jk)-1._wp)*999._wp 500 tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 503 501 END DO 504 502 END DO … … 513 511 tabin(:) = 0._wp 514 512 DO jk=k1,k2 !k2=jpk of child grid 515 IF( tabres(ji,jj,jk,2) < -900) EXIT 513 !jc_alt 514 ! IF( tabres(ji,jj,jk,2) < -900._wp) EXIT 515 IF( tabres(ji,jj,jk,2) == 0.) EXIT 516 516 N_in = N_in + 1 517 tabin(jk) 517 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 518 518 h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj) 519 519 ENDDO … … 525 525 ENDDO 526 526 IF (N_in * N_out > 0) THEN 527 h_diff = SUM( h_out(1:N_out) ) - SUM( h_in(1:N_in) ) 527 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 528 excess = 0._wp 528 529 IF (h_diff < -1.e-4) THEN 529 530 !Even if bathy at T points match it's possible for the U points to be deeper in the child grid. 530 531 !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 531 excess = 0._wp532 532 DO jk=N_in,1,-1 533 533 thick = MIN(-1*h_diff, h_in(jk)) … … 542 542 ENDDO 543 543 ENDIF 544 CALL reconstructandremap( tabin(1:N_in) , h_in(1:N_in), tabres_child(ji,jj,1:N_out), & 545 & h_out(1:N_out), N_in , N_out ) 546 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/( e2u(ji,jj)*h_out(N_out) ) 544 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 545 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out)) 547 546 ENDIF 548 547 ENDDO 549 548 ENDDO 550 549 ! 551 550 DO jk=1,jpk 552 551 DO jj=j1,j2 553 552 DO ji=i1,i2 554 553 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 555 uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) & 556 & + atfp * ( tabres_child(ji,jj,jk) - uu(ji,jj,jk,Kmm_a) ) * umask(ji,jj,jk) 554 zub = uu(ji,jj,jk,Kbb_a) * e3u(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 555 zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a) 556 zunu = tabres_child(ji,jj,jk) * e3u(ji,jj,jk,Kmm_a) 557 uu(ji,jj,jk,Kbb_a) = ( zub + atfp * ( zunu - zuno) ) & 558 & * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a) 557 559 ENDIF 558 560 ! … … 561 563 END DO 562 564 END DO 565 ! 566 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 567 uu(i1:i2,j1:j2,1:jpkm1,Kbb_a) = uu(i1:i2,j1:j2,1:jpkm1,Kmm_a) 568 ENDIF 569 ! 563 570 ENDIF 564 571 ! … … 582 589 zrhoy = Agrif_Rhoy() 583 590 DO jk = k1, k2 584 tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) & 585 & * uu(i1:i2,j1:j2,jk,Kmm_a) 591 tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) * uu(i1:i2,j1:j2,jk,Kmm_a) 586 592 END DO 587 593 ELSE … … 595 601 zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a) 596 602 zunu = tabres(ji,jj,jk,1) 597 uu(ji,jj,jk,Kbb_a) = ( zub + atfp * ( zunu - zuno ) )&598 &* umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a)603 uu(ji,jj,jk,Kbb_a) = ( zub + atfp * ( zunu - zuno) ) & 604 & * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a) 599 605 ENDIF 600 606 ! … … 669 675 ! 670 676 INTEGER :: ji, jj, jk 671 REAL(wp) :: zrhox 677 REAL(wp) :: zrhox, zvb, zvnu, zvno 672 678 ! VERTICAL REFINEMENT BEGIN 673 679 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child … … 682 688 IF( before ) THEN 683 689 zrhox = Agrif_Rhox() 684 AGRIF_SpecialValue = -999._wp 690 !jc_alt 691 ! AGRIF_SpecialValue = -999._wp 685 692 DO jk=k1,k2 686 693 DO jj=j1,j2 687 694 DO ji=i1,i2 688 tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) & 689 & * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a) & 690 & + (vmask(ji,jj,jk)-1)*999._wp 691 tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) & 692 + (vmask(ji,jj,jk)-1)*999._wp 695 !jc_alt 696 ! tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a) & 697 ! & + (vmask(ji,jj,jk)-1._wp) * 999._wp 698 tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a) 699 !jc_alt 700 ! tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) & 701 ! & + (vmask(ji,jj,jk)-1._wp) * 999._wp 702 tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) 693 703 END DO 694 704 END DO … … 701 711 N_in = 0 702 712 DO jk=k1,k2 703 IF (tabres(ji,jj,jk,2) < -900) EXIT 713 !jc_alt 714 ! IF (tabres(ji,jj,jk,2) < -900._wp) EXIT 715 IF (tabres(ji,jj,jk,2) == 0) EXIT 704 716 N_in = N_in + 1 705 717 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) … … 713 725 ENDDO 714 726 IF (N_in * N_out > 0) THEN 715 h_diff = SUM( h_out(1:N_out) ) - SUM( h_in(1:N_in) ) 727 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 728 excess = 0._wp 716 729 IF (h_diff < -1.e-4) then 717 !Even if bathy at T points match it's possible for the Upoints to be deeper in the child grid.730 !Even if bathy at T points match it's possible for the V points to be deeper in the child grid. 718 731 !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 719 excess = 0._wp720 732 DO jk=N_in,1,-1 721 thick = MIN( -1._wp*h_diff, h_in(jk))733 thick = MIN(-1*h_diff, h_in(jk)) 722 734 excess = excess + tabin(jk)*thick*e2u(ji,jj) 723 tabin(jk) = tabin(jk)*(1. _wp- thick/h_in(jk))735 tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) 724 736 h_diff = h_diff + thick 725 737 IF ( h_diff == 0) THEN … … 730 742 ENDDO 731 743 ENDIF 732 CALL reconstructandremap( tabin(1:N_in) , h_in(1:N_in), tabres_child(ji,jj,1:N_out), & 733 & h_out(1:N_out), N_in , N_out ) 744 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 734 745 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out)) 735 746 ENDIF 736 747 ENDDO 737 748 ENDDO 738 739 DO jk=1,jpk 749 ! 750 DO jk=1,jpkm1 740 751 DO jj=j1,j2 741 752 DO ji=i1,i2 742 ! 743 IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN ! Add asselin part 744 vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) & 745 & + atfp * ( tabres_child(ji,jj,jk) - vv(ji,jj,jk,Kmm_a) ) & 746 & * vmask(ji,jj,jk) 753 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 754 zvb = vv(ji,jj,jk,Kbb_a) * e3v(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 755 zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a) 756 zvnu = tabres_child(ji,jj,jk) * e3v(ji,jj,jk,Kmm_a) 757 vv(ji,jj,jk,Kbb_a) = ( zvb + atfp * ( zvnu - zvno) ) & 758 & * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a) 747 759 ENDIF 748 760 ! … … 751 763 END DO 752 764 END DO 765 ! 766 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 767 vv(i1:i2,j1:j2,1:jpkm1,Kbb_a) = vv(i1:i2,j1:j2,1:jpkm1,Kmm_a) 768 ENDIF 769 ! 753 770 ENDIF 754 771 ! … … 788 805 zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a) 789 806 zvnu = tabres(ji,jj,jk,1) 790 vv(ji,jj,jk,Kbb_a) = ( zvb + atfp * ( zvnu - zvno) ) 791 &* vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a)807 vv(ji,jj,jk,Kbb_a) = ( zvb + atfp * ( zvnu - zvno) ) & 808 & * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a) 792 809 ENDIF 793 810 ! … … 890 907 ! Update barotropic velocities: 891 908 IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 892 IF ( .NOT.( lk_agrif_fstep .AND. (neuler==0) ) ) THEN! Add asselin part909 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 893 910 zcorr = (tabres(ji,jj) - uu_b(ji,jj,Kmm_a) * hu(ji,jj,Krhs_a)) * r1_hu(ji,jj,Kbb_a) 894 911 uu_b(ji,jj,Kbb_a) = uu_b(ji,jj,Kbb_a) + atfp * zcorr * umask(ji,jj,1) … … 955 972 ! 956 973 ! Update barotropic velocities: 957 IF ( .NOT.ln_dynspg_ts .OR. ( ln_dynspg_ts .AND. ( .NOT.ln_bt_fw )) ) THEN958 IF ( .NOT. ( lk_agrif_fstep .AND. ( neuler==0 ) ) ) THEN! Add asselin part974 IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 975 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 959 976 zcorr = (tabres(ji,jj) - vv_b(ji,jj,Kmm_a) * hv(ji,jj,Krhs_a)) * r1_hv(ji,jj,Kbb_a) 960 vv_b(ji,jj,Kbb_a) = vv_b(ji,jj,Kbb_a) + atfp * zcorr *vmask(ji,jj,1)977 vv_b(ji,jj,Kbb_a) = vv_b(ji,jj,Kbb_a) + atfp * zcorr * vmask(ji,jj,1) 961 978 END IF 962 979 ENDIF … … 1007 1024 DO jj=j1,j2 1008 1025 DO ji=i1,i2 1009 ssh(ji,jj,Kbb_a) = ssh(ji,jj,Kbb_a) & 1010 & + atfp * ( tabres(ji,jj) - ssh(ji,jj,Kmm_a) ) & 1011 & * tmask(ji,jj,1) 1026 ssh(ji,jj,Kbb_a) = ssh(ji,jj,Kbb_a) & 1027 & + atfp * ( tabres(ji,jj) - ssh(ji,jj,Kmm_a) ) * tmask(ji,jj,1) 1012 1028 END DO 1013 1029 END DO … … 1103 1119 zcor = rdt * r1_e1e2t(i1 ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj)) 1104 1120 ssh(i1 ,jj,Kmm_a) = ssh(i1 ,jj,Kmm_a) + zcor 1105 IF ( .NOT. ( lk_agrif_fstep .AND. ( neuler==0 ) ) ) & 1106 & ssh(i1 ,jj,Kbb_a) = ssh(i1 ,jj,Kbb_a) + atfp * zcor 1121 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(i1 ,jj,Kbb_a) = ssh(i1 ,jj,Kbb_a) + atfp * zcor 1107 1122 END DO 1108 1123 ENDIF … … 1111 1126 zcor = - rdt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj)) 1112 1127 ssh(i2+1,jj,Kmm_a) = ssh(i2+1,jj,Kmm_a) + zcor 1113 IF ( .NOT. ( lk_agrif_fstep .AND. ( neuler==0 ) ) ) & 1114 & ssh(i2+1,jj,Kbb_a) = ssh(i2+1,jj,Kbb_a) + atfp * zcor 1128 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(i2+1,jj,Kbb_a) = ssh(i2+1,jj,Kbb_a) + atfp * zcor 1115 1129 END DO 1116 1130 ENDIF … … 1191 1205 IF (southern_side) THEN 1192 1206 DO ji=i1,i2 1193 zcor = rdt * r1_e1e2t(ji,j1 ) * e1v(ji,j1 ) * ( vb2_b(ji,j1)-tabres(ji,j1))1207 zcor = rdt * r1_e1e2t(ji,j1 ) * e1v(ji,j1 ) * (vb2_b(ji,j1)-tabres(ji,j1)) 1194 1208 ssh(ji,j1 ,Kmm_a) = ssh(ji,j1 ,Kmm_a) + zcor 1195 IF ( .NOT. ( lk_agrif_fstep .AND. ( neuler==0 ) ) ) & 1196 & ssh(ji,j1 ,Kbb_a) = ssh(ji,j1,Kbb_a) + atfp * zcor 1209 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(ji,j1 ,Kbb_a) = ssh(ji,j1,Kbb_a) + atfp * zcor 1197 1210 END DO 1198 1211 ENDIF 1199 1212 IF (northern_side) THEN 1200 1213 DO ji=i1,i2 1201 zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2 ) * ( vb2_b(ji,j2)-tabres(ji,j2))1214 zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2 ) * (vb2_b(ji,j2)-tabres(ji,j2)) 1202 1215 ssh(ji,j2+1,Kmm_a) = ssh(ji,j2+1,Kmm_a) + zcor 1203 IF ( .NOT. ( lk_agrif_fstep .AND. ( neuler==0 ) ) ) & 1204 & ssh(ji,j2+1,Kbb_a) = ssh(ji,j2+1,Kbb_a) + atfp * zcor 1216 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(ji,j2+1,Kbb_a) = ssh(ji,j2+1,Kbb_a) + atfp * zcor 1205 1217 END DO 1206 1218 ENDIF … … 1235 1247 END DO 1236 1248 END DO 1237 tabres(:,:,:,1) =tabres(:,:,:,1)*Agrif_Rhox()*Agrif_Rhoy()1238 tabres(:,:,:,2) =tabres(:,:,:,2)*Agrif_Rhox()1239 tabres(:,:,:,3) =tabres(:,:,:,3)*Agrif_Rhoy()1249 tabres(:,:,:,1)=tabres(:,:,:,1)*Agrif_Rhox()*Agrif_Rhoy() 1250 tabres(:,:,:,2)=tabres(:,:,:,2)*Agrif_Rhox() 1251 tabres(:,:,:,3)=tabres(:,:,:,3)*Agrif_Rhoy() 1240 1252 ELSE 1241 1253 DO jk=k1,k2 … … 1246 1258 print *,'VAL2 = ',ji,jj,jk,tabres(ji,jj,jk,2),e1t(ji,jj)*tmask(ji,jj,jk) 1247 1259 print *,'VAL3 = ',ji,jj,jk,tabres(ji,jj,jk,3),e2t(ji,jj)*tmask(ji,jj,jk) 1248 ztemp = SQRT( tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3)))1260 ztemp = sqrt(tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3))) 1249 1261 print *,'CORR = ',ztemp-1. 1250 1262 print *,'NEW VALS = ',tabres(ji,jj,jk,2)*ztemp,tabres(ji,jj,jk,3)*ztemp, & 1251 1263 tabres(ji,jj,jk,2)*ztemp*tabres(ji,jj,jk,3)*ztemp 1252 e1t(ji,jj) = tabres(ji,jj,jk,2) *ztemp1253 e2t(ji,jj) = tabres(ji,jj,jk,3) *ztemp1264 e1t(ji,jj) = tabres(ji,jj,jk,2)*ztemp 1265 e2t(ji,jj) = tabres(ji,jj,jk,3)*ztemp 1254 1266 END IF 1255 1267 END DO … … 1331 1343 DO jj=j1,j2 1332 1344 DO ji=i1,i2 1333 ptab(ji,jj,jk) = e3t_0(ji,jj,jk) & 1334 & * ( 1._wp + ssh(ji,jj,Kmm_a) & 1335 & * ssmask(ji,jj) & 1336 & / ( ht_0(ji,jj)-1._wp + ssmask(ji,jj) ) ) 1345 ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + ssh(ji,jj,Kmm_a) & 1346 & *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj))) 1337 1347 END DO 1338 1348 END DO … … 1353 1363 DO jj=j1,j2 1354 1364 DO ji=i1,i2 1355 e3t(ji,jj,jk,Kbb_a) = e3t(ji,jj,jk,Kbb_a)&1356 &+ atfp * ( ptab(ji,jj,jk) - e3t(ji,jj,jk,Kmm_a) )1365 e3t(ji,jj,jk,Kbb_a) = e3t(ji,jj,jk,Kbb_a) & 1366 & + atfp * ( ptab(ji,jj,jk) - e3t(ji,jj,jk,Kmm_a) ) 1357 1367 END DO 1358 1368 END DO … … 1367 1377 DO ji = i1,i2 1368 1378 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 1369 e3w(ji,jj,jk,Kbb_a) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) 1370 & *( e3t(ji,jj,jk-1,Kbb_a) - e3t_0(ji,jj,jk-1) ) &1371 & + 0.5_wp * tmask(ji,jj,jk)&1372 & *( e3t(ji,jj,jk ,Kbb_a) - e3t_0(ji,jj,jk ) )1373 gdepw(ji,jj,jk,Kbb_a) = gdepw(ji,jj,jk-1,Kbb_a) +e3t(ji,jj,jk-1,Kbb_a)1374 gdept(ji,jj,jk,Kbb_a) = zcoef * ( gdepw(ji,jj,jk ,Kbb_a) + 0.5 * e3w(ji,jj,jk ,Kbb_a)) &1375 & + (1._wp - zcoef) * ( gdept(ji,jj,jk-1,Kbb_a) + e3w(ji,jj,jk ,Kbb_a))1379 e3w(ji,jj,jk,Kbb_a) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * & 1380 & ( e3t(ji,jj,jk-1,Kbb_a) - e3t_0(ji,jj,jk-1) ) & 1381 & + 0.5_wp * tmask(ji,jj,jk) * & 1382 & ( e3t(ji,jj,jk ,Kbb_a) - e3t_0(ji,jj,jk ) ) 1383 gdepw(ji,jj,jk,Kbb_a) = gdepw(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk-1,Kbb_a) 1384 gdept(ji,jj,jk,Kbb_a) = zcoef * ( gdepw(ji,jj,jk ,Kbb_a) + 0.5 * e3w(ji,jj,jk,Kbb_a)) & 1385 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb_a) + e3w(ji,jj,jk,Kbb_a)) 1376 1386 END DO 1377 1387 END DO … … 1393 1403 ! 1394 1404 ! Update vertical scale factor at W-points and depths: 1395 e3w(i1:i2,j1:j2,1,Kmm_a) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kmm_a) - e3t_0(i1:i2,j1:j2,1)1405 e3w (i1:i2,j1:j2,1,Kmm_a) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kmm_a) - e3t_0(i1:i2,j1:j2,1) 1396 1406 gdept(i1:i2,j1:j2,1,Kmm_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kmm_a) 1397 1407 gdepw(i1:i2,j1:j2,1,Kmm_a) = 0.0_wp 1398 gde3w(i1:i2,j1:j2,1 ) = gdept(i1:i2,j1:j2,1,Kmm_a) - ( ht(i1:i2,j1:j2) - ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh1408 gde3w(i1:i2,j1:j2,1) = gdept(i1:i2,j1:j2,1,Kmm_a) - (ht(i1:i2,j1:j2)-ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh 1399 1409 ! 1400 1410 DO jk = 2, jpk … … 1402 1412 DO ji = i1,i2 1403 1413 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 1404 e3w(ji,jj,jk,Kmm_a) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) & 1405 & * ( e3t(ji,jj,jk-1,Kmm_a) - e3t_0(ji,jj,jk-1) ) & 1406 & + 0.5_wp * tmask(ji,jj,jk) & 1407 & * ( e3t(ji,jj,jk ,Kmm_a) - e3t_0(ji,jj,jk ) ) 1408 gdepw(ji,jj,jk,Kmm_a) = gdepw(ji,jj,jk-1,Kmm_a) + e3t(ji,jj,jk-1,Kmm_a) 1409 gdept(ji,jj,jk,Kmm_a) = zcoef * ( gdepw(ji,jj,jk ,Kmm_a) + 0.5 * e3w(ji,jj,jk ,Kmm_a) ) & 1410 & + ( 1._wp - zcoef ) * ( gdept(ji,jj,jk-1,Kmm_a) + e3w(ji,jj,jk ,Kmm_a) ) 1411 gde3w(ji,jj,jk ) = gdept(ji,jj,jk ,Kmm_a) - ( ht(ji,jj)-ht_0(ji,jj) ) ! Last term in the rhs is ssh 1414 e3w(ji,jj,jk,Kmm_a) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * ( e3t(ji,jj,jk-1,Kmm_a) - e3t_0(ji,jj,jk-1) ) & 1415 & + 0.5_wp * tmask(ji,jj,jk) * ( e3t(ji,jj,jk ,Kmm_a) - e3t_0(ji,jj,jk ) ) 1416 gdepw(ji,jj,jk,Kmm_a) = gdepw(ji,jj,jk-1,Kmm_a) + e3t(ji,jj,jk-1,Kmm_a) 1417 gdept(ji,jj,jk,Kmm_a) = zcoef * ( gdepw(ji,jj,jk ,Kmm_a) + 0.5 * e3w(ji,jj,jk,Kmm_a)) & 1418 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm_a) + e3w(ji,jj,jk,Kmm_a)) 1419 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm_a) - (ht(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh 1412 1420 END DO 1413 1421 END DO … … 1415 1423 ! 1416 1424 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 1417 e3t(i1:i2,j1:j2,1:jpk,Kbb_a) = e3t(i1:i2,j1:j2,1:jpk,Kmm_a)1418 e3w(i1:i2,j1:j2,1:jpk,Kbb_a) = e3w(i1:i2,j1:j2,1:jpk,Kmm_a)1425 e3t (i1:i2,j1:j2,1:jpk,Kbb_a) = e3t (i1:i2,j1:j2,1:jpk,Kmm_a) 1426 e3w (i1:i2,j1:j2,1:jpk,Kbb_a) = e3w (i1:i2,j1:j2,1:jpk,Kmm_a) 1419 1427 gdepw(i1:i2,j1:j2,1:jpk,Kbb_a) = gdepw(i1:i2,j1:j2,1:jpk,Kmm_a) 1420 1428 gdept(i1:i2,j1:j2,1:jpk,Kbb_a) = gdept(i1:i2,j1:j2,1:jpk,Kmm_a) -
NEMO/branches/2019/dev_r11943_MERGE_2019/src/NST/agrif_top_interp.F90
r11949 r12229 18 18 USE par_trc 19 19 USE trc 20 USE vremap 20 21 ! 21 22 USE lib_mpp ! MPP library … … 48 49 END SUBROUTINE Agrif_trc 49 50 50 SUBROUTINE interptrn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before , nb, ndir)51 SUBROUTINE interptrn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before ) 51 52 !!---------------------------------------------------------------------- 52 53 !! *** ROUTINE interptrn *** … … 55 56 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2 56 57 LOGICAL , INTENT(in ) :: before 57 INTEGER , INTENT(in ) :: nb , ndir58 58 ! 59 INTEGER :: ji, jj, jk, jn, i ref, jref, ibdy, jbdy ! dummy loop indices59 INTEGER :: ji, jj, jk, jn, ibdy, jbdy ! dummy loop indices 60 60 INTEGER :: imin, imax, jmin, jmax, N_in, N_out 61 61 REAL(wp) :: zrho, z1, z2, z3, z4, z5, z6, z7 62 LOGICAL :: western_side, eastern_side,northern_side,southern_side 62 63 63 ! vertical interpolation: 64 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk, n1:n2) :: ptab_child65 REAL(wp), DIMENSION(k1:k2, n1:n2-1) :: tabin64 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,1:jptra) :: ptab_child 65 REAL(wp), DIMENSION(k1:k2,1:jptra) :: tabin 66 66 REAL(wp), DIMENSION(k1:k2) :: h_in 67 67 REAL(wp), DIMENSION(1:jpk) :: h_out 68 REAL(wp) :: h_diff68 !!---------------------------------------------------------------------- 69 69 70 70 IF( before ) THEN … … 90 90 ELSE 91 91 92 western_side = (nb == 1).AND.(ndir == 1) ; eastern_side = (nb == 1).AND.(ndir == 2) 93 southern_side = (nb == 2).AND.(ndir == 1) ; northern_side = (nb == 2).AND.(ndir == 2) 94 95 # if defined key_vertical 92 # if defined key_vertical 96 93 DO jj=j1,j2 97 94 DO ji=i1,i2 98 iref = ji 99 jref = jj 100 if(western_side) iref=MAX(2,ji) 101 if(eastern_side) iref=MIN(nlci-1,ji) 102 if(southern_side) jref=MAX(2,jj) 103 if(northern_side) jref=MIN(nlcj-1,jj) 95 ptab_child(ji,jj,:) = 0._wp 104 96 N_in = 0 105 97 DO jk=k1,k2 !k2 = jpk of parent grid … … 111 103 N_out = 0 112 104 DO jk=1,jpk ! jpk of child grid 113 IF (tmask( iref,jref,jk) == 0) EXIT105 IF (tmask(ji,jj,jk) == 0) EXIT 114 106 N_out = N_out + 1 115 h_out(jk) = e3t( iref,jref,jk,Kmm_a)107 h_out(jk) = e3t(ji,jj,jk,Krhs_a) 116 108 ENDDO 117 109 IF (N_in > 0) THEN 118 DO jn=1,jptra 119 call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 120 ENDDO 110 CALL reconstructandremap(tabin(1:N_in,1:jptra),h_in,ptab_child(ji,jj,1:N_out,1:jptra),h_out,N_in,N_out,jptra) 121 111 ENDIF 122 112 ENDDO … … 129 119 tr(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab_child(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 130 120 END DO 131 132 IF ( .NOT.lk_agrif_clp ) THEN133 !134 imin = i1 ; imax = i2135 jmin = j1 ; jmax = j2136 !137 ! Remove CORNERS138 IF((nbondj == -1).OR.(nbondj == 2)) jmin = 2 + nbghostcells139 IF((nbondj == +1).OR.(nbondj == 2)) jmax = nlcj - nbghostcells - 1140 IF((nbondi == -1).OR.(nbondi == 2)) imin = 2 + nbghostcells141 IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci - nbghostcells - 1142 !143 IF( eastern_side ) THEN144 zrho = Agrif_Rhox()145 z1 = ( zrho - 1._wp ) * 0.5_wp146 z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )147 z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )148 z7 = - ( zrho - 1._wp ) / ( zrho + 3._wp )149 z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7150 !151 ibdy = nlci-nbghostcells152 DO jn = 1, jptra153 tr(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs_a) = z1 * ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) &154 & + z2 * ptab_child(ibdy ,jmin:jmax,1:jpkm1,jn)155 DO jk = 1, jpkm1156 DO jj = jmin,jmax157 IF( umask(ibdy-1,jj,jk) == 0._wp ) THEN158 tr(ibdy,jj,jk,jn,Krhs_a) = tr(ibdy+1,jj,jk,jn,Krhs_a) * tmask(ibdy,jj,jk)159 ELSE160 tr(ibdy,jj,jk,jn,Krhs_a) = ( z4 * tr(ibdy+1,jj,jk,jn,Krhs_a) &161 & + z3 * tr(ibdy-1,jj,jk,jn,Krhs_a) &162 & ) * tmask(ibdy ,jj,jk)163 IF( uu(ibdy-1,jj,jk,Kmm_a) > 0._wp ) THEN164 tr(ibdy,jj,jk,jn,Krhs_a) = ( z6 * tr(ibdy-1,jj,jk,jn,Krhs_a) &165 & + z5 * tr(ibdy+1,jj,jk,jn,Krhs_a) &166 & + z7 * tr(ibdy-2,jj,jk,jn,Krhs_a) &167 & ) * tmask(ibdy ,jj,jk)168 ENDIF169 ENDIF170 END DO171 END DO172 ! Restore ghost points:173 tr(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs_a) = ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) &174 & * tmask(ibdy+1,jmin:jmax,1:jpkm1)175 END DO176 ENDIF177 !178 IF( northern_side ) THEN179 zrho = Agrif_Rhoy()180 z1 = ( zrho - 1._wp ) * 0.5_wp181 z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )182 z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )183 z7 = - ( zrho - 1._wp ) / ( zrho + 3._wp )184 z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7185 !186 jbdy = nlcj-nbghostcells187 DO jn = 1, jptra188 tr(imin:imax,jbdy+1,1:jpkm1,jn,Krhs_a) = z1 * ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) &189 & + z2 * ptab_child(imin:imax,jbdy ,1:jpkm1,jn)190 DO jk = 1, jpkm1191 DO ji = imin,imax192 IF( vmask(ji,jbdy-1,jk) == 0._wp ) THEN193 tr(ji,jbdy,jk,jn,Krhs_a) = tr(ji,jbdy+1,jk,jn,Krhs_a) * tmask(ji,jbdy,jk)194 ELSE195 tr(ji,jbdy,jk,jn,Krhs_a) = ( z4 * tr(ji,jbdy+1,jk,jn,Krhs_a)196 & + z3 * tr(ji,jbdy-1,jk,jn,Krhs_a) ) * tmask(ji,jbdy,jk)197 IF (vv(ji,jbdy-1,jk,Kmm_a) > 0._wp ) THEN198 tr(ji,jbdy,jk,jn,Krhs_a) = ( z6 * tr(ji,jbdy-1,jk,jn,Krhs_a) &199 & + z5 * tr(ji,jbdy+1,jk,jn,Krhs_a) &200 & + z7 * tr(ji,jbdy-2,jk,jn,Krhs_a) ) * tmask(ji,jbdy,jk)201 ENDIF202 ENDIF203 END DO204 END DO205 ! Restore ghost points:206 tr(imin:imax,jbdy+1,1:jpkm1,jn,Krhs_a) = ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) &207 & * tmask(imin:imax,jbdy+1,1:jpkm1)208 END DO209 ENDIF210 !211 IF( western_side ) THEN212 zrho = Agrif_Rhox()213 z1 = ( zrho - 1._wp ) * 0.5_wp214 z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )215 z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )216 z7 = - ( zrho - 1._wp ) / ( zrho + 3._wp )217 z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7218 !219 ibdy = 1+nbghostcells220 DO jn = 1, jptra221 tr(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs_a) = z1 * ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) &222 & + z2 * ptab_child(ibdy ,jmin:jmax,1:jpkm1,jn)223 DO jk = 1, jpkm1224 DO jj = jmin,jmax225 IF( umask(ibdy,jj,jk) == 0._wp ) THEN226 tr(ibdy,jj,jk,jn,Krhs_a) = tr(ibdy-1,jj,jk,jn,Krhs_a) * tmask(ibdy,jj,jk)227 ELSE228 tr(ibdy,jj,jk,jn,Krhs_a) = ( z4 * tr(ibdy-1,jj,jk,jn,Krhs_a) &229 & + z3 * tr(ibdy+1,jj,jk,jn,Krhs_a) ) * tmask(ibdy,jj,jk)230 IF( uu(ibdy,jj,jk,Kmm_a) < 0._wp ) THEN231 tr(ibdy,jj,jk,jn,Krhs_a) = ( z6 * tr(ibdy+1,jj,jk,jn,Krhs_a) &232 & + z5 * tr(ibdy-1,jj,jk,jn,Krhs_a) &233 & + z7 * tr(ibdy+2,jj,jk,jn,Krhs_a) ) * tmask(ibdy,jj,jk)234 ENDIF235 ENDIF236 END DO237 END DO238 ! Restore ghost points:239 tr(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs_a) = ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) &240 & * tmask(ibdy-1,jmin:jmax,1:jpkm1)241 END DO242 ENDIF243 !244 IF( southern_side ) THEN245 zrho = Agrif_Rhoy()246 z1 = ( zrho - 1._wp ) * 0.5_wp247 z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )248 z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )249 z7 = - ( zrho - 1._wp ) / ( zrho + 3._wp )250 z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7251 !252 jbdy=1+nbghostcells253 DO jn = 1, jptra254 tr(imin:imax,jbdy-1,1:jpkm1,jn,Krhs_a) = z1 * ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) &255 & + z2 * ptab_child(imin:imax,jbdy ,1:jpkm1,jn)256 DO jk = 1, jpkm1257 DO ji = imin,imax258 IF( vmask(ji,jbdy,jk) == 0._wp ) THEN259 tr(ji,jbdy,jk,jn,Krhs_a) = tr(ji,jbdy-1,jk,jn,Krhs_a) * tmask(ji,jbdy,jk)260 ELSE261 tr(ji,jbdy,jk,jn,Krhs_a) = ( z4 * tr(ji,jbdy-1,jk,jn,Krhs_a) &262 & + z3 * tr(ji,jbdy+1,jk,jn,Krhs_a) ) * tmask(ji,jbdy,jk)263 IF( vv(ji,jbdy,jk,Kmm_a) < 0._wp ) THEN264 tr(ji,jbdy,jk,jn,Krhs_a) = ( z6 * tr(ji,jbdy+1,jk,jn,Krhs_a) &265 & + z5 * tr(ji,jbdy-1,jk,jn,Krhs_a) &266 & + z7 * tr(ji,jbdy+2,jk,jn,Krhs_a) ) * tmask(ji,jbdy,jk)267 ENDIF268 ENDIF269 END DO270 END DO271 ! Restore ghost points:272 tr(imin:imax,jbdy-1,1:jpkm1,jn,Krhs_a) = ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) &273 & * tmask(imin:imax,jbdy-1,1:jpkm1)274 END DO275 ENDIF276 !277 ENDIF278 121 279 122 ENDIF -
NEMO/branches/2019/dev_r11943_MERGE_2019/src/NST/agrif_top_sponge.F90
r11949 r12229 20 20 USE agrif_oce 21 21 USE agrif_oce_sponge 22 USE vremap 22 23 ! 23 24 USE in_out_manager … … 66 67 ! 67 68 INTEGER :: ji, jj, jk, jn ! dummy loop indices 68 REAL(wp) :: zabe1, zabe2 69 REAL(wp), DIMENSION(i1:i2,j1:j2) :: ztu, ztv70 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2, n1:n2) :: trbdiff69 REAL(wp) :: zabe1, zabe2, ztrelax 70 REAL(wp), DIMENSION(i1:i2,j1:j2) :: ztu, ztv 71 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,1:jptra) :: trbdiff 71 72 ! vertical interpolation: 72 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk, n1:n2) ::tabres_child73 REAL(wp), DIMENSION(k1:k2, n1:n2-1) :: tabin73 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,1:jptra) ::tabres_child 74 REAL(wp), DIMENSION(k1:k2,1:jptra) :: tabin 74 75 REAL(wp), DIMENSION(k1:k2) :: h_in 75 76 REAL(wp), DIMENSION(1:jpk) :: h_out … … 83 84 DO jj=j1,j2 84 85 DO ji=i1,i2 85 tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kbb )86 tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kbb_a) 86 87 END DO 87 88 END DO … … 93 94 DO jj=j1,j2 94 95 DO ji=i1,i2 95 tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,K mm)96 tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kbb_a) 96 97 END DO 97 98 END DO … … 108 109 N_in = N_in + 1 109 110 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1) 110 h_in(N_in) 111 h_in(N_in) = tabres(ji,jj,jk,n2) 111 112 END DO 112 113 N_out = 0 113 114 DO jk=1,jpk ! jpk of child grid 114 115 IF (tmask(ji,jj,jk) == 0) EXIT 115 N_out 116 h_out(jk) = e3t(ji,jj,jk,K mm) !Child grid scale factors. Could multiply by e1e2t here instead of division above116 N_out = N_out + 1 117 h_out(jk) = e3t(ji,jj,jk,Kbb_a) !Child grid scale factors. Could multiply by e1e2t here instead of division above 117 118 ENDDO 118 119 IF (N_in > 0) THEN 119 h_diff = SUM( h_out(1:N_out) ) - SUM( h_in(1:N_in) ) 120 tabres(ji,jj,k2,:) = tabres(ji,jj,k2-1,:) !what is this line for????? 121 DO jn=1,jptra 122 call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 123 ENDDO 120 CALL reconstructandremap(tabin(1:N_in,1:jptra),h_in,tabres_child(ji,jj,1:N_out,1:jptra),h_out,N_in,N_out,jptra) 124 121 ENDIF 125 122 ENDDO … … 131 128 DO jk=1,jpkm1 132 129 # if defined key_vertical 133 trbdiff(ji,jj,jk,1:jptra) = tr(ji,jj,jk,1:jptra,Kbb ) - tabres_child(ji,jj,jk,1:jptra)130 trbdiff(ji,jj,jk,1:jptra) = tr(ji,jj,jk,1:jptra,Kbb_a) - tabres_child(ji,jj,jk,1:jptra) 134 131 # else 135 trbdiff(ji,jj,jk,1:jptra) = tr(ji,jj,jk,1:jptra,Kbb ) -tabres(ji,jj,jk,1:jptra)132 trbdiff(ji,jj,jk,1:jptra) = tr(ji,jj,jk,1:jptra,Kbb_a) - tabres(ji,jj,jk,1:jptra) 136 133 # endif 137 134 ENDDO … … 139 136 ENDDO 140 137 138 !* set relaxation time scale 139 IF( neuler == 0 .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_tra / ( rdt ) 140 ELSE ; ztrelax = rn_trelax_tra / (2._wp * rdt ) 141 ENDIF 142 141 143 DO jn = 1, jptra 142 144 DO jk = 1, jpkm1 143 145 DO jj = j1,j2-1 144 146 DO ji = i1,i2-1 145 zabe1 = fsaht_spu(ji,jj) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk)146 zabe2 = fsaht_spv(ji,jj) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk)147 zabe1 = rn_sponge_tra * fspu(ji,jj) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) 148 zabe2 = rn_sponge_tra * fspv(ji,jj) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) 147 149 ztu(ji,jj) = zabe1 * ( trbdiff(ji+1,jj ,jk,jn) - trbdiff(ji,jj,jk,jn) ) 148 150 ztv(ji,jj) = zabe2 * ( trbdiff(ji ,jj+1,jk,jn) - trbdiff(ji,jj,jk,jn) ) … … 153 155 DO ji = i1+1,i2-1 154 156 IF( .NOT. tabspongedone_trn(ji,jj) ) THEN 155 tr(ji,jj,jk,jn,Krhs) = tr(ji,jj,jk,jn,Krhs) + ( ztu(ji,jj) - ztu(ji-1,jj ) & 156 & + ztv(ji,jj) - ztv(ji ,jj-1) ) & 157 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 157 tr(ji,jj,jk,jn,Krhs_a) = tr(ji,jj,jk,jn,Krhs_a) + ( ztu(ji,jj) - ztu(ji-1,jj ) & 158 & + ztv(ji,jj) - ztv(ji ,jj-1) ) & 159 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a) & 160 & - ztrelax * fspt(ji,jj) * trbdiff(ji,jj,jk,jn) 158 161 ENDIF 159 162 END DO -
NEMO/branches/2019/dev_r11943_MERGE_2019/src/NST/agrif_top_update.F90
r11949 r12229 1 #define TWO_WAY2 1 #undef DECAL_FEEDBACK 3 2 … … 20 19 USE par_trc 21 20 USE trc 21 USE vremap 22 22 23 23 IMPLICIT NONE … … 40 40 IF (Agrif_Root()) RETURN 41 41 ! 42 #if defined TWO_WAY43 42 Agrif_UseSpecialValueInUpdate = .TRUE. 44 43 Agrif_SpecialValueFineGrid = 0._wp … … 53 52 ! 54 53 Agrif_UseSpecialValueInUpdate = .FALSE. 55 !56 #endif57 54 ! 58 55 END SUBROUTINE Agrif_Update_Trc … … 68 65 !! 69 66 INTEGER :: ji,jj,jk,jn 70 REAL(wp) , DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child67 REAL(wp) :: ztb, ztnu, ztno 71 68 REAL(wp) :: h_in(k1:k2) 72 69 REAL(wp) :: h_out(1:jpk) 73 70 INTEGER :: N_in, N_out 74 71 REAL(wp) :: h_diff 75 REAL(wp) :: zrho_xy76 REAL(wp) :: tabin(k1:k2,n1:n2)72 REAL(wp) :: tabin(k1:k2,1:jptra) 73 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,1:jptra) :: tabres_child 77 74 !!--------------------------------------------- 78 75 ! 79 76 IF (before) THEN 80 77 AGRIF_SpecialValue = -999._wp 81 zrho_xy = Agrif_rhox() * Agrif_rhoy()82 78 DO jn = n1,n2-1 83 79 DO jk=k1,k2 84 80 DO jj=j1,j2 85 81 DO ji=i1,i2 86 tabres(ji,jj,jk,jn) = ( tr(ji,jj,jk,jn,Kmm) * e3t(ji,jj,jk,Kmm) ) &87 &* tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp82 tabres(ji,jj,jk,jn) = (tr(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) ) & 83 * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp 88 84 END DO 89 85 END DO … … 93 89 DO jj=j1,j2 94 90 DO ji=i1,i2 95 tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) &96 &+ (tmask(ji,jj,jk)-1)*999._wp91 tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) & 92 + (tmask(ji,jj,jk)-1)*999._wp 97 93 END DO 98 94 END DO … … 114 110 IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF 115 111 N_out = N_out + 1 116 h_out(N_out) = e3t(ji,jj,jk,Kmm ) !Parent grid scale factors. Could multiply by e1e2t here instead of division above112 h_out(N_out) = e3t(ji,jj,jk,Kmm_a) !Parent grid scale factors. Could multiply by e1e2t here instead of division above 117 113 ENDDO 118 114 IF (N_in > 0) THEN !Remove this? … … 124 120 STOP 125 121 ENDIF 126 DO jn=1,jptra 127 CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out) 128 ENDDO 122 CALL reconstructandremap(tabin(1:N_in,1:jptra),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jptra),h_out(1:N_out),N_in,N_out,jptra) 129 123 ENDIF 130 124 ENDDO 131 125 ENDDO 132 126 ! 133 127 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 134 128 ! Add asselin part 135 129 DO jn = 1,jptra 136 DO jk=1,jpk 130 DO jk=1,jpkm1 137 131 DO jj=j1,j2 138 132 DO ji=i1,i2 139 133 IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 140 tr(ji,jj,jk,jn,Kbb) = ts(ji,jj,jk,jn,Kbb) & 141 & + atfp * ( tabres_child(ji,jj,jk,jn) & 142 & - tr(ji,jj,jk,jn,Kmm) & 143 & ) * tmask(ji,jj,jk) 134 ztb = tr(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 135 ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a) 136 ztno = tr(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 137 tr(ji,jj,jk,jn,Kbb_a) = ( ztb + atfp * ( ztnu - ztno) ) & 138 & * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 144 139 ENDIF 145 140 ENDDO … … 149 144 ENDIF 150 145 DO jn = 1,jptra 151 DO jk=1,jpk 146 DO jk=1,jpkm1 152 147 DO jj=j1,j2 153 148 DO ji=i1,i2 154 149 IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 155 tr(ji,jj,jk,jn,Kmm ) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk)150 tr(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn) 156 151 END IF 157 152 END DO … … 159 154 END DO 160 155 END DO 156 ! 157 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 158 tr(i1:i2,j1:j2,1:jpkm1,1:jptra,Kbb_a) = tr(i1:i2,j1:j2,1:jpkm1,1:jptra,Kmm_a) 159 ENDIF 160 ! 161 161 162 ENDIF 162 163 ! … … 184 185 DO ji=i1,i2 185 186 !> jc tmp 186 tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm ) * e3t(ji,jj,jk,Kmm) / e3t_0(ji,jj,jk)187 ! tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm) * e3t(ji,jj,jk,Kmm)187 tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk) 188 ! tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) 188 189 !< jc tmp 189 190 END DO … … 195 196 DO jn = n1,n2 196 197 tabres(i1:i2,j1:j2,k1:k2,jn) = tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) & 197 &* tmask(i1:i2,j1:j2,k1:k2)198 & * tmask(i1:i2,j1:j2,k1:k2) 198 199 ENDDO 199 200 !< jc tmp … … 205 206 DO ji=i1,i2 206 207 IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN 207 ztb = tr(ji,jj,jk,jn,Kbb) * e3t(ji,jj,jk,Kbb) ! fse3t_b prior update should be used208 ztb = tr(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 208 209 ztnu = tabres(ji,jj,jk,jn) 209 ztno = tr(ji,jj,jk,jn,Kmm) * e3t(ji,jj,jk,Krhs)210 tr(ji,jj,jk,jn,Kbb ) = ( ztb + atfp * ( ztnu - ztno) ) * tmask(ji,jj,jk) &211 & / e3t(ji,jj,jk,Kbb)210 ztno = tr(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 211 tr(ji,jj,jk,jn,Kbb_a) = ( ztb + atfp * ( ztnu - ztno) ) & 212 & * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 212 213 ENDIF 213 214 ENDDO … … 221 222 DO ji=i1,i2 222 223 IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN 223 tr(ji,jj,jk,jn,Kmm ) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm)224 tr(ji,jj,jk,jn,Kmm_a) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm_a) 224 225 END IF 225 226 END DO … … 229 230 ! 230 231 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 231 tr(i1:i2,j1:j2,k1:k2,n1:n2,Kbb ) = tr(i1:i2,j1:j2,k1:k2,n1:n2,Kmm)232 tr(i1:i2,j1:j2,k1:k2,n1:n2,Kbb_a) = tr(i1:i2,j1:j2,k1:k2,n1:n2,Kmm_a) 232 233 ENDIF 233 234 ! -
NEMO/branches/2019/dev_r11943_MERGE_2019/src/NST/agrif_user.F90
r11960 r12229 6 6 !! Software governed by the CeCILL license (see ./LICENSE) 7 7 !!---------------------------------------------------------------------- 8 SUBROUTINE agrif_user 9 END SUBROUTINE agrif_user 10 11 SUBROUTINE agrif_before_regridding 12 END SUBROUTINE agrif_before_regridding 13 14 SUBROUTINE Agrif_InitWorkspace 15 !!---------------------------------------------------------------------- 16 !! *** ROUTINE Agrif_InitWorkspace *** 17 !!---------------------------------------------------------------------- 18 USE par_oce 19 USE dom_oce 20 USE nemogcm 21 USE mppini 22 !! 23 IMPLICIT NONE 24 !!---------------------------------------------------------------------- 25 ! 26 IF( .NOT. Agrif_Root() ) THEN 27 ! no more static variables 28 !!$! JC: change to allow for different vertical levels 29 !!$! jpk is already set 30 !!$! keep it jpk possibly different from jpkglo which 31 !!$! hold parent grid vertical levels number (set earlier) 32 !!$! jpk = jpkglo 33 ENDIF 34 ! 35 END SUBROUTINE Agrif_InitWorkspace 36 37 38 SUBROUTINE Agrif_InitValues 8 SUBROUTINE agrif_user 9 END SUBROUTINE agrif_user 10 11 SUBROUTINE agrif_before_regridding 12 END SUBROUTINE agrif_before_regridding 13 14 SUBROUTINE Agrif_InitWorkspace 15 END SUBROUTINE Agrif_InitWorkspace 16 17 SUBROUTINE Agrif_InitValues 39 18 !!---------------------------------------------------------------------- 40 19 !! *** ROUTINE Agrif_InitValues *** 41 !! 42 !! ** Purpose :: Declaration of variables to be interpolated 43 !!---------------------------------------------------------------------- 44 USE Agrif_Util 45 USE oce 46 USE dom_oce 47 USE nemogcm 48 USE tradmp 49 USE bdy_oce , ONLY: ln_bdy 50 !! 51 IMPLICIT NONE 52 !!---------------------------------------------------------------------- 53 ! 54 CALL nemo_init !* Initializations of each fine grid 55 Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs ! agrif_oce module copies of time level indices 56 57 ! !* Agrif initialization 58 CALL agrif_nemo_init 59 CALL Agrif_InitValues_cont_dom 60 CALL Agrif_InitValues_cont 20 !!---------------------------------------------------------------------- 21 USE nemogcm 22 !!---------------------------------------------------------------------- 23 ! 24 CALL nemo_init !* Initializations of each fine grid 25 Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs ! agrif_oce module copies of time level indices 26 ! 27 ! !* Agrif initialization 28 CALL agrif_nemo_init 29 CALL Agrif_InitValues_cont_dom 30 CALL Agrif_InitValues_cont 61 31 # if defined key_top 62 CALL Agrif_InitValues_cont_top32 CALL Agrif_InitValues_cont_top 63 33 # endif 64 34 # if defined key_si3 65 CALL Agrif_InitValues_cont_ice 66 # endif 67 ! 68 END SUBROUTINE Agrif_initvalues 69 70 71 SUBROUTINE Agrif_InitValues_cont_dom 72 !!---------------------------------------------------------------------- 73 !! *** ROUTINE Agrif_InitValues_cont *** 74 !! 75 !! ** Purpose :: Declaration of variables to be interpolated 76 !!---------------------------------------------------------------------- 77 USE Agrif_Util 78 USE oce 79 USE dom_oce 80 USE nemogcm 81 USE in_out_manager 82 USE agrif_oce_update 83 USE agrif_oce_interp 84 USE agrif_oce_sponge 85 ! 86 IMPLICIT NONE 87 !!---------------------------------------------------------------------- 88 ! 89 ! Declaration of the type of variable which have to be interpolated 90 ! 91 CALL agrif_declare_var_dom 92 ! 93 END SUBROUTINE Agrif_InitValues_cont_dom 94 95 96 SUBROUTINE agrif_declare_var_dom 97 !!---------------------------------------------------------------------- 98 !! *** ROUTINE agrif_declare_var *** 99 !! 100 !! ** Purpose :: Declaration of variables to be interpolated 101 !!---------------------------------------------------------------------- 102 USE agrif_util 103 USE par_oce 104 USE oce 105 ! 106 IMPLICIT NONE 107 ! 108 INTEGER :: ind1, ind2, ind3 35 CALL Agrif_InitValues_cont_ice 36 # endif 37 ! 38 END SUBROUTINE Agrif_initvalues 39 40 SUBROUTINE Agrif_InitValues_cont_dom 41 !!---------------------------------------------------------------------- 42 !! *** ROUTINE Agrif_InitValues_cont_dom *** 43 !!---------------------------------------------------------------------- 44 ! 45 CALL agrif_declare_var_dom 46 ! 47 END SUBROUTINE Agrif_InitValues_cont_dom 48 49 SUBROUTINE agrif_declare_var_dom 50 !!---------------------------------------------------------------------- 51 !! *** ROUTINE agrif_declare_var_dom *** 52 !!---------------------------------------------------------------------- 53 USE par_oce, ONLY: nbghostcells 54 ! 55 IMPLICIT NONE 56 ! 57 INTEGER :: ind1, ind2, ind3 109 58 !!---------------------------------------------------------------------- 110 59 111 60 ! 1. Declaration of the type of variable which have to be interpolated 112 61 !--------------------------------------------------------------------- 113 ind1 = nbghostcells114 ind2 = 1 + nbghostcells115 ind3 = 2 + nbghostcells116 CALL agrif_declare_variable((/1,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),e1u_id)117 CALL agrif_declare_variable((/2,1/),(/ind3,ind2/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),e2v_id)62 ind1 = nbghostcells 63 ind2 = 1 + nbghostcells 64 ind3 = 2 + nbghostcells 65 CALL agrif_declare_variable((/1,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),e1u_id) 66 CALL agrif_declare_variable((/2,1/),(/ind3,ind2/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),e2v_id) 118 67 119 68 ! 2. Type of interpolation 120 69 !------------------------- 121 CALL Agrif_Set_bcinterp( e1u_id, interp1=Agrif_linear, interp2=AGRIF_ppm )122 CALL Agrif_Set_bcinterp( e2v_id, interp1=AGRIF_ppm , interp2=Agrif_linear )70 CALL Agrif_Set_bcinterp( e1u_id, interp1=Agrif_linear, interp2=AGRIF_ppm ) 71 CALL Agrif_Set_bcinterp( e2v_id, interp1=AGRIF_ppm , interp2=Agrif_linear ) 123 72 124 73 ! 3. Location of interpolation 125 74 !----------------------------- 126 CALL Agrif_Set_bc(e1u_id,(/0,ind1-1/))127 CALL Agrif_Set_bc(e2v_id,(/0,ind1-1/))75 CALL Agrif_Set_bc(e1u_id,(/0,ind1-1/)) 76 CALL Agrif_Set_bc(e2v_id,(/0,ind1-1/)) 128 77 129 78 ! 4. Update type 130 79 !--------------- 131 80 # if defined UPD_HIGH 132 CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Average, update2=Agrif_Update_Full_Weighting)133 CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Full_Weighting, update2=Agrif_Update_Average)81 CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Average, update2=Agrif_Update_Full_Weighting) 82 CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Full_Weighting, update2=Agrif_Update_Average) 134 83 #else 135 CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Copy, update2=Agrif_Update_Average)136 CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Average, update2=Agrif_Update_Copy)84 CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Copy, update2=Agrif_Update_Average) 85 CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Average, update2=Agrif_Update_Copy) 137 86 #endif 138 87 139 END SUBROUTINE agrif_declare_var_dom 140 141 142 SUBROUTINE Agrif_InitValues_cont 88 END SUBROUTINE agrif_declare_var_dom 89 90 SUBROUTINE Agrif_InitValues_cont 143 91 !!---------------------------------------------------------------------- 144 92 !! *** ROUTINE Agrif_InitValues_cont *** 145 !! 146 !! ** Purpose :: Declaration of variables to be interpolated 147 !!---------------------------------------------------------------------- 148 USE agrif_oce_update 149 USE agrif_oce_interp 150 USE agrif_oce_sponge 151 USE Agrif_Util 152 USE oce 153 USE dom_oce 154 USE zdf_oce 155 USE nemogcm 156 ! 157 USE lib_mpp 158 USE in_out_manager 159 ! 160 IMPLICIT NONE 161 ! 162 LOGICAL :: check_namelist 163 CHARACTER(len=15) :: cl_check1, cl_check2, cl_check3, cl_check4 164 !!---------------------------------------------------------------------- 165 166 ! 1. Declaration of the type of variable which have to be interpolated 167 !--------------------------------------------------------------------- 168 CALL agrif_declare_var 169 170 ! 2. First interpolations of potentially non zero fields 171 !------------------------------------------------------- 172 Agrif_SpecialValue = 0._wp 173 Agrif_UseSpecialValue = .TRUE. 174 CALL Agrif_Bc_variable(tsn_id,calledweight=1.,procname=interptsn) 175 CALL Agrif_Sponge 176 tabspongedone_tsn = .FALSE. 177 CALL Agrif_Bc_variable(tsn_sponge_id,calledweight=1.,procname=interptsn_sponge) 178 ! reset ts(:,:,:,:,Krhs_a) to zero 179 ts(:,:,:,:,Krhs_a) = 0. 180 181 Agrif_UseSpecialValue = ln_spc_dyn 182 CALL Agrif_Bc_variable(un_interp_id,calledweight=1.,procname=interpun) 183 CALL Agrif_Bc_variable(vn_interp_id,calledweight=1.,procname=interpvn) 184 tabspongedone_u = .FALSE. 185 tabspongedone_v = .FALSE. 186 CALL Agrif_Bc_variable(un_sponge_id,calledweight=1.,procname=interpun_sponge) 187 tabspongedone_u = .FALSE. 188 tabspongedone_v = .FALSE. 189 CALL Agrif_Bc_variable(vn_sponge_id,calledweight=1.,procname=interpvn_sponge) 190 191 Agrif_UseSpecialValue = .TRUE. 192 CALL Agrif_Bc_variable(sshn_id,calledweight=1., procname=interpsshn ) 193 hbdy_w(:,:) = 0.e0 ; hbdy_e(:,:) = 0.e0 ; hbdy_n(:,:) = 0.e0 ; hbdy_s(:,:) = 0.e0 194 ssh(:,:,Krhs_a) = 0.e0 195 196 IF ( ln_dynspg_ts ) THEN 93 !!---------------------------------------------------------------------- 94 USE agrif_oce 95 USE agrif_oce_interp 96 USE agrif_oce_sponge 97 USE dom_oce 98 USE oce 99 USE lib_mpp 100 USE lbclnk 101 ! 102 IMPLICIT NONE 103 ! 104 INTEGER :: ji, jj 105 LOGICAL :: check_namelist 106 CHARACTER(len=15) :: cl_check1, cl_check2, cl_check3, cl_check4 107 #if defined key_vertical 108 REAL(wp), DIMENSION(jpi,jpj) :: zk ! workspace 109 #endif 110 !!---------------------------------------------------------------------- 111 112 ! 1. Declaration of the type of variable which have to be interpolated 113 !--------------------------------------------------------------------- 114 CALL agrif_declare_var 115 116 ! 2. First interpolations of potentially non zero fields 117 !------------------------------------------------------- 118 119 #if defined key_vertical 120 ! Build consistent parent bathymetry and number of levels 121 ! on the child grid 122 Agrif_UseSpecialValue = .FALSE. 123 ht0_parent(:,:) = 0._wp 124 mbkt_parent(:,:) = 0 125 ! 126 CALL Agrif_Bc_variable(ht0_id ,calledweight=1.,procname=interpht0 ) 127 CALL Agrif_Bc_variable(mbkt_id,calledweight=1.,procname=interpmbkt) 128 ! 129 ! Assume step wise change of bathymetry near interface 130 ! TODO: Switch to linear interpolation of bathymetry in the s-coordinate case 131 ! and no refinement 132 DO jj = 1, jpjm1 133 DO ji = 1, jpim1 134 mbku_parent(ji,jj) = MIN( mbkt_parent(ji+1,jj ) , mbkt_parent(ji,jj) ) 135 mbkv_parent(ji,jj) = MIN( mbkt_parent(ji ,jj+1) , mbkt_parent(ji,jj) ) 136 END DO 137 END DO 138 IF ( ln_sco.AND.Agrif_Parent(ln_sco) ) THEN 139 DO jj = 1, jpjm1 140 DO ji = 1, jpim1 141 hu0_parent(ji,jj) = 0.5_wp * ( ht0_parent(ji,jj)+ht0_parent(ji+1,jj) ) 142 hv0_parent(ji,jj) = 0.5_wp * ( ht0_parent(ji,jj)+ht0_parent(ji,jj+1) ) 143 END DO 144 END DO 145 ELSE 146 DO jj = 1, jpjm1 147 DO ji = 1, jpim1 148 hu0_parent(ji,jj) = MIN( ht0_parent(ji,jj), ht0_parent(ji+1,jj)) 149 hv0_parent(ji,jj) = MIN( ht0_parent(ji,jj), ht0_parent(ji,jj+1)) 150 END DO 151 END DO 152 153 ENDIF 154 ! 155 CALL lbc_lnk( 'Agrif_InitValues_cont', hu0_parent, 'U', 1. ) 156 CALL lbc_lnk( 'Agrif_InitValues_cont', hv0_parent, 'V', 1. ) 157 zk(:,:) = REAL( mbku_parent(:,:), wp ) ; CALL lbc_lnk( 'Agrif_InitValues_cont', zk, 'U', 1. ) 158 mbku_parent(:,:) = MAX( NINT( zk(:,:) ), 1 ) 159 zk(:,:) = REAL( mbkv_parent(:,:), wp ) ; CALL lbc_lnk( 'Agrif_InitValues_cont', zk, 'V', 1. ) 160 mbkv_parent(:,:) = MAX( NINT( zk(:,:) ), 1 ) 161 #endif 162 163 Agrif_SpecialValue = 0._wp 164 Agrif_UseSpecialValue = .TRUE. 165 CALL Agrif_Bc_variable(tsn_id,calledweight=1.,procname=interptsn) 166 CALL Agrif_Sponge 167 tabspongedone_tsn = .FALSE. 168