Changeset 12191
- Timestamp:
- 2019-12-11T16:56:06+01:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019
- Files:
-
- 13 edited
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/cfgs/SHARED/field_def_nemo-oce.xml
r12184 r12191 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) --> … … 489 491 <field id="uoces" long_name="ocean transport along i-axis times salinity (CRS)" unit="1e-3*m/s" grid_ref="grid_U_3D" /> 490 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> 491 493 <!-- AGRIF sponge --> 494 <field id="agrif_spu" long_name=" AGRIF u-sponge coefficient" unit=" " /> 492 495 <!-- u-eddy diffusivity coefficients (available if ln_traldf_OFF=F) --> 493 496 <field id="ahtu_2d" long_name=" surface u-eddy diffusivity coefficient" unit="m2/s or m4/s" /> … … 547 550 <field id="voces" long_name="ocean transport along j-axis times salinity (CRS)" unit="1e-3*m/s" grid_ref="grid_V_3D" /> 548 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> 549 552 <!-- AGRIF sponge --> 553 <field id="agrif_spv" long_name=" AGRIF v-sponge coefficient" unit=" " /> 550 554 <!-- v-eddy diffusivity coefficients (available if ln_traldf_OFF=F) --> 551 555 <field id="ahtv_2d" long_name=" surface v-eddy diffusivity coefficient" unit="m2/s or (m4/s)^1/2" /> … … 636 640 637 641 <!-- F grid --> 642 <!-- AGRIF sponge --> 643 <field id="agrif_spf" long_name=" AGRIF f-sponge coefficient" unit=" " /> 638 644 <!-- f-eddy viscosity coefficients (ldfdyn) --> 639 645 <field id="ahmf_2d" long_name=" surface f-eddy viscosity coefficient" unit="m2/s or m4/s" /> -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/cfgs/SHARED/namelist_ref
r12184 r12191 637 637 &namagrif ! AGRIF zoom ("key_agrif") 638 638 !----------------------------------------------------------------------- 639 ln_agrif_2way = .true. ! activate two way nesting 639 640 ln_spc_dyn = .true. ! use 0 as special value for dynamics 640 641 rn_sponge_tra = 2880. ! coefficient for tracer sponge layer [m2/s] 641 642 rn_sponge_dyn = 2880. ! coefficient for dynamics sponge layer [m2/s] 643 rn_trelax_tra = 0.01 ! inverse of relaxation time (in steps) for tracers [] 644 rn_trelax_dyn = 0.01 ! inverse of relaxation time (in steps) for dynamics [] 642 645 ln_chk_bathy = .false. ! =T check the parent bathymetry 643 646 / -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/NST/agrif_all_update.F90
r10069 r12191 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_r12072_MERGE_OPTION2_2019/src/NST/agrif_ice_update.F90
r10069 r12191 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_r12072_MERGE_OPTION2_2019/src/NST/agrif_oce.F90
r10425 r12191 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 50 # if defined key_vertical 51 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht0_parent, hu0_parent, hv0_parent 52 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt_parent, mbku_parent, mbkv_parent 53 # endif 53 54 54 55 INTEGER, PUBLIC :: tsn_id ! AGRIF profile for tracers interpolation and update … … 64 65 INTEGER, PUBLIC :: scales_t_id 65 66 INTEGER, PUBLIC :: avt_id, avm_id, en_id ! TKE related identificators 66 INTEGER, PUBLIC :: umsk_id, vmsk_id67 INTEGER, PUBLIC :: mbkt_id, ht0_id 67 68 INTEGER, PUBLIC :: kindic_agr 68 69 … … 82 83 ierr(:) = 0 83 84 ! 84 ALLOCATE( fsaht_spu(jpi,jpj), fsaht_spv(jpi,jpj), & 85 & fsahm_spt(jpi,jpj), fsahm_spf(jpi,jpj), & 86 & tabspongedone_tsn(jpi,jpj), & 85 ALLOCATE( fspu(jpi,jpj), fspv(jpi,jpj), & 86 & fspt(jpi,jpj), fspf(jpi,jpj), & 87 & tabspongedone_tsn(jpi,jpj), & 88 & utint_stage(jpi,jpj), vtint_stage(jpi,jpj), & 87 89 # if defined key_top 88 90 & tabspongedone_trn(jpi,jpj), & 89 # endif 91 # endif 92 # if defined key_vertical 93 & ht0_parent(jpi,jpj), mbkt_parent(jpi,jpj), & 94 & hu0_parent(jpi,jpj), mbku_parent(jpi,jpj), & 95 & hv0_parent(jpi,jpj), mbkv_parent(jpi,jpj), & 96 # endif 90 97 & tabspongedone_u (jpi,jpj), & 91 98 & tabspongedone_v (jpi,jpj), STAT = ierr(1) ) 92 99 93 ALLOCATE( ubdy_w(nbghostcells,jpj), vbdy_w(nbghostcells,jpj), hbdy_w(nbghostcells,jpj), & 94 & ubdy_e(nbghostcells,jpj), vbdy_e(nbghostcells,jpj), hbdy_e(nbghostcells,jpj), & 95 & ubdy_n(jpi,nbghostcells), vbdy_n(jpi,nbghostcells), hbdy_n(jpi,nbghostcells), & 96 & ubdy_s(jpi,nbghostcells), vbdy_s(jpi,nbghostcells), hbdy_s(jpi,nbghostcells), STAT = ierr(2) ) 100 ALLOCATE( ubdy(jpi,jpj), vbdy(jpi,jpj), hbdy(jpi,jpj), STAT = ierr(2) ) 97 101 98 102 agrif_oce_alloc = MAXVAL(ierr) … … 100 104 END FUNCTION agrif_oce_alloc 101 105 102 #if defined key_vertical103 SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout)104 !!----------------------------------------------------------------------105 !! *** FUNCTION reconstructandremap ***106 !!----------------------------------------------------------------------107 IMPLICIT NONE108 INTEGER N, Nout109 REAL(wp) tabin(N), tabout(Nout)110 REAL(wp) hin(N), hout(Nout)111 REAL(wp) coeffremap(N,3),zwork(N,3)112 REAL(wp) zwork2(N+1,3)113 INTEGER jk114 DOUBLE PRECISION, PARAMETER :: dsmll=1.0d-8115 REAL(wp) q,q01,q02,q001,q002,q0116 REAL(wp) z_win(1:N+1), z_wout(1:Nout+1)117 REAL(wp),PARAMETER :: dpthin = 1.D-3118 INTEGER :: k1, kbox, ktop, ka, kbot119 REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop120 121 z_win(1)=0.; z_wout(1)= 0.122 DO jk=1,N123 z_win(jk+1)=z_win(jk)+hin(jk)124 ENDDO125 126 DO jk=1,Nout127 z_wout(jk+1)=z_wout(jk)+hout(jk)128 ENDDO129 130 DO jk=2,N131 zwork(jk,1)=1./(hin(jk-1)+hin(jk))132 ENDDO133 134 DO jk=2,N-1135 q0 = 1./(hin(jk-1)+hin(jk)+hin(jk+1))136 zwork(jk,2)=hin(jk-1)+2.*hin(jk)+hin(jk+1)137 zwork(jk,3)=q0138 ENDDO139 140 DO jk= 2,N141 zwork2(jk,1)=zwork(jk,1)*(tabin(jk)-tabin(jk-1))142 ENDDO143 144 coeffremap(:,1) = tabin(:)145 146 DO jk=2,N-1147 q001 = hin(jk)*zwork2(jk+1,1)148 q002 = hin(jk)*zwork2(jk,1)149 IF (q001*q002 < 0) then150 q001 = 0.151 q002 = 0.152 ENDIF153 q=zwork(jk,2)154 q01=q*zwork2(jk+1,1)155 q02=q*zwork2(jk,1)156 IF (abs(q001) > abs(q02)) q001 = q02157 IF (abs(q002) > abs(q01)) q002 = q01158 159 q=(q001-q002)*zwork(jk,3)160 q001=q001-q*hin(jk+1)161 q002=q002+q*hin(jk-1)162 163 coeffremap(jk,3)=coeffremap(jk,1)+q001164 coeffremap(jk,2)=coeffremap(jk,1)-q002165 166 zwork2(jk,1)=(2.*q001-q002)**2167 zwork2(jk,2)=(2.*q002-q001)**2168 ENDDO169 170 DO jk=1,N171 IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN172 coeffremap(jk,3) = coeffremap(jk,1)173 coeffremap(jk,2) = coeffremap(jk,1)174 zwork2(jk,1) = 0.175 zwork2(jk,2) = 0.176 ENDIF177 ENDDO178 179 DO jk=2,N180 q002=max(zwork2(jk-1,2),dsmll)181 q001=max(zwork2(jk,1),dsmll)182 zwork2(jk,3)=(q001*coeffremap(jk-1,3)+q002*coeffremap(jk,2))/(q001+q002)183 ENDDO184 185 zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3)186 zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3)187 188 DO jk=1,N189 q01=zwork2(jk+1,3)-coeffremap(jk,1)190 q02=coeffremap(jk,1)-zwork2(jk,3)191 q001=2.*q01192 q002=2.*q02193 IF (q01*q02<0) then194 q01=0.195 q02=0.196 ELSEIF (abs(q01)>abs(q002)) then197 q01=q002198 ELSEIF (abs(q02)>abs(q001)) then199 q02=q001200 ENDIF201 coeffremap(jk,2)=coeffremap(jk,1)-q02202 coeffremap(jk,3)=coeffremap(jk,1)+q01203 ENDDO204 205 zbot=0.0206 kbot=1207 DO jk=1,Nout208 ztop=zbot !top is bottom of previous layer209 ktop=kbot210 IF (ztop.GE.z_win(ktop+1)) then211 ktop=ktop+1212 ENDIF213 214 zbot=z_wout(jk+1)215 zthk=zbot-ztop216 217 IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(Nout+1)) THEN218 219 kbot=ktop220 DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.N)221 kbot=kbot+1222 ENDDO223 zbox=zbot224 DO k1= jk+1,Nout225 IF (z_wout(k1+1)-z_wout(k1).GT.dpthin) THEN226 exit !thick layer227 ELSE228 zbox=z_wout(k1+1) !include thin adjacent layers229 IF(zbox.EQ.z_wout(Nout+1)) THEN230 exit !at bottom231 ENDIF232 ENDIF233 ENDDO234 zthk=zbox-ztop235 236 kbox=ktop237 DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.N)238 kbox=kbox+1239 ENDDO240 241 IF(ktop.EQ.kbox) THEN242 IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN243 IF(hin(kbox).GT.dpthin) THEN244 q001 = (zbox-z_win(kbox))/hin(kbox)245 q002 = (ztop-z_win(kbox))/hin(kbox)246 q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002)247 q02=q01-1.+(q001+q002)248 q0=1.-q01-q02249 ELSE250 q0 = 1.0251 q01 = 0.252 q02 = 0.253 ENDIF254 tabout(jk)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3)255 ELSE256 tabout(jk) = tabin(kbox)257 ENDIF258 ELSE259 IF(ktop.LE.jk .AND. kbox.GE.jk) THEN260 ka = jk261 ELSEIF (kbox-ktop.GE.3) THEN262 ka = (kbox+ktop)/2263 ELSEIF (hin(ktop).GE.hin(kbox)) THEN264 ka = ktop265 ELSE266 ka = kbox267 ENDIF !choose ka268 269 offset=coeffremap(ka,1)270 271 qtop = z_win(ktop+1)-ztop !partial layer thickness272 IF(hin(ktop).GT.dpthin) THEN273 q=(ztop-z_win(ktop))/hin(ktop)274 q01=q*(q-1.)275 q02=q01+q276 q0=1-q01-q02277 ELSE278 q0 = 1.279 q01 = 0.280 q02 = 0.281 ENDIF282 283 tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop284 285 DO k1= ktop+1,kbox-1286 tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1)287 ENDDO !k1288 289 qbot = zbox-z_win(kbox) !partial layer thickness290 IF(hin(kbox).GT.dpthin) THEN291 q=qbot/hin(kbox)292 q01=(q-1.)**2293 q02=q01-1.+q294 q0=1-q01-q02295 ELSE296 q0 = 1.0297 q01 = 0.298 q02 = 0.299 ENDIF300 301 tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot302 303 rpsum=1.0d0/zthk304 tabout(jk)=offset+tsum*rpsum305 306 ENDIF !single or multiple layers307 ELSE308 IF (jk==1) THEN309 write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1)310 ENDIF311 tabout(jk) = tabout(jk-1)312 313 ENDIF !normal:thin layer314 ENDDO !jk315 316 return317 end subroutine reconstructandremap318 #endif319 320 106 #endif 321 107 !!====================================================================== -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/NST/agrif_oce_interp.F90
r10068 r12191 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 ua_b(ibdy1:ibdy2,:) = 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 ua_b(ji,:) = 0._wp 104 110 105 DO jk = 1, jpkm1 111 106 DO jj = 1, jpj 112 ua_b( ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) &113 & + e3u_a(ibdy1:ibdy2,jj,jk) * ua(ibdy1:ibdy2,jj,jk) * umask(ibdy1:ibdy2,jj,jk)114 115 END DO 107 ua_b(ji,jj) = ua_b(ji,jj) + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 108 END DO 109 END DO 110 116 111 DO jj = 1, jpj 117 ua_b(ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj) 118 END DO 119 ENDIF 120 ! 121 IF( .NOT.lk_agrif_clp ) THEN 122 DO jk=1,jpkm1 ! Smooth 123 DO jj=j1,j2 124 ua(ibdy2,jj,jk) = 0.25_wp*(ua(ibdy2-1,jj,jk)+2._wp*ua(ibdy2,jj,jk)+ua(ibdy2+1,jj,jk)) 125 END DO 126 END DO 127 ENDIF 128 ! 129 zub(ibdy1:ibdy2,:) = 0._wp ! Correct transport 112 ua_b(ji,jj) = ua_b(ji,jj) * r1_hu_a(ji,jj) 113 END DO 114 END DO 115 ENDIF 116 ! 117 DO ji = mi0(ibdy1), mi1(ibdy2) 118 zub(ji,:) = 0._wp ! Correct transport 130 119 DO jk = 1, jpkm1 131 120 DO jj = 1, jpj 132 zub( ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) &133 & + e3u_a( ibdy1:ibdy2,jj,jk) * ua(ibdy1:ibdy2,jj,jk)*umask(ibdy1:ibdy2,jj,jk)121 zub(ji,jj) = zub(ji,jj) & 122 & + e3u_a(ji,jj,jk) * ua(ji,jj,jk)*umask(ji,jj,jk) 134 123 END DO 135 124 END DO 136 125 DO jj=1,jpj 137 zub( ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj)126 zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj) 138 127 END DO 139 128 140 129 DO jk = 1, jpkm1 141 130 DO jj = 1, jpj 142 ua( ibdy1:ibdy2,jj,jk) = ( ua(ibdy1:ibdy2,jj,jk) &143 & + ua_b(ibdy1:ibdy2,jj)-zub(ibdy1:ibdy2,jj)) * umask(ibdy1:ibdy2,jj,jk)144 145 131 ua(ji,jj,jk) = ( ua(ji,jj,jk) + ua_b(ji,jj)-zub(ji,jj)) * umask(ji,jj,jk) 132 END DO 133 END DO 134 END DO 146 135 147 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 148 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 149 139 DO jk = 1, jpkm1 150 140 DO jj = 1, jpj 151 zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) & 152 & + e3v_a(ibdy1:ibdy2,jj,jk) * va(ibdy1:ibdy2,jj,jk) * vmask(ibdy1:ibdy2,jj,jk) 141 zvb(ji,jj) = zvb(ji,jj) + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 153 142 END DO 154 143 END DO 155 144 DO jj = 1, jpj 156 zvb( ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv_a(ibdy1:ibdy2,jj)145 zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj) 157 146 END DO 158 147 DO jk = 1, jpkm1 159 148 DO jj = 1, jpj 160 va(ibdy1:ibdy2,jj,jk) = ( va(ibdy1:ibdy2,jj,jk) & 161 & + va_b(ibdy1:ibdy2,jj)-zvb(ibdy1:ibdy2,jj))*vmask(ibdy1:ibdy2,jj,jk) 162 END DO 163 END DO 164 ENDIF 165 ! 166 DO jk = 1, jpkm1 ! Mask domain edges 167 DO jj = 1, jpj 168 ua(1,jj,jk) = 0._wp 169 va(1,jj,jk) = 0._wp 170 END DO 171 END DO 149 va(ji,jj,jk) = ( va(ji,jj,jk) + va_b(ji,jj)-zvb(ji,jj))*vmask(ji,jj,jk) 150 END DO 151 END DO 152 END DO 172 153 ENDIF 173 154 174 155 ! --- East --- ! 175 IF( nbondi == 1 .OR. nbondi == 2 ) THEN176 ibdy1 = nlci-1-nbghostcells177 ibdy2 = nlci-2178 !179 IF( .NOT.ln_dynspg_ts ) THEN ! Store transport180 ua_b( ibdy1:ibdy2,:) = 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 ua_b(ji,:) = 0._wp 181 162 DO jk = 1, jpkm1 182 163 DO jj = 1, jpj 183 ua_b( ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) &184 & + e3u_a( ibdy1:ibdy2,jj,jk) * ua(ibdy1:ibdy2,jj,jk) * umask(ibdy1:ibdy2,jj,jk)164 ua_b(ji,jj) = ua_b(ji,jj) & 165 & + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 185 166 END DO 186 167 END DO 187 168 DO jj = 1, jpj 188 ua_b(ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj) 189 END DO 190 ENDIF 191 ! 192 IF( .NOT.lk_agrif_clp ) THEN 193 DO jk=1,jpkm1 ! Smooth 194 DO jj=j1,j2 195 ua(ibdy1,jj,jk) = 0.25_wp*(ua(ibdy1-1,jj,jk)+2._wp*ua(ibdy1,jj,jk)+ua(ibdy1+1,jj,jk)) 196 END DO 197 END DO 198 ENDIF 199 ! 200 zub(ibdy1:ibdy2,:) = 0._wp ! Correct transport 169 ua_b(ji,jj) = ua_b(ji,jj) * r1_hu_a(ji,jj) 170 END DO 171 END DO 172 ENDIF 173 ! 174 DO ji = mi0(ibdy1), mi1(ibdy2) 175 zub(ji,:) = 0._wp ! Correct transport 201 176 DO jk = 1, jpkm1 202 177 DO jj = 1, jpj 203 zub( ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) &204 & + e3u_a( ibdy1:ibdy2,jj,jk) * ua(ibdy1:ibdy2,jj,jk) * umask(ibdy1:ibdy2,jj,jk)178 zub(ji,jj) = zub(ji,jj) & 179 & + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 205 180 END DO 206 181 END DO 207 182 DO jj=1,jpj 208 zub( ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj)183 zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj) 209 184 END DO 210 185 211 186 DO jk = 1, jpkm1 212 187 DO jj = 1, jpj 213 ua(ibdy1:ibdy2,jj,jk) = ( ua(ibdy1:ibdy2,jj,jk) & 214 & + ua_b(ibdy1:ibdy2,jj)-zub(ibdy1:ibdy2,jj))*umask(ibdy1:ibdy2,jj,jk) 215 END DO 216 END DO 188 ua(ji,jj,jk) = ( ua(ji,jj,jk) & 189 & + ua_b(ji,jj)-zub(ji,jj))*umask(ji,jj,jk) 190 END DO 191 END DO 192 END DO 217 193 218 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 219 ibdy1 = ibdy1 + 1 220 ibdy2 = ibdy2 + 1 221 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 222 199 DO jk = 1, jpkm1 223 200 DO jj = 1, jpj 224 zvb( ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) &225 & + e3v_a( ibdy1:ibdy2,jj,jk) * va(ibdy1:ibdy2,jj,jk) * vmask(ibdy1:ibdy2,jj,jk)201 zvb(ji,jj) = zvb(ji,jj) & 202 & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 226 203 END DO 227 204 END DO 228 205 DO jj = 1, jpj 229 zvb( ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv_a(ibdy1:ibdy2,jj)206 zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj) 230 207 END DO 231 208 DO jk = 1, jpkm1 232 209 DO jj = 1, jpj 233 va(ibdy1:ibdy2,jj,jk) = ( va(ibdy1:ibdy2,jj,jk) & 234 & + va_b(ibdy1:ibdy2,jj)-zvb(ibdy1:ibdy2,jj)) * vmask(ibdy1:ibdy2,jj,jk) 235 END DO 236 END DO 237 ENDIF 238 ! 239 DO jk = 1, jpkm1 ! Mask domain edges 240 DO jj = 1, jpj 241 ua(nlci-1,jj,jk) = 0._wp 242 va(nlci ,jj,jk) = 0._wp 243 END DO 244 END DO 210 va(ji,jj,jk) = ( va(ji,jj,jk) & 211 & + va_b(ji,jj)-zvb(ji,jj)) * vmask(ji,jj,jk) 212 END DO 213 END DO 214 END DO 245 215 ENDIF 246 216 247 217 ! --- South --- ! 248 IF( nbondj == -1 .OR. nbondj == 2 ) THEN249 jbdy1 = 2250 jbdy2 = 1+nbghostcells251 !252 IF( .NOT.ln_dynspg_ts ) THEN ! Store transport253 va_b(:,j bdy1:jbdy2) = 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 va_b(:,jj) = 0._wp 254 224 DO jk = 1, jpkm1 255 225 DO ji = 1, jpi 256 va_b(ji,j bdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) &257 & + e3v_a(ji,j bdy1:jbdy2,jk) * va(ji,jbdy1:jbdy2,jk) * vmask(ji,jbdy1:jbdy2,jk)226 va_b(ji,jj) = va_b(ji,jj) & 227 & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 258 228 END DO 259 229 END DO 260 230 DO ji=1,jpi 261 va_b(ji,jbdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2) 262 END DO 263 ENDIF 264 ! 265 IF ( .NOT.lk_agrif_clp ) THEN 266 DO jk = 1, jpkm1 ! Smooth 267 DO ji = i1, i2 268 va(ji,jbdy2,jk) = 0.25_wp*(va(ji,jbdy2-1,jk)+2._wp*va(ji,jbdy2,jk)+va(ji,jbdy2+1,jk)) 269 END DO 270 END DO 271 ENDIF 272 ! 273 zvb(:,jbdy1:jbdy2) = 0._wp ! Correct transport 231 va_b(ji,jj) = va_b(ji,jj) * r1_hv_a(ji,jj) 232 END DO 233 END DO 234 ENDIF 235 ! 236 DO jj = mj0(jbdy1), mj1(jbdy2) 237 zvb(:,jj) = 0._wp ! Correct transport 274 238 DO jk=1,jpkm1 275 239 DO ji=1,jpi 276 zvb(ji,j bdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) &277 & + e3v_a(ji,j bdy1:jbdy2,jk) * va(ji,jbdy1:jbdy2,jk) * vmask(ji,jbdy1:jbdy2,jk)240 zvb(ji,jj) = zvb(ji,jj) & 241 & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 278 242 END DO 279 243 END DO 280 244 DO ji = 1, jpi 281 zvb(ji,j bdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2)245 zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj) 282 246 END DO 283 247 284 248 DO jk = 1, jpkm1 285 249 DO ji = 1, jpi 286 va(ji,jbdy1:jbdy2,jk) = ( va(ji,jbdy1:jbdy2,jk) & 287 & + va_b(ji,jbdy1:jbdy2) - zvb(ji,jbdy1:jbdy2) ) * vmask(ji,jbdy1:jbdy2,jk) 288 END DO 289 END DO 250 va(ji,jj,jk) = ( va(ji,jj,jk) & 251 & + va_b(ji,jj) - zvb(ji,jj) ) * vmask(ji,jj,jk) 252 END DO 253 END DO 254 END DO 290 255 291 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 292 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 293 259 DO jk = 1, jpkm1 294 260 DO ji = 1, jpi 295 zub(ji,j bdy1:jbdy2) = zub(ji,jbdy1:jbdy2) &296 & + e3u_a(ji,j bdy1:jbdy2,jk) * ua(ji,jbdy1:jbdy2,jk) * umask(ji,jbdy1:jbdy2,jk)261 zub(ji,jj) = zub(ji,jj) & 262 & + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 297 263 END DO 298 264 END DO 299 265 DO ji = 1, jpi 300 zub(ji,j bdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu_a(ji,jbdy1:jbdy2)266 zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj) 301 267 END DO 302 268 303 269 DO jk = 1, jpkm1 304 270 DO ji = 1, jpi 305 ua(ji,jbdy1:jbdy2,jk) = ( ua(ji,jbdy1:jbdy2,jk) & 306 & + ua_b(ji,jbdy1:jbdy2) - zub(ji,jbdy1:jbdy2) ) * umask(ji,jbdy1:jbdy2,jk) 307 END DO 308 END DO 309 ENDIF 310 ! 311 DO jk = 1, jpkm1 ! Mask domain edges 312 DO ji = 1, jpi 313 ua(ji,1,jk) = 0._wp 314 va(ji,1,jk) = 0._wp 315 END DO 316 END DO 271 ua(ji,jj,jk) = ( ua(ji,jj,jk) & 272 & + ua_b(ji,jj) - zub(ji,jj) ) * umask(ji,jj,jk) 273 END DO 274 END DO 275 END DO 317 276 ENDIF 318 277 319 278 ! --- North --- ! 320 IF( nbondj == 1 .OR. nbondj == 2 ) THEN321 jbdy1 = nlcj-1-nbghostcells322 jbdy2 = nlcj-2323 !324 IF( .NOT.ln_dynspg_ts ) THEN ! Store transport325 va_b(:,j bdy1:jbdy2) = 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 va_b(:,jj) = 0._wp 326 285 DO jk = 1, jpkm1 327 286 DO ji = 1, jpi 328 va_b(ji,j bdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) &329 & + e3v_a(ji,j bdy1:jbdy2,jk) * va(ji,jbdy1:jbdy2,jk) * vmask(ji,jbdy1:jbdy2,jk)287 va_b(ji,jj) = va_b(ji,jj) & 288 & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 330 289 END DO 331 290 END DO 332 291 DO ji=1,jpi 333 va_b(ji,jbdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2) 334 END DO 335 ENDIF 336 ! 337 IF ( .NOT.lk_agrif_clp ) THEN 338 DO jk = 1, jpkm1 ! Smooth 339 DO ji = i1, i2 340 va(ji,jbdy1,jk) = 0.25_wp*(va(ji,jbdy1-1,jk)+2._wp*va(ji,jbdy1,jk)+va(ji,jbdy1+1,jk)) 341 END DO 342 END DO 343 ENDIF 344 ! 345 zvb(:,jbdy1:jbdy2) = 0._wp ! Correct transport 292 va_b(ji,jj) = va_b(ji,jj) * r1_hv_a(ji,jj) 293 END DO 294 END DO 295 ENDIF 296 ! 297 DO jj = mj0(jbdy1), mj1(jbdy2) 298 zvb(:,jj) = 0._wp ! Correct transport 346 299 DO jk=1,jpkm1 347 300 DO ji=1,jpi 348 zvb(ji,j bdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) &349 & + e3v_a(ji,j bdy1:jbdy2,jk) * va(ji,jbdy1:jbdy2,jk) * vmask(ji,jbdy1:jbdy2,jk)301 zvb(ji,jj) = zvb(ji,jj) & 302 & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 350 303 END DO 351 304 END DO 352 305 DO ji = 1, jpi 353 zvb(ji,j bdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2)306 zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj) 354 307 END DO 355 308 356 309 DO jk = 1, jpkm1 357 310 DO ji = 1, jpi 358 va(ji,jbdy1:jbdy2,jk) = ( va(ji,jbdy1:jbdy2,jk) & 359 & + va_b(ji,jbdy1:jbdy2) - zvb(ji,jbdy1:jbdy2) ) * vmask(ji,jbdy1:jbdy2,jk) 360 END DO 361 END DO 311 va(ji,jj,jk) = ( va(ji,jj,jk) & 312 & + va_b(ji,jj) - zvb(ji,jj) ) * vmask(ji,jj,jk) 313 END DO 314 END DO 315 END DO 362 316 363 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 364 jbdy1 = jbdy1 + 1 365 jbdy2 = jbdy2 + 1 366 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 367 322 DO jk = 1, jpkm1 368 323 DO ji = 1, jpi 369 zub(ji,j bdy1:jbdy2) = zub(ji,jbdy1:jbdy2) &370 & + e3u_a(ji,j bdy1:jbdy2,jk) * ua(ji,jbdy1:jbdy2,jk) * umask(ji,jbdy1:jbdy2,jk)324 zub(ji,jj) = zub(ji,jj) & 325 & + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 371 326 END DO 372 327 END DO 373 328 DO ji = 1, jpi 374 zub(ji,j bdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu_a(ji,jbdy1:jbdy2)329 zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj) 375 330 END DO 376 331 377 332 DO jk = 1, jpkm1 378 333 DO ji = 1, jpi 379 ua(ji,jbdy1:jbdy2,jk) = ( ua(ji,jbdy1:jbdy2,jk) & 380 & + ua_b(ji,jbdy1:jbdy2) - zub(ji,jbdy1:jbdy2) ) * umask(ji,jbdy1:jbdy2,jk) 381 END DO 382 END DO 383 ENDIF 384 ! 385 DO jk = 1, jpkm1 ! Mask domain edges 386 DO ji = 1, jpi 387 ua(ji,nlcj ,jk) = 0._wp 388 va(ji,nlcj-1,jk) = 0._wp 389 END DO 390 END DO 334 ua(ji,jj,jk) = ( ua(ji,jj,jk) & 335 & + ua_b(ji,jj) - zub(ji,jj) ) * umask(ji,jj,jk) 336 END DO 337 END DO 338 END DO 391 339 ENDIF 392 340 ! … … 401 349 !! 402 350 INTEGER :: ji, jj 351 INTEGER :: istart, iend, jstart, jend 403 352 !!---------------------------------------------------------------------- 404 353 ! 405 354 IF( Agrif_Root() ) RETURN 406 355 ! 407 IF((nbondi == -1).OR.(nbondi == 2)) THEN 356 !--- West ---! 357 istart = 2 358 iend = nbghostcells+1 359 DO ji = mi0(istart), mi1(iend) 408 360 DO jj=1,jpj 409 va_e(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * hvr_e(2:nbghostcells+1,jj) 410 ! Specified fluxes: 411 ua_e(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * hur_e(2:nbghostcells+1,jj) 412 ! Characteristics method (only if ghostcells=1): 413 !alt ua_e(2,jj) = 0.5_wp * ( ubdy_w(jj) * hur_e(2,jj) + ua_e(3,jj) & 414 !alt & - sqrt(grav * hur_e(2,jj)) * (sshn_e(3,jj) - hbdy_w(jj)) ) 415 END DO 416 ENDIF 417 ! 418 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) 419 370 DO jj=1,jpj 420 va_e(nlci-nbghostcells:nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * hvr_e(nlci-nbghostcells:nlci-1,jj) 421 ! Specified fluxes: 422 ua_e(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * hur_e(nlci-nbghostcells-1:nlci-2,jj) 423 ! Characteristics method (only if ghostcells=1): 424 !alt ua_e(nlci-2,jj) = 0.5_wp * ( ubdy_e(jj) * hur_e(nlci-2,jj) + ua_e(nlci-3,jj) & 425 !alt & + sqrt(grav * hur_e(nlci-2,jj)) * (sshn_e(nlci-2,jj) - hbdy_e(jj)) ) 426 END DO 427 ENDIF 428 ! 429 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) 430 386 DO ji=1,jpi 431 ua_e(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * hur_e(ji,2:nbghostcells+1) 432 ! Specified fluxes: 433 va_e(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * hvr_e(ji,2:nbghostcells+1) 434 ! Characteristics method (only if ghostcells=1): 435 !alt va_e(ji,2) = 0.5_wp * ( vbdy_s(ji) * hvr_e(ji,2) + va_e(ji,3) & 436 !alt & - sqrt(grav * hvr_e(ji,2)) * (sshn_e(ji,3) - hbdy_s(ji)) ) 437 END DO 438 ENDIF 439 ! 440 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) 441 396 DO ji=1,jpi 442 ua_e(ji,nlcj-nbghostcells:nlcj-1) = ubdy_n(ji,1:nbghostcells) * hur_e(ji,nlcj-nbghostcells:nlcj-1) 443 ! Specified fluxes: 444 va_e(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * hvr_e(ji,nlcj-nbghostcells-1:nlcj-2) 445 ! Characteristics method (only if ghostcells=1): 446 !alt va_e(ji,nlcj-2) = 0.5_wp * ( vbdy_n(ji) * hvr_e(ji,nlcj-2) + va_e(ji,nlcj-3) & 447 !alt & + sqrt(grav * hvr_e(ji,nlcj-2)) * (sshn_e(ji,nlcj-2) - hbdy_n(ji)) ) 448 END DO 449 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 450 407 ! 451 408 END SUBROUTINE Agrif_dyn_ts 452 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 453 476 454 477 SUBROUTINE Agrif_dta_ts( kt ) … … 470 493 ! 471 494 ! Interpolate barotropic fluxes 472 Agrif_SpecialValue =0._wp495 Agrif_SpecialValue = 0._wp 473 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 474 501 ! 475 502 IF( ll_int_cons ) THEN ! Conservative interpolation 476 503 ! order matters here !!!!!! 477 504 CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b ) ! Time integrated 478 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 ! 479 507 bdy_tinterp = 1 480 508 CALL Agrif_Bc_variable( unb_id , calledweight=1._wp, procname=interpunb ) ! After 481 CALL Agrif_Bc_variable( vnb_id , calledweight=1._wp, procname=interpvnb ) 509 CALL Agrif_Bc_variable( vnb_id , calledweight=1._wp, procname=interpvnb ) 510 ! 482 511 bdy_tinterp = 2 483 512 CALL Agrif_Bc_variable( unb_id , calledweight=0._wp, procname=interpunb ) ! Before 484 CALL Agrif_Bc_variable( vnb_id , calledweight=0._wp, procname=interpvnb ) 513 CALL Agrif_Bc_variable( vnb_id , calledweight=0._wp, procname=interpvnb ) 485 514 ELSE ! Linear interpolation 486 bdy_tinterp = 0 487 ubdy_w(:,:) = 0._wp ; vbdy_w(:,:) = 0._wp 488 ubdy_e(:,:) = 0._wp ; vbdy_e(:,:) = 0._wp 489 ubdy_n(:,:) = 0._wp ; vbdy_n(:,:) = 0._wp 490 ubdy_s(:,:) = 0._wp ; vbdy_s(:,:) = 0._wp 515 ! 516 ubdy(:,:) = 0._wp ; vbdy(:,:) = 0._wp 491 517 CALL Agrif_Bc_variable( unb_id, procname=interpunb ) 492 518 CALL Agrif_Bc_variable( vnb_id, procname=interpvnb ) … … 503 529 INTEGER, INTENT(in) :: kt 504 530 ! 505 INTEGER :: ji, jj, indx, indy 531 INTEGER :: ji, jj 532 INTEGER :: istart, iend, jstart, jend 506 533 !!---------------------------------------------------------------------- 507 534 ! … … 516 543 ! 517 544 ! --- West --- ! 518 IF((nbondi == -1).OR.(nbondi == 2)) THEN 519 indx = 1+nbghostcells 545 istart = 2 546 iend = 1 + nbghostcells 547 DO ji = mi0(istart), mi1(iend) 520 548 DO jj = 1, jpj 521 DO ji = 2, indx 522 ssha(ji,jj) = hbdy_w(ji-1,jj) 523 ENDDO 549 ssha(ji,jj) = hbdy(ji,jj) 524 550 ENDDO 525 END IF551 ENDDO 526 552 ! 527 553 ! --- East --- ! 528 IF((nbondi == 1).OR.(nbondi == 2)) THEN 529 indx = nlci-nbghostcells 554 istart = jpiglo - nbghostcells 555 iend = jpiglo - 1 556 DO ji = mi0(istart), mi1(iend) 530 557 DO jj = 1, jpj 531 DO ji = indx, nlci-1 532 ssha(ji,jj) = hbdy_e(ji-indx+1,jj) 533 ENDDO 558 ssha(ji,jj) = hbdy(ji,jj) 534 559 ENDDO 535 END IF560 ENDDO 536 561 ! 537 562 ! --- South --- ! 538 IF((nbondj == -1).OR.(nbondj == 2)) THEN 539 indy = 1+nbghostcells 540 DO jj = 2, indy 541 DO ji = 1, jpi 542 ssha(ji,jj) = hbdy_s(ji,jj-1) 543 ENDDO 563 jstart = 2 564 jend = 1 + nbghostcells 565 DO jj = mj0(jstart), mj1(jend) 566 DO ji = 1, jpi 567 ssha(ji,jj) = hbdy(ji,jj) 544 568 ENDDO 545 END IF569 ENDDO 546 570 ! 547 571 ! --- North --- ! 548 IF((nbondj == 1).OR.(nbondj == 2)) THEN 549 indy = nlcj-nbghostcells 550 DO jj = indy, nlcj-1 551 DO ji = 1, jpi 552 ssha(ji,jj) = hbdy_n(ji,jj-indy+1) 553 ENDDO 572 jstart = jpjglo - nbghostcells 573 jend = jpjglo - 1 574 DO jj = mj0(jstart), mj1(jend) 575 DO ji = 1, jpi 576 ssha(ji,jj) = hbdy(ji,jj) 554 577 ENDDO 555 END IF578 ENDDO 556 579 ! 557 580 END SUBROUTINE Agrif_ssh … … 564 587 INTEGER, INTENT(in) :: jn 565 588 !! 566 INTEGER :: ji, jj , indx, indy567 !!----------------------------------------------------------------------568 !! 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 !!---------------------------------------------------------------------- 569 592 ! 570 593 IF( Agrif_Root() ) RETURN 571 594 ! 572 595 ! --- West --- ! 573 IF((nbondi == -1).OR.(nbondi == 2)) THEN 574 indx = 1+nbghostcells 596 istart = 2 597 iend = 1+nbghostcells 598 DO ji = mi0(istart), mi1(iend) 575 599 DO jj = 1, jpj 576 DO ji = 2, indx 577 ssha_e(ji,jj) = hbdy_w(ji-1,jj) 578 ENDDO 600 ssha_e(ji,jj) = hbdy(ji,jj) 579 601 ENDDO 580 END IF602 ENDDO 581 603 ! 582 604 ! --- East --- ! 583 IF((nbondi == 1).OR.(nbondi == 2)) THEN 584 indx = nlci-nbghostcells 605 istart = jpiglo - nbghostcells 606 iend = jpiglo - 1 607 DO ji = mi0(istart), mi1(iend) 585 608 DO jj = 1, jpj 586 DO ji = indx, nlci-1 587 ssha_e(ji,jj) = hbdy_e(ji-indx+1,jj) 588 ENDDO 609 ssha_e(ji,jj) = hbdy(ji,jj) 589 610 ENDDO 590 END IF611 ENDDO 591 612 ! 592 613 ! --- South --- ! 593 IF((nbondj == -1).OR.(nbondj == 2)) THEN 594 indy = 1+nbghostcells 595 DO jj = 2, indy 596 DO ji = 1, jpi 597 ssha_e(ji,jj) = hbdy_s(ji,jj-1) 598 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) 599 619 ENDDO 600 END IF620 ENDDO 601 621 ! 602 622 ! --- North --- ! 603 IF((nbondj == 1).OR.(nbondj == 2)) THEN 604 indy = nlcj-nbghostcells 605 DO jj = indy, nlcj-1 606 DO ji = 1, jpi 607 ssha_e(ji,jj) = hbdy_n(ji,jj-indy+1) 608 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) 609 628 ENDDO 610 END IF629 ENDDO 611 630 ! 612 631 END SUBROUTINE Agrif_ssh_ts … … 634 653 635 654 636 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 ) 637 656 !!---------------------------------------------------------------------- 638 657 !! *** ROUTINE interptsn *** … … 641 660 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2 642 661 LOGICAL , INTENT(in ) :: before 643 INTEGER , INTENT(in ) :: nb , ndir 644 ! 645 INTEGER :: ji, jj, jk, jn, iref, jref, ibdy, jbdy ! dummy loop indices 646 INTEGER :: imin, imax, jmin, jmax, N_in, N_out 647 REAL(wp) :: zrho, z1, z2, z3, z4, z5, z6, z7 648 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 649 665 ! vertical interpolation: 650 REAL(wp) , DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child651 REAL(wp), DIMENSION(k1:k2, n1:n2-1) :: tabin666 REAL(wp) :: zhtot 667 REAL(wp), DIMENSION(k1:k2,1:jpts) :: tabin 652 668 REAL(wp), DIMENSION(k1:k2) :: h_in 653 669 REAL(wp), DIMENSION(1:jpk) :: h_out 654 REAL(wp) :: h_diff670 !!---------------------------------------------------------------------- 655 671 656 672 IF( before ) THEN … … 666 682 667 683 # if defined key_vertical 684 ! Interpolate thicknesses 685 ! Warning: these are masked, hence extrapolated prior interpolation. 668 686 DO jk=k1,k2 669 687 DO jj=j1,j2 670 688 DO ji=i1,i2 671 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk)689 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 672 690 END DO 673 691 END DO 674 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) = sshn(i1:i2,j1:j2)*tmask(i1:i2,j1:j2,1) 708 ELSE 709 ptab(i1:i2,j1:j2,k2,jpts+1) = 0._wp 710 END IF 675 711 # endif 676 712 ELSE 677 713 678 western_side = (nb == 1).AND.(ndir == 1) ; eastern_side = (nb == 1).AND.(ndir == 2) 679 southern_side = (nb == 2).AND.(ndir == 1) ; northern_side = (nb == 2).AND.(ndir == 2) 680 681 # if defined key_vertical 714 # if defined key_vertical 715 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp 716 682 717 DO jj=j1,j2 683 718 DO ji=i1,i2 684 iref = ji685 jref = jj686 if(western_side) iref=MAX(2,ji)687 if(eastern_side) iref=MIN(nlci-1,ji)688 if(southern_side) jref=MAX(2,jj)689 if(northern_side) jref=MIN(nlcj-1,jj)690 N_in = 0691 DO jk=k1,k2 !k2 = jpk of parent grid692 IF (ptab(ji,jj,jk,n2) == 0) EXIT693 N_in = N_in + 1719 tsa(ji,jj,:,:) = 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) 694 729 tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1) 695 h_in(N_in) = ptab(ji,jj,jk,n2)696 730 END DO 697 731 N_out = 0 698 732 DO jk=1,jpk ! jpk of child grid 699 IF (tmask( iref,jref,jk) == 0) EXIT733 IF (tmask(ji,jj,jk) == 0._wp) EXIT 700 734 N_out = N_out + 1 701 h_out(jk) = e3t_ n(iref,jref,jk)735 h_out(jk) = e3t_a(ji,jj,jk) 702 736 ENDDO 703 IF (N_in > 0) THEN 704 DO jn=1,jpts 705 call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 706 ENDDO 737 IF (N_in*N_out > 0) THEN 738 CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tsa(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts) 707 739 ENDIF 708 740 ENDDO 709 741 ENDDO 710 742 # else 711 ptab_child(i1:i2,j1:j2,1:jpk,1:jpts) = ptab(i1:i2,j1:j2,1:jpk,1:jpts)712 # endif713 743 ! 714 744 DO jn=1, jpts 715 tsa(i1:i2,j1:j2,1:jpk,jn)=ptab_child(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 716 END DO 717 718 IF ( .NOT.lk_agrif_clp ) THEN 719 ! 720 imin = i1 ; imax = i2 721 jmin = j1 ; jmax = j2 722 ! 723 ! Remove CORNERS 724 IF((nbondj == -1).OR.(nbondj == 2)) jmin = 2 + nbghostcells 725 IF((nbondj == +1).OR.(nbondj == 2)) jmax = nlcj - nbghostcells - 1 726 IF((nbondi == -1).OR.(nbondi == 2)) imin = 2 + nbghostcells 727 IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci - nbghostcells - 1 728 ! 729 IF( eastern_side ) THEN 730 zrho = Agrif_Rhox() 731 z1 = ( zrho - 1._wp ) * 0.5_wp 732 z3 = ( zrho - 1._wp ) / ( zrho + 1._wp ) 733 z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp ) 734 z7 = - ( zrho - 1._wp ) / ( zrho + 3._wp ) 735 z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7 736 ! 737 ibdy = nlci-nbghostcells 738 DO jn = 1, jpts 739 tsa(ibdy+1,jmin:jmax,1:jpkm1,jn) = z1 * ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn) 740 DO jk = 1, jpkm1 741 DO jj = jmin,jmax 742 IF( umask(ibdy-1,jj,jk) == 0._wp ) THEN 743 tsa(ibdy,jj,jk,jn) = tsa(ibdy+1,jj,jk,jn) * tmask(ibdy,jj,jk) 744 ELSE 745 tsa(ibdy,jj,jk,jn)=(z4*tsa(ibdy+1,jj,jk,jn)+z3*tsa(ibdy-1,jj,jk,jn))*tmask(ibdy,jj,jk) 746 IF( un(ibdy-1,jj,jk) > 0._wp ) THEN 747 tsa(ibdy,jj,jk,jn)=( z6*tsa(ibdy-1,jj,jk,jn)+z5*tsa(ibdy+1,jj,jk,jn) & 748 + z7*tsa(ibdy-2,jj,jk,jn) ) * tmask(ibdy,jj,jk) 749 ENDIF 750 ENDIF 751 END DO 752 END DO 753 ! Restore ghost points: 754 tsa(ibdy+1,jmin:jmax,1:jpkm1,jn) = ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy+1,jmin:jmax,1:jpkm1) 755 END DO 756 ENDIF 757 ! 758 IF( northern_side ) THEN 759 zrho = Agrif_Rhoy() 760 z1 = ( zrho - 1._wp ) * 0.5_wp 761 z3 = ( zrho - 1._wp ) / ( zrho + 1._wp ) 762 z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp ) 763 z7 = - ( zrho - 1._wp ) / ( zrho + 3._wp ) 764 z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7 765 ! 766 jbdy = nlcj-nbghostcells 767 DO jn = 1, jpts 768 tsa(imin:imax,jbdy+1,1:jpkm1,jn) = z1 * ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn) 769 DO jk = 1, jpkm1 770 DO ji = imin,imax 771 IF( vmask(ji,jbdy-1,jk) == 0._wp ) THEN 772 tsa(ji,jbdy,jk,jn) = tsa(ji,jbdy+1,jk,jn) * tmask(ji,jbdy,jk) 773 ELSE 774 tsa(ji,jbdy,jk,jn)=(z4*tsa(ji,jbdy+1,jk,jn)+z3*tsa(ji,jbdy-1,jk,jn))*tmask(ji,jbdy,jk) 775 IF (vn(ji,jbdy-1,jk) > 0._wp ) THEN 776 tsa(ji,jbdy,jk,jn)=( z6*tsa(ji,jbdy-1,jk,jn)+z5*tsa(ji,jbdy+1,jk,jn) & 777 + z7*tsa(ji,jbdy-2,jk,jn) ) * tmask(ji,jbdy,jk) 778 ENDIF 779 ENDIF 780 END DO 781 END DO 782 ! Restore ghost points: 783 tsa(imin:imax,jbdy+1,1:jpkm1,jn) = ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) * tmask(imin:imax,jbdy+1,1:jpkm1) 784 END DO 785 ENDIF 786 ! 787 IF( western_side ) THEN 788 zrho = Agrif_Rhox() 789 z1 = ( zrho - 1._wp ) * 0.5_wp 790 z3 = ( zrho - 1._wp ) / ( zrho + 1._wp ) 791 z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp ) 792 z7 = - ( zrho - 1._wp ) / ( zrho + 3._wp ) 793 z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7 794 ! 795 ibdy = 1+nbghostcells 796 DO jn = 1, jpts 797 tsa(ibdy-1,jmin:jmax,1:jpkm1,jn) = z1 * ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn) 798 DO jk = 1, jpkm1 799 DO jj = jmin,jmax 800 IF( umask(ibdy,jj,jk) == 0._wp ) THEN 801 tsa(ibdy,jj,jk,jn) = tsa(ibdy-1,jj,jk,jn) * tmask(ibdy,jj,jk) 802 ELSE 803 tsa(ibdy,jj,jk,jn)=(z4*tsa(ibdy-1,jj,jk,jn)+z3*tsa(ibdy+1,jj,jk,jn))*tmask(ibdy,jj,jk) 804 IF( un(ibdy,jj,jk) < 0._wp ) THEN 805 tsa(ibdy,jj,jk,jn)=( z6*tsa(ibdy+1,jj,jk,jn)+z5*tsa(ibdy-1,jj,jk,jn) & 806 + z7*tsa(ibdy+2,jj,jk,jn) ) * tmask(ibdy,jj,jk) 807 ENDIF 808 ENDIF 809 END DO 810 END DO 811 ! Restore ghost points: 812 tsa(ibdy-1,jmin:jmax,1:jpkm1,jn) = ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy-1,jmin:jmax,1:jpkm1) 813 END DO 814 ENDIF 815 ! 816 IF( southern_side ) THEN 817 zrho = Agrif_Rhoy() 818 z1 = ( zrho - 1._wp ) * 0.5_wp 819 z3 = ( zrho - 1._wp ) / ( zrho + 1._wp ) 820 z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp ) 821 z7 = - ( zrho - 1._wp ) / ( zrho + 3._wp ) 822 z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7 823 ! 824 jbdy=1+nbghostcells 825 DO jn = 1, jpts 826 tsa(imin:imax,jbdy-1,1:jpkm1,jn) = z1 * ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn) 827 DO jk = 1, jpkm1 828 DO ji = imin,imax 829 IF( vmask(ji,jbdy,jk) == 0._wp ) THEN 830 tsa(ji,jbdy,jk,jn)=tsa(ji,jbdy-1,jk,jn) * tmask(ji,jbdy,jk) 831 ELSE 832 tsa(ji,jbdy,jk,jn)=(z4*tsa(ji,jbdy-1,jk,jn)+z3*tsa(ji,jbdy+1,jk,jn))*tmask(ji,jbdy,jk) 833 IF( vn(ji,jbdy,jk) < 0._wp ) THEN 834 tsa(ji,jbdy,jk,jn)=( z6*tsa(ji,jbdy+1,jk,jn)+z5*tsa(ji,jbdy-1,jk,jn) & 835 + z7*tsa(ji,jbdy+2,jk,jn) ) * tmask(ji,jbdy,jk) 836 ENDIF 837 ENDIF 838 END DO 839 END DO 840 ! Restore ghost points: 841 tsa(imin:imax,jbdy-1,1:jpkm1,jn) = ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) * tmask(imin:imax,jbdy-1,1:jpkm1) 842 END DO 843 ENDIF 844 ! 845 ENDIF 745 tsa(i1:i2,j1:j2,1:jpk,jn)=ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 746 END DO 747 # endif 748 846 749 ENDIF 847 750 ! 848 751 END SUBROUTINE interptsn 849 752 850 SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before , nb, ndir)753 SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before ) 851 754 !!---------------------------------------------------------------------- 852 755 !! *** ROUTINE interpsshn *** … … 855 758 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 856 759 LOGICAL , INTENT(in ) :: before 857 INTEGER , INTENT(in ) :: nb , ndir 858 ! 859 LOGICAL :: western_side, eastern_side,northern_side,southern_side 760 ! 860 761 !!---------------------------------------------------------------------- 861 762 ! … … 863 764 ptab(i1:i2,j1:j2) = sshn(i1:i2,j1:j2) 864 765 ELSE 865 western_side = (nb == 1).AND.(ndir == 1) 866 eastern_side = (nb == 1).AND.(ndir == 2) 867 southern_side = (nb == 2).AND.(ndir == 1) 868 northern_side = (nb == 2).AND.(ndir == 2) 869 !! clem ghost 870 IF(western_side) hbdy_w(1:nbghostcells,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1) 871 IF(eastern_side) hbdy_e(1:nbghostcells,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1) 872 IF(southern_side) hbdy_s(i1:i2,1:nbghostcells) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1) 873 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) 874 767 ENDIF 875 768 ! 876 769 END SUBROUTINE interpsshn 877 770 878 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 ) 879 772 !!---------------------------------------------------------------------- 880 773 !! *** ROUTINE interpun *** … … 884 777 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 885 778 LOGICAL, INTENT(in) :: before 886 INTEGER, INTENT(in) :: nb , ndir887 779 !! 888 780 INTEGER :: ji,jj,jk 889 REAL(wp) :: zrhoy 781 REAL(wp) :: zrhoy, zhtot 890 782 ! vertical interpolation: 891 783 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 892 784 REAL(wp), DIMENSION(1:jpk) :: h_out 893 INTEGER :: N_in, N_out , iref785 INTEGER :: N_in, N_out 894 786 REAL(wp) :: h_diff 895 LOGICAL :: western_side, eastern_side896 787 !!--------------------------------------------- 897 788 ! … … 902 793 ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk)) 903 794 # if defined key_vertical 904 ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)) 795 ! Interpolate thicknesses (masked for subsequent extrapolation) 796 ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk) 905 797 # endif 906 798 END DO 907 799 END DO 908 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_n(i1:i2,j1:j2,jk) * 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 ! 909 823 ELSE 910 824 zrhoy = Agrif_rhoy() 911 825 # if defined key_vertical 912 826 ! VERTICAL REFINEMENT BEGIN 913 western_side = (nb == 1).AND.(ndir == 1) 914 eastern_side = (nb == 1).AND.(ndir == 2)827 828 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 915 829 916 830 DO ji=i1,i2 917 iref = ji918 IF (western_side) iref = MAX(2,ji)919 IF (eastern_side) iref = MIN(nlci-2,ji)920 831 DO jj=j1,j2 921 N_in = 0 922 DO jk=k1,k2 923 IF (ptab(ji,jj,jk,2) == 0) EXIT 924 N_in = N_in + 1 925 tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2) 926 h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 832 ua(ji,jj,:) = 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)) 927 843 ENDDO 928 929 IF (N_in == 0) THEN 930 ua(ji,jj,:) = 0._wp 931 CYCLE 932 ENDIF 933 844 934 845 N_out = 0 935 846 DO jk=1,jpk 936 if (umask( iref,jj,jk) == 0) EXIT847 if (umask(ji,jj,jk) == 0) EXIT 937 848 N_out = N_out + 1 938 h_out(N_out) = e3u_a( iref,jj,jk)849 h_out(N_out) = e3u_a(ji,jj,jk) 939 850 ENDDO 940 941 IF (N_out == 0) THEN 942 ua(ji,jj,:) = 0._wp 943 CYCLE 851 IF (N_in*N_out > 0) THEN 852 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),ua(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 944 853 ENDIF 945 946 IF (N_in * N_out > 0) THEN947 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))948 ! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly949 if (h_diff < -1.e4) then950 print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))951 ! stop952 endif953 ENDIF954 call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ua(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)955 854 ENDDO 956 855 ENDDO … … 968 867 END SUBROUTINE interpun 969 868 970 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 ) 971 870 !!---------------------------------------------------------------------- 972 871 !! *** ROUTINE interpvn *** … … 976 875 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 977 876 LOGICAL, INTENT(in) :: before 978 INTEGER, INTENT(in) :: nb , ndir979 877 ! 980 878 INTEGER :: ji,jj,jk … … 983 881 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 984 882 REAL(wp), DIMENSION(1:jpk) :: h_out 985 INTEGER :: N_in, N_out, jref 986 REAL(wp) :: h_diff 987 LOGICAL :: northern_side,southern_side 883 INTEGER :: N_in, N_out 884 REAL(wp) :: h_diff, zhtot 988 885 !!--------------------------------------------- 989 886 ! … … 994 891 ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)*vmask(ji,jj,jk)) 995 892 # if defined key_vertical 893 ! Interpolate thicknesses (masked for subsequent extrapolation) 996 894 ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk) 997 895 # endif … … 999 897 END DO 1000 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_n(i1:i2,j1:j2,jk) * 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 1001 920 ELSE 1002 921 zrhox = Agrif_rhox() 1003 922 # if defined key_vertical 1004 923 1005 southern_side = (nb == 2).AND.(ndir == 1) 1006 northern_side = (nb == 2).AND.(ndir == 2) 924 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 1007 925 1008 926 DO jj=j1,j2 1009 jref = jj1010 IF (southern_side) jref = MAX(2,jj)1011 IF (northern_side) jref = MIN(nlcj-2,jj)1012 927 DO ji=i1,i2 1013 N_in = 0 1014 DO jk=k1,k2 1015 if (ptab(ji,jj,jk,2) == 0) EXIT 1016 N_in = N_in + 1 1017 tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2) 1018 h_in(N_in) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 1019 END DO 1020 IF (N_in == 0) THEN 1021 va(ji,jj,:) = 0._wp 1022 CYCLE 1023 ENDIF 928 va(ji,jj,:) = 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 1024 940 1025 941 N_out = 0 1026 942 DO jk=1,jpk 1027 if (vmask(ji,j ref,jk) == 0) EXIT943 if (vmask(ji,jj,jk) == 0) EXIT 1028 944 N_out = N_out + 1 1029 h_out(N_out) = e3v_a(ji,jref,jk) 1030 END DO 1031 IF (N_out == 0) THEN 1032 va(ji,jj,:) = 0._wp 1033 CYCLE 945 h_out(N_out) = e3v_a(ji,jj,jk) 946 END DO 947 IF (N_in*N_out > 0) THEN 948 call reconstructandremap(tabin(1:N_in),h_in(1:N_in),va(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 1034 949 ENDIF 1035 call reconstructandremap(tabin(1:N_in),h_in(1:N_in),va(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)1036 950 END DO 1037 951 END DO … … 1045 959 END SUBROUTINE interpvn 1046 960 1047 SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before , nb, ndir)961 SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before) 1048 962 !!---------------------------------------------------------------------- 1049 963 !! *** ROUTINE interpunb *** … … 1052 966 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 1053 967 LOGICAL , INTENT(in ) :: before 1054 INTEGER , INTENT(in ) :: nb , ndir1055 968 ! 1056 969 INTEGER :: ji, jj 1057 970 REAL(wp) :: zrhoy, zrhot, zt0, zt1, ztcoeff 1058 LOGICAL :: western_side, eastern_side,northern_side,southern_side1059 971 !!---------------------------------------------------------------------- 1060 972 ! … … 1062 974 ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu_n(i1:i2,j1:j2) * un_b(i1:i2,j1:j2) 1063 975 ELSE 1064 western_side = (nb == 1).AND.(ndir == 1)1065 eastern_side = (nb == 1).AND.(ndir == 2)1066 southern_side = (nb == 2).AND.(ndir == 1)1067 northern_side = (nb == 2).AND.(ndir == 2)1068 976 zrhoy = Agrif_Rhoy() 1069 977 zrhot = Agrif_rhot() … … 1071 979 zt0 = REAL(Agrif_NbStepint() , wp) / zrhot 1072 980 zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot 1073 ! Polynomial interpolation coefficients: 1074 IF( bdy_tinterp == 1 ) THEN 1075 ztcoeff = zrhot * ( zt1**2._wp * ( zt1 - 1._wp) & 1076 & - zt0**2._wp * ( zt0 - 1._wp) ) 1077 ELSEIF( bdy_tinterp == 2 ) THEN 1078 ztcoeff = zrhot * ( zt1 * ( zt1 - 1._wp)**2._wp & 1079 & - zt0 * ( zt0 - 1._wp)**2._wp ) 1080 ELSE 1081 ztcoeff = 1 1082 ENDIF 1083 ! 1084 IF(western_side) ubdy_w(1:nbghostcells,j1:j2) = ubdy_w(1:nbghostcells,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2) 1085 IF(eastern_side) ubdy_e(1:nbghostcells,j1:j2) = ubdy_e(1:nbghostcells,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2) 1086 IF(southern_side) ubdy_s(i1:i2,1:nbghostcells) = ubdy_s(i1:i2,1:nbghostcells) + ztcoeff * ptab(i1:i2,j1:j2) 1087 IF(northern_side) ubdy_n(i1:i2,1:nbghostcells) = ubdy_n(i1:i2,1:nbghostcells) + ztcoeff * ptab(i1:i2,j1:j2) 1088 ! 1089 IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN 1090 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) 1091 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) 1092 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) 1093 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) 1094 ENDIF 1095 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 1096 1008 ! 1097 1009 END SUBROUTINE interpunb 1098 1010 1099 1011 1100 SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before , nb, ndir)1012 SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before ) 1101 1013 !!---------------------------------------------------------------------- 1102 1014 !! *** ROUTINE interpvnb *** … … 1105 1017 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 1106 1018 LOGICAL , INTENT(in ) :: before 1107 INTEGER , INTENT(in ) :: nb , ndir 1108 ! 1109 INTEGER :: ji,jj 1019 ! 1020 INTEGER :: ji, jj 1110 1021 REAL(wp) :: zrhox, zrhot, zt0, zt1, ztcoeff 1111 LOGICAL :: western_side, eastern_side,northern_side,southern_side1112 1022 !!---------------------------------------------------------------------- 1113 1023 ! … … 1115 1025 ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv_n(i1:i2,j1:j2) * vn_b(i1:i2,j1:j2) 1116 1026 ELSE 1117 western_side = (nb == 1).AND.(ndir == 1)1118 eastern_side = (nb == 1).AND.(ndir == 2)1119 southern_side = (nb == 2).AND.(ndir == 1)1120 northern_side = (nb == 2).AND.(ndir == 2)1121 1027 zrhox = Agrif_Rhox() 1122 1028 zrhot = Agrif_rhot() 1123 1029 ! Time indexes bounds for integration 1124 1030 zt0 = REAL(Agrif_NbStepint() , wp) / zrhot 1125 zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot 1126 IF( bdy_tinterp == 1 ) THEN 1127 ztcoeff = zrhot * ( zt1**2._wp * ( zt1 - 1._wp) & 1128 & - zt0**2._wp * ( zt0 - 1._wp) ) 1129 ELSEIF( bdy_tinterp == 2 ) THEN 1130 ztcoeff = zrhot * ( zt1 * ( zt1 - 1._wp)**2._wp & 1131 & - zt0 * ( zt0 - 1._wp)**2._wp ) 1132 ELSE 1133 ztcoeff = 1 1134 ENDIF 1135 !! clem ghost 1136 IF(western_side) vbdy_w(1:nbghostcells,j1:j2) = vbdy_w(1:nbghostcells,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2) 1137 IF(eastern_side) vbdy_e(1:nbghostcells,j1:j2) = vbdy_e(1:nbghostcells,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2) 1138 IF(southern_side) vbdy_s(i1:i2,1:nbghostcells) = vbdy_s(i1:i2,1:nbghostcells) + ztcoeff * ptab(i1:i2,j1:j2) 1139 IF(northern_side) vbdy_n(i1:i2,1:nbghostcells) = vbdy_n(i1:i2,1:nbghostcells) + ztcoeff * ptab(i1:i2,j1:j2) 1140 ! 1141 IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN 1142 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) 1143 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) 1144 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) 1145 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) 1146 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 1147 1058 ENDIF 1148 1059 ! … … 1150 1061 1151 1062 1152 SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before , nb, ndir)1063 SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before ) 1153 1064 !!---------------------------------------------------------------------- 1154 1065 !! *** ROUTINE interpub2b *** … … 1157 1068 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 1158 1069 LOGICAL , INTENT(in ) :: before 1159 INTEGER , INTENT(in ) :: nb , ndir1160 1070 ! 1161 1071 INTEGER :: ji,jj 1162 REAL(wp) :: zrhot, zt0, zt1,zat 1163 LOGICAL :: western_side, eastern_side,northern_side,southern_side 1072 REAL(wp) :: zrhot, zt0, zt1, zat 1164 1073 !!---------------------------------------------------------------------- 1165 1074 IF( before ) THEN … … 1170 1079 ENDIF 1171 1080 ELSE 1172 western_side = (nb == 1).AND.(ndir == 1)1173 eastern_side = (nb == 1).AND.(ndir == 2)1174 southern_side = (nb == 2).AND.(ndir == 1)1175 northern_side = (nb == 2).AND.(ndir == 2)1176 zrhot = Agrif_rhot()1177 ! Time indexes bounds for integration1178 zt0 = REAL(Agrif_NbStepint() , wp) / zrhot1179 zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot1180 ! Polynomial interpolation coefficients:1181 zat = zrhot * ( zt1**2._wp * (-2._wp*zt1 + 3._wp) &1182 & - zt0**2._wp * (-2._wp*zt0 + 3._wp) )1183 !! clem ghost1184 IF(western_side ) ubdy_w(1:nbghostcells,j1:j2) = zat * ptab(i1:i2,j1:j2)1185 IF(eastern_side ) ubdy_e(1:nbghostcells,j1:j2) = zat * ptab(i1:i2,j1:j2)1186 IF(southern_side) ubdy_s(i1:i2,1:nbghostcells) = zat * ptab(i1:i2,j1:j2)1187 IF(northern_side) ubdy_n(i1:i2,1:nbghostcells) = zat * ptab(i1:i2,j1:j2)1188 ENDIF1189 !1190 END SUBROUTINE interpub2b1191 1192 1193 SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before, nb, ndir )1194 !!----------------------------------------------------------------------1195 !! *** ROUTINE interpvb2b ***1196 !!----------------------------------------------------------------------1197 INTEGER , INTENT(in ) :: i1, i2, j1, j21198 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab1199 LOGICAL , INTENT(in ) :: before1200 INTEGER , INTENT(in ) :: nb , ndir1201 !1202 INTEGER :: ji,jj1203 REAL(wp) :: zrhot, zt0, zt1,zat1204 LOGICAL :: western_side, eastern_side,northern_side,southern_side1205 !!----------------------------------------------------------------------1206 !1207 IF( before ) THEN1208 IF ( ln_bt_fw ) THEN1209 ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)1210 ELSE1211 ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2)1212 ENDIF1213 ELSE1214 western_side = (nb == 1).AND.(ndir == 1)1215 eastern_side = (nb == 1).AND.(ndir == 2)1216 southern_side = (nb == 2).AND.(ndir == 1)1217 northern_side = (nb == 2).AND.(ndir == 2)1218 1081 zrhot = Agrif_rhot() 1219 1082 ! Time indexes bounds for integration … … 1224 1087 & - zt0**2._wp * (-2._wp*zt0 + 3._wp) ) 1225 1088 ! 1226 IF(western_side ) vbdy_w(1:nbghostcells,j1:j2) = zat * ptab(i1:i2,j1:j2) 1227 IF(eastern_side ) vbdy_e(1:nbghostcells,j1:j2) = zat * ptab(i1:i2,j1:j2) 1228 IF(southern_side) vbdy_s(i1:i2,1:nbghostcells) = zat * ptab(i1:i2,j1:j2) 1229 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 1230 1129 ENDIF 1231 1130 ! … … 1233 1132 1234 1133 1235 SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before , nb, ndir)1134 SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before ) 1236 1135 !!---------------------------------------------------------------------- 1237 1136 !! *** ROUTINE interpe3t *** … … 1240 1139 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab 1241 1140 LOGICAL , INTENT(in ) :: before 1242 INTEGER , INTENT(in ) :: nb , ndir1243 1141 ! 1244 1142 INTEGER :: ji, jj, jk 1245 LOGICAL :: western_side, eastern_side, northern_side, southern_side1246 1143 !!---------------------------------------------------------------------- 1247 1144 ! … … 1249 1146 ptab(i1:i2,j1:j2,k1:k2) = tmask(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2) 1250 1147 ELSE 1251 western_side = (nb == 1).AND.(ndir == 1)1252 eastern_side = (nb == 1).AND.(ndir == 2)1253 southern_side = (nb == 2).AND.(ndir == 1)1254 northern_side = (nb == 2).AND.(ndir == 2)1255 1148 ! 1256 1149 DO jk = k1, k2 1257 1150 DO jj = j1, j2 1258 1151 DO ji = i1, i2 1259 !1260 1152 IF( ABS( ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk) ) > 1.D-2) THEN 1261 IF (western_side.AND.(ptab(i1+nbghostcells-1,jj,jk)>0._wp)) THEN 1262 WRITE(numout,*) 'ERROR bathymetry merge at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk 1263 WRITE(numout,*) ptab(ji,jj,jk), e3t_0(ji,jj,jk) 1264 kindic_agr = kindic_agr + 1 1265 ELSEIF (eastern_side.AND.(ptab(i2-nbghostcells+1,jj,jk)>0._wp)) THEN 1266 WRITE(numout,*) 'ERROR bathymetry merge at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk 1267 WRITE(numout,*) ptab(ji,jj,jk), e3t_0(ji,jj,jk) 1268 kindic_agr = kindic_agr + 1 1269 ELSEIF (southern_side.AND.(ptab(ji,j1+nbghostcells-1,jk)>0._wp)) THEN 1270 WRITE(numout,*) 'ERROR bathymetry merge at the southern border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk 1271 WRITE(numout,*) ptab(ji,jj,jk), e3t_0(ji,jj,jk) 1272 kindic_agr = kindic_agr + 1 1273 ELSEIF (northern_side.AND.(ptab(ji,j2-nbghostcells+1,jk)>0._wp)) THEN 1274 WRITE(numout,*) 'ERROR bathymetry merge at the northen border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk 1275 WRITE(numout,*) ptab(ji,jj,jk), e3t_0(ji,jj,jk) 1276 kindic_agr = kindic_agr + 1 1277 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 1278 1157 ENDIF 1279 1158 END DO … … 1284 1163 ! 1285 1164 END SUBROUTINE interpe3t 1286 1287 1288 SUBROUTINE interpumsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )1289 !!----------------------------------------------------------------------1290 !! *** ROUTINE interpumsk ***1291 !!----------------------------------------------------------------------1292 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k21293 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab1294 LOGICAL , INTENT(in ) :: before1295 INTEGER , INTENT(in ) :: nb , ndir1296 !1297 INTEGER :: ji, jj, jk1298 LOGICAL :: western_side, eastern_side1299 !!----------------------------------------------------------------------1300 !1301 IF( before ) THEN1302 ptab(i1:i2,j1:j2,k1:k2) = umask(i1:i2,j1:j2,k1:k2)1303 ELSE1304 western_side = (nb == 1).AND.(ndir == 1)1305 eastern_side = (nb == 1).AND.(ndir == 2)1306 DO jk = k1, k21307 DO jj = j1, j21308 DO ji = i1, i21309 ! Velocity mask at boundary edge points:1310 IF (ABS(ptab(ji,jj,jk) - umask(ji,jj,jk)) > 1.D-2) THEN1311 IF (western_side) THEN1312 WRITE(numout,*) 'ERROR with umask at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk1313 WRITE(numout,*) ' masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)1314 kindic_agr = kindic_agr + 11315 ELSEIF (eastern_side) THEN1316 WRITE(numout,*) 'ERROR with umask at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk1317 WRITE(numout,*) ' masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)1318 kindic_agr = kindic_agr + 11319 ENDIF1320 ENDIF1321 END DO1322 END DO1323 END DO1324 !1325 ENDIF1326 !1327 END SUBROUTINE interpumsk1328 1329 1330 SUBROUTINE interpvmsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )1331 !!----------------------------------------------------------------------1332 !! *** ROUTINE interpvmsk ***1333 !!----------------------------------------------------------------------1334 INTEGER , INTENT(in ) :: i1,i2,j1,j2,k1,k21335 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab1336 LOGICAL , INTENT(in ) :: before1337 INTEGER , INTENT(in ) :: nb , ndir1338 !1339 INTEGER :: ji, jj, jk1340 LOGICAL :: northern_side, southern_side1341 !!----------------------------------------------------------------------1342 !1343 IF( before ) THEN1344 ptab(i1:i2,j1:j2,k1:k2) = vmask(i1:i2,j1:j2,k1:k2)1345 ELSE1346 southern_side = (nb == 2).AND.(ndir == 1)1347 northern_side = (nb == 2).AND.(ndir == 2)1348 DO jk = k1, k21349 DO jj = j1, j21350 DO ji = i1, i21351 ! Velocity mask at boundary edge points:1352 IF (ABS(ptab(ji,jj,jk) - vmask(ji,jj,jk)) > 1.D-2) THEN1353 IF (southern_side) THEN1354 WRITE(numout,*) 'ERROR with vmask at the southern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk1355 WRITE(numout,*) ' masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)1356 kindic_agr = kindic_agr + 11357 ELSEIF (northern_side) THEN1358 WRITE(numout,*) 'ERROR with vmask at the northern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk1359 WRITE(numout,*) ' masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)1360 kindic_agr = kindic_agr + 11361 ENDIF1362 ENDIF1363 END DO1364 END DO1365 END DO1366 !1367 ENDIF1368 !1369 END SUBROUTINE interpvmsk1370 1165 1371 1166 … … 1377 1172 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 1378 1173 LOGICAL , INTENT(in ) :: before 1379 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 1380 REAL(wp), DIMENSION(1:jpk) :: h_out 1381 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 1382 1179 !!---------------------------------------------------------------------- 1383 1180 ! … … 1390 1187 END DO 1391 1188 END DO 1392 #ifdef key_vertical 1189 1190 # if defined key_vertical 1191 ! Interpolate thicknesses 1192 ! Warning: these are masked, hence extrapolated prior interpolation. 1393 1193 DO jk=k1,k2 1394 1194 DO jj=j1,j2 1395 1195 DO ji=i1,i2 1396 ptab(ji,jj,jk,2) = wmask(ji,jj,jk) * e3w_n(ji,jj,jk)1196 ptab(ji,jj,jk,2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 1397 1197 END DO 1398 1198 END DO 1399 1199 END DO 1400 #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) = sshn(i1:i2,j1:j2)*tmask(i1:i2,j1:j2,1) 1215 ELSE 1216 ptab(i1:i2,j1:j2,k2,2) = 0._wp 1217 END IF 1218 # endif 1401 1219 ELSE 1402 1220 #ifdef key_vertical 1403 avm_k(i1:i2,j1:j2,1:jpk) = 0.1404 DO jj=j1,j21405 DO ji=i1,i21406 N_in = 01407 DO jk=k1,k2 !k2 = jpk of parent grid1408 IF (ptab(ji,jj,jk,2) == 0) EXIT1409 N_in = N_in + 11410 tabin(jk) = ptab(ji,jj,jk,1)1411 h_in(N_in) = ptab(ji,jj,jk,2)1412 END DO1413 N_out = 01414 DO jk=1,jpk ! jpk of child grid1415 IF (wmask(ji,jj,jk) == 0) EXIT1416 N_out = N_out + 11417 h_out(jk) = e3t_n(ji,jj,jk)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_n(ji,jj,jk) 1418 1236 ENDDO 1419 IF (N_in > 0) THEN1420 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) 1421 1239 ENDIF 1422 1240 ENDDO … … 1428 1246 ! 1429 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 1430 1286 1431 1287 #else -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/NST/agrif_oce_sponge.F90
r10425 r12191 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 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(:,:) ) 324 #endif 188 325 ! 189 326 #endif … … 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_n(ji,jj,jk) 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_b(ji,jj,jk) 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) = sshb(i1:i2,j1:j2)*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_ n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above416 h_out(jk) = e3t_b(ji,jj,jk) !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) = tsb(ji,jj,jk,1:jpts) - tabres_child(ji,jj,jk,1:jpts)436 tsbdiff(ji,jj,jk,1:jpts) = (tsb(ji,jj,jk,1:jpts) - tabres_child(ji,jj,jk,1:jpts)) * tmask(ji,jj,jk) 271 437 # else 272 tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres(ji,jj,jk,1:jpts)438 tsbdiff(ji,jj,jk,1:jpts) = (tsb(ji,jj,jk,1:jpts) - 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_n(ji,jj,jk)454 zabe1 = rn_sponge_tra * fspu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) 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_n(ji,jj,jk)461 zabe2 = rn_sponge_tra * fspv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) 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_n(ji,jj,jk) 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 tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + 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) = ub(ji,jj,jk) 356 528 # if defined key_vertical 357 tabres(ji,jj,jk,m2) = e3u_ n(ji,jj,jk)*umask(ji,jj,jk)529 tabres(ji,jj,jk,m2) = e3u_b(ji,jj,jk)*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_b(i1:i2,j1:j2,jk) * 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_n(ji,jj,jk) 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_b(ji,jj,jk) 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 … … 407 599 ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - 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_ n(ji,jj,jk) * fsahm_spt(ji,jj)419 hdivdiff(ji,jj,jk) = ( e2u(ji ,jj)*e3u_ n(ji ,jj,jk) * ubdiff(ji ,jj,jk) &420 & -e2u(ji-1,jj)*e3u_ n(ji-1,jj,jk) * ubdiff(ji-1,jj,jk) ) * zbtr614 zbtr = r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk) * rn_sponge_dyn * fspt(ji,jj) 615 hdivdiff(ji,jj,jk) = ( e2u(ji ,jj)*e3u_b(ji ,jj,jk) * ubdiff(ji ,jj,jk) & 616 & -e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * 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_n(ji,jj,jk) * fsahm_spf(ji,jj)622 zbtr = r1_e1e2f(ji,jj) * e3f_n(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_n(ji,jj,jk) ) & 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 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 446 642 ua(ji,jj,jk) = ua(ji,jj,jk) + 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) = vb(ji,jj,jk) 509 705 # if defined key_vertical 510 tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v_ n(ji,jj,jk)706 tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v_b(ji,jj,jk) 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_b(i1:i2,j1:j2,jk) * 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_b(ji,jj,jk) 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_n(ji,jj,jk) 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 … … 553 775 vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - 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_ n(ji,jj,jk) * fsahm_spt(ji,jj)565 hdivdiff(ji,jj,jk) = ( e1v(ji,jj ) * e3v_ n(ji,jj ,jk) * vbdiff(ji,jj ,jk) &566 & -e1v(ji,jj-1) * e3v_ n(ji,jj-1,jk) * vbdiff(ji,jj-1,jk) ) * zbtr790 zbtr = r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk) * rn_sponge_dyn * fspt(ji,jj) 791 hdivdiff(ji,jj,jk) = ( e1v(ji,jj ) * e3v_b(ji,jj ,jk) * vbdiff(ji,jj ,jk) & 792 & -e1v(ji,jj-1) * e3v_b(ji,jj-1,jk) * 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_n(ji,jj,jk) * fsahm_spf(ji,jj)797 zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * rn_sponge_dyn * fspf(ji,jj) 572 798 rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) & 573 799 & -e2v(ji ,jj) * vbdiff(ji ,jj,jk) ) * fmask(ji,jj,jk) * zbtr … … 602 828 va(ji,jj,jk) = va(ji,jj,jk) & 603 829 & + ( rotdiff (ji,jj ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) ) & 604 & + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji ,jj,jk) ) * r1_e2v(ji,jj) 830 & + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji ,jj,jk) ) * r1_e2v(ji,jj) & 831 & - ztrelax * fspv(ji,jj) * vbdiff(ji,jj,jk) 605 832 END DO 606 833 ENDIF -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/NST/agrif_oce_update.F90
r10068 r12191 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) = (tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) ) & 318 * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp 309 !jc_alt 310 ! tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) ) & 311 ! & * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1._wp) * 999._wp 312 tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) 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_n(ji,jj,jk) & 327 + (tmask(ji,jj,jk)-1)*999._wp 320 !jc_alt 321 ! tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) & 322 ! & + (tmask(ji,jj,jk) - 1._wp) * 999._wp 323 tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 328 324 END DO 329 325 END DO 330 326 END DO 331 327 ELSE 332 tabres_child(:,:,:,:) = 0. 328 tabres_child(:,:,:,:) = 0._wp 333 329 AGRIF_SpecialValue = 0._wp 334 330 DO jj=j1,j2 … … 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_n(ji,jj,jk) 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),h_out(1:N_out),N_in,N_out) 359 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) 360 349 ENDIF 361 350 ENDDO … … 364 353 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 365 354 ! Add asselin part 366 DO jn = n1,n2-1 367 DO jk=1,jpk 368 DO jj=j1,j2 369 DO ji=i1,i2 370 IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 371 tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 372 & + atfp * ( tabres_child(ji,jj,jk,jn) & 373 & - tsn(ji,jj,jk,jn) ) * 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 = tsb(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! fse3t_b prior update should be used 361 ztnu = tabres_child(ji,jj,jk,jn) * e3t_n(ji,jj,jk) 362 ztno = tsn(ji,jj,jk,jn) * e3t_a(ji,jj,jk) 363 tsb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) ) & 364 & * tmask(ji,jj,jk) / e3t_b(ji,jj,jk) 374 365 ENDIF 375 END DO376 END DO377 END DO378 END DO379 ENDIF 380 DO jn = n1,n2-1381 DO jk =1,jpk382 DO jj =j1,j2383 DO ji =i1,i2384 IF( tabres_child(ji,jj,jk,jn) .NE. 0.) THEN385 tsn(ji,jj,jk,jn) = 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 tsn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) 386 377 END IF 387 378 END DO … … 389 380 END DO 390 381 END DO 382 ! 383 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 384 tsb(i1:i2,j1:j2,1:jpkm1,1:jpts) = tsn(i1:i2,j1:j2,1:jpkm1,1:jpts) 385 ENDIF 391 386 ENDIF 392 387 ! … … 478 473 ! 479 474 INTEGER :: ji, jj, jk 480 REAL(wp):: zrhoy 475 REAL(wp):: zrhoy, zub, zunu, zuno 481 476 ! VERTICAL REFINEMENT BEGIN 482 477 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child … … 491 486 IF( before ) THEN 492 487 zrhoy = Agrif_Rhoy() 493 AGRIF_SpecialValue = -999._wp 488 !jc_alt 489 ! AGRIF_SpecialValue = -999._wp 494 490 DO jk=k1,k2 495 491 DO jj=j1,j2 496 492 DO ji=i1,i2 497 tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk) & 498 + (umask(ji,jj,jk)-1)*999._wp 499 tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk) & 500 + (umask(ji,jj,jk)-1)*999._wp 493 !jc_alt 494 ! tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk) & 495 ! & + (umask(ji,jj,jk)-1._wp)*999._wp 496 tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk) 497 !jc_alt 498 ! tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk) & 499 ! & + (umask(ji,jj,jk)-1._wp)*999._wp 500 tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk) 501 501 END DO 502 502 END DO … … 511 511 tabin(:) = 0._wp 512 512 DO jk=k1,k2 !k2=jpk of child grid 513 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 514 516 N_in = N_in + 1 515 517 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) … … 524 526 IF (N_in * N_out > 0) THEN 525 527 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 528 excess = 0._wp 526 529 IF (h_diff < -1.e-4) THEN 527 530 !Even if bathy at T points match it's possible for the U points to be deeper in the child grid. 528 531 !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 529 excess = 0._wp530 532 DO jk=N_in,1,-1 531 533 thick = MIN(-1*h_diff, h_in(jk)) … … 540 542 ENDDO 541 543 ENDIF 542 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 )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) 543 545 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out)) 544 546 ENDIF 545 547 ENDDO 546 548 ENDDO 547 549 ! 548 550 DO jk=1,jpk 549 551 DO jj=j1,j2 550 552 DO ji=i1,i2 551 553 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 552 ub(ji,jj,jk) = ub(ji,jj,jk) & 553 & + atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk) 554 zub = ub(ji,jj,jk) * e3u_b(ji,jj,jk) ! fse3t_b prior update should be used 555 zuno = un(ji,jj,jk) * e3u_a(ji,jj,jk) 556 zunu = tabres_child(ji,jj,jk) * e3u_n(ji,jj,jk) 557 ub(ji,jj,jk) = ( zub + atfp * ( zunu - zuno) ) & 558 & * umask(ji,jj,jk) / e3u_b(ji,jj,jk) 554 559 ENDIF 555 560 ! … … 558 563 END DO 559 564 END DO 565 ! 566 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 567 ub(i1:i2,j1:j2,1:jpkm1) = un(i1:i2,j1:j2,1:jpkm1) 568 ENDIF 569 ! 560 570 ENDIF 561 571 ! … … 665 675 ! 666 676 INTEGER :: ji, jj, jk 667 REAL(wp) :: zrhox 677 REAL(wp) :: zrhox, zvb, zvnu, zvno 668 678 ! VERTICAL REFINEMENT BEGIN 669 679 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child … … 678 688 IF( before ) THEN 679 689 zrhox = Agrif_Rhox() 680 AGRIF_SpecialValue = -999._wp 690 !jc_alt 691 ! AGRIF_SpecialValue = -999._wp 681 692 DO jk=k1,k2 682 693 DO jj=j1,j2 683 694 DO ji=i1,i2 684 tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) & 685 + (vmask(ji,jj,jk)-1)*999._wp 686 tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) & 687 + (vmask(ji,jj,jk)-1)*999._wp 695 !jc_alt 696 ! tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) & 697 ! & + (vmask(ji,jj,jk)-1._wp) * 999._wp 698 tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) 699 !jc_alt 700 ! tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) & 701 ! & + (vmask(ji,jj,jk)-1._wp) * 999._wp 702 tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) 688 703 END DO 689 704 END DO … … 696 711 N_in = 0 697 712 DO jk=k1,k2 698 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 699 716 N_in = N_in + 1 700 717 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) … … 709 726 IF (N_in * N_out > 0) THEN 710 727 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 728 excess = 0._wp 711 729 IF (h_diff < -1.e-4) then 712 !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. 713 731 !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 714 excess = 0._wp715 732 DO jk=N_in,1,-1 716 733 thick = MIN(-1*h_diff, h_in(jk)) … … 725 742 ENDDO 726 743 ENDIF 727 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 )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) 728 745 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out)) 729 746 ENDIF 730 747 ENDDO 731 748 ENDDO 732 733 DO jk=1,jpk 749 ! 750 DO jk=1,jpkm1 734 751 DO jj=j1,j2 735 752 DO ji=i1,i2 736 ! 737 IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN ! Add asselin part 738 vb(ji,jj,jk) = vb(ji,jj,jk) & 739 & + atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 753 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 754 zvb = vb(ji,jj,jk) * e3v_b(ji,jj,jk) ! fse3t_b prior update should be used 755 zvno = vn(ji,jj,jk) * e3v_a(ji,jj,jk) 756 zvnu = tabres_child(ji,jj,jk) * e3v_n(ji,jj,jk) 757 vb(ji,jj,jk) = ( zvb + atfp * ( zvnu - zvno) ) & 758 & * vmask(ji,jj,jk) / e3v_b(ji,jj,jk) 740 759 ENDIF 741 760 ! … … 744 763 END DO 745 764 END DO 765 ! 766 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 767 vb(i1:i2,j1:j2,1:jpkm1) = vn(i1:i2,j1:j2,1:jpkm1) 768 ENDIF 769 ! 746 770 ENDIF 747 771 ! -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/NST/agrif_top_interp.F90
r10068 r12191 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_ n(iref,jref,jk)107 h_out(jk) = e3t_a(ji,jj,jk) 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 tra(i1:i2,j1:j2,1:jpk,jn)=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 tra(ibdy+1,jmin:jmax,1:jpkm1,jn) = z1 * ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)154 DO jk = 1, jpkm1155 DO jj = jmin,jmax156 IF( umask(ibdy-1,jj,jk) == 0._wp ) THEN157 tra(ibdy,jj,jk,jn) = tra(ibdy+1,jj,jk,jn) * tmask(ibdy,jj,jk)158 ELSE159 tra(ibdy,jj,jk,jn)=(z4*tra(ibdy+1,jj,jk,jn)+z3*tra(ibdy-1,jj,jk,jn))*tmask(ibdy,jj,jk)160 IF( un(ibdy-1,jj,jk) > 0._wp ) THEN161 tra(ibdy,jj,jk,jn)=( z6*tra(ibdy-1,jj,jk,jn)+z5*tra(ibdy+1,jj,jk,jn) &162 + z7*tra(ibdy-2,jj,jk,jn) ) * tmask(ibdy,jj,jk)163 ENDIF164 ENDIF165 END DO166 END DO167 ! Restore ghost points:168 tra(ibdy+1,jmin:jmax,1:jpkm1,jn) = ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy+1,jmin:jmax,1:jpkm1)169 END DO170 ENDIF171 !172 IF( northern_side ) THEN173 zrho = Agrif_Rhoy()174 z1 = ( zrho - 1._wp ) * 0.5_wp175 z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )176 z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )177 z7 = - ( zrho - 1._wp ) / ( zrho + 3._wp )178 z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7179 !180 jbdy = nlcj-nbghostcells181 DO jn = 1, jptra182 tra(imin:imax,jbdy+1,1:jpkm1,jn) = z1 * ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)183 DO jk = 1, jpkm1184 DO ji = imin,imax185 IF( vmask(ji,jbdy-1,jk) == 0._wp ) THEN186 tra(ji,jbdy,jk,jn) = tra(ji,jbdy+1,jk,jn) * tmask(ji,jbdy,jk)187 ELSE188 tra(ji,jbdy,jk,jn)=(z4*tra(ji,jbdy+1,jk,jn)+z3*tra(ji,jbdy-1,jk,jn))*tmask(ji,jbdy,jk)189 IF (vn(ji,jbdy-1,jk) > 0._wp ) THEN190 tra(ji,jbdy,jk,jn)=( z6*tra(ji,jbdy-1,jk,jn)+z5*tra(ji,jbdy+1,jk,jn) &191 + z7*tra(ji,jbdy-2,jk,jn) ) * tmask(ji,jbdy,jk)192 ENDIF193 ENDIF194 END DO195 END DO196 ! Restore ghost points:197 tra(imin:imax,jbdy+1,1:jpkm1,jn) = ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) * tmask(imin:imax,jbdy+1,1:jpkm1)198 END DO199 ENDIF200 !201 IF( western_side ) THEN202 zrho = Agrif_Rhox()203 z1 = ( zrho - 1._wp ) * 0.5_wp204 z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )205 z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )206 z7 = - ( zrho - 1._wp ) / ( zrho + 3._wp )207 z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7208 !209 ibdy = 1+nbghostcells210 DO jn = 1, jptra211 tra(ibdy-1,jmin:jmax,1:jpkm1,jn) = z1 * ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)212 DO jk = 1, jpkm1213 DO jj = jmin,jmax214 IF( umask(ibdy,jj,jk) == 0._wp ) THEN215 tra(ibdy,jj,jk,jn) = tra(ibdy-1,jj,jk,jn) * tmask(ibdy,jj,jk)216 ELSE217 tra(ibdy,jj,jk,jn)=(z4*tra(ibdy-1,jj,jk,jn)+z3*tra(ibdy+1,jj,jk,jn))*tmask(ibdy,jj,jk)218 IF( un(ibdy,jj,jk) < 0._wp ) THEN219 tra(ibdy,jj,jk,jn)=( z6*tra(ibdy+1,jj,jk,jn)+z5*tra(ibdy-1,jj,jk,jn) &220 + z7*tra(ibdy+2,jj,jk,jn) ) * tmask(ibdy,jj,jk)221 ENDIF222 ENDIF223 END DO224 END DO225 ! Restore ghost points:226 tra(ibdy-1,jmin:jmax,1:jpkm1,jn) = ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy-1,jmin:jmax,1:jpkm1)227 END DO228 ENDIF229 !230 IF( southern_side ) THEN231 zrho = Agrif_Rhoy()232 z1 = ( zrho - 1._wp ) * 0.5_wp233 z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )234 z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )235 z7 = - ( zrho - 1._wp ) / ( zrho + 3._wp )236 z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7237 !238 jbdy=1+nbghostcells239 DO jn = 1, jptra240 tra(imin:imax,jbdy-1,1:jpkm1,jn) = z1 * ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)241 DO jk = 1, jpkm1242 DO ji = imin,imax243 IF( vmask(ji,jbdy,jk) == 0._wp ) THEN244 tra(ji,jbdy,jk,jn)=tra(ji,jbdy-1,jk,jn) * tmask(ji,jbdy,jk)245 ELSE246 tra(ji,jbdy,jk,jn)=(z4*tra(ji,jbdy-1,jk,jn)+z3*tra(ji,jbdy+1,jk,jn))*tmask(ji,jbdy,jk)247 IF( vn(ji,jbdy,jk) < 0._wp ) THEN248 tra(ji,jbdy,jk,jn)=( z6*tra(ji,jbdy+1,jk,jn)+z5*tra(ji,jbdy-1,jk,jn) &249 + z7*tra(ji,jbdy+2,jk,jn) ) * tmask(ji,jbdy,jk)250 ENDIF251 ENDIF252 END DO253 END DO254 ! Restore ghost points:255 tra(imin:imax,jbdy-1,1:jpkm1,jn) = ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) * tmask(imin:imax,jbdy-1,1:jpkm1)256 END DO257 ENDIF258 !259 ENDIF260 121 261 122 ENDIF -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/NST/agrif_top_sponge.F90
r10068 r12191 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 … … 93 94 DO jj=j1,j2 94 95 DO ji=i1,i2 95 tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_ n(ji,jj,jk)96 tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_b(ji,jj,jk) 96 97 END DO 97 98 END DO … … 114 115 IF (tmask(ji,jj,jk) == 0) EXIT 115 116 N_out = N_out + 1 116 h_out(jk) = e3t_ n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above117 h_out(jk) = e3t_b(ji,jj,jk) !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 … … 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_n(ji,jj,jk) * umask(ji,jj,jk)146 zabe2 = fsaht_spv(ji,jj) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk)147 zabe1 = rn_sponge_tra * fspu(ji,jj) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) 148 zabe2 = rn_sponge_tra * fspv(ji,jj) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) * 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) ) … … 155 157 tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ( ztu(ji,jj) - ztu(ji-1,jj ) & 156 158 & + ztv(ji,jj) - ztv(ji ,jj-1) ) & 157 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 159 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) & 160 & - ztrelax * fspt(ji,jj) * trbdiff(ji,jj,jk,jn) 158 161 ENDIF 159 162 END DO -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/NST/agrif_top_update.F90
r11078 r12191 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 … … 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 trb(ji,jj,jk,jn) = trb(ji,jj,jk,jn) & 141 & + atfp * ( tabres_child(ji,jj,jk,jn) & 142 & - trn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 134 ztb = trb(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! fse3t_b prior update should be used 135 ztnu = tabres_child(ji,jj,jk,jn) * e3t_n(ji,jj,jk) 136 ztno = trn(ji,jj,jk,jn) * e3t_a(ji,jj,jk) 137 trb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) ) & 138 & * tmask(ji,jj,jk) / e3t_b(ji,jj,jk) 143 139 ENDIF 144 140 ENDDO … … 148 144 ENDIF 149 145 DO jn = 1,jptra 150 DO jk=1,jpk 146 DO jk=1,jpkm1 151 147 DO jj=j1,j2 152 148 DO ji=i1,i2 153 149 IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 154 trn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk)150 trn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) 155 151 END IF 156 152 END DO … … 158 154 END DO 159 155 END DO 156 ! 157 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 158 trb(i1:i2,j1:j2,1:jpkm1,1:jptra) = trn(i1:i2,j1:j2,1:jpkm1,1:jptra) 159 ENDIF 160 ! 161 160 162 ENDIF 161 163 ! -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/NST/agrif_user.F90
r11536 r12191 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 56 ! !* Agrif initialization 57 CALL agrif_nemo_init 58 CALL Agrif_InitValues_cont_dom 59 CALL Agrif_InitValues_cont 20 !!---------------------------------------------------------------------- 21 USE nemogcm 22 !!---------------------------------------------------------------------- 23 ! 24 CALL nemo_init !* Initializations of each fine grid 25 ! 26 ! !* Agrif initialization 27 CALL agrif_nemo_init 28 CALL Agrif_InitValues_cont_dom 29 CALL Agrif_InitValues_cont 60 30 # if defined key_top 61 CALL Agrif_InitValues_cont_top31 CALL Agrif_InitValues_cont_top 62 32 # endif 63 33 # if defined key_si3 64 CALL Agrif_InitValues_cont_ice 65 # endif 66 ! 67 END SUBROUTINE Agrif_initvalues 68 69 70 SUBROUTINE Agrif_InitValues_cont_dom 71 !!---------------------------------------------------------------------- 72 !! *** ROUTINE Agrif_InitValues_cont *** 73 !! 74 !! ** Purpose :: Declaration of variables to be interpolated 75 !!---------------------------------------------------------------------- 76 USE Agrif_Util 77 USE oce 78 USE dom_oce 79 USE nemogcm 80 USE in_out_manager 81 USE agrif_oce_update 82 USE agrif_oce_interp 83 USE agrif_oce_sponge 84 ! 85 IMPLICIT NONE 86 !!---------------------------------------------------------------------- 87 ! 88 ! Declaration of the type of variable which have to be interpolated 89 ! 90 CALL agrif_declare_var_dom 91 ! 92 END SUBROUTINE Agrif_InitValues_cont_dom 93 94 95 SUBROUTINE agrif_declare_var_dom 96 !!---------------------------------------------------------------------- 97 !! *** ROUTINE agrif_declare_var *** 98 !! 99 !! ** Purpose :: Declaration of variables to be interpolated 100 !!---------------------------------------------------------------------- 101 USE agrif_util 102 USE par_oce 103 USE oce 104 ! 105 IMPLICIT NONE 106 ! 107 INTEGER :: ind1, ind2, ind3 34 CALL Agrif_InitValues_cont_ice 35 # endif 36 ! 37 END SUBROUTINE Agrif_initvalues 38 39 SUBROUTINE Agrif_InitValues_cont_dom 40 !!---------------------------------------------------------------------- 41 !! *** ROUTINE Agrif_InitValues_cont_dom *** 42 !!---------------------------------------------------------------------- 43 ! 44 CALL agrif_declare_var_dom 45 ! 46 END SUBROUTINE Agrif_InitValues_cont_dom 47 48 SUBROUTINE agrif_declare_var_dom 49 !!---------------------------------------------------------------------- 50 !! *** ROUTINE agrif_declare_var_dom *** 51 !!---------------------------------------------------------------------- 52 USE par_oce, ONLY: nbghostcells 53 ! 54 IMPLICIT NONE 55 ! 56 INTEGER :: ind1, ind2, ind3 108 57 !!---------------------------------------------------------------------- 109 58 110 59 ! 1. Declaration of the type of variable which have to be interpolated 111 60 !--------------------------------------------------------------------- 112 ind1 = nbghostcells113 ind2 = 1 + nbghostcells114 ind3 = 2 + nbghostcells115 CALL agrif_declare_variable((/1,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),e1u_id)116 CALL agrif_declare_variable((/2,1/),(/ind3,ind2/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),e2v_id)61 ind1 = nbghostcells 62 ind2 = 1 + nbghostcells 63 ind3 = 2 + nbghostcells 64 CALL agrif_declare_variable((/1,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),e1u_id) 65 CALL agrif_declare_variable((/2,1/),(/ind3,ind2/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),e2v_id) 117 66 118 67 ! 2. Type of interpolation 119 68 !------------------------- 120 CALL Agrif_Set_bcinterp( e1u_id, interp1=Agrif_linear, interp2=AGRIF_ppm )121 CALL Agrif_Set_bcinterp( e2v_id, interp1=AGRIF_ppm , interp2=Agrif_linear )69 CALL Agrif_Set_bcinterp( e1u_id, interp1=Agrif_linear, interp2=AGRIF_ppm ) 70 CALL Agrif_Set_bcinterp( e2v_id, interp1=AGRIF_ppm , interp2=Agrif_linear ) 122 71 123 72 ! 3. Location of interpolation 124 73 !----------------------------- 125 CALL Agrif_Set_bc(e1u_id,(/0,ind1-1/))126 CALL Agrif_Set_bc(e2v_id,(/0,ind1-1/))74 CALL Agrif_Set_bc(e1u_id,(/0,ind1-1/)) 75 CALL Agrif_Set_bc(e2v_id,(/0,ind1-1/)) 127 76 128 77 ! 4. Update type 129 78 !--------------- 130 79 # if defined UPD_HIGH 131 CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Average, update2=Agrif_Update_Full_Weighting)132 CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Full_Weighting, update2=Agrif_Update_Average)80 CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Average, update2=Agrif_Update_Full_Weighting) 81 CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Full_Weighting, update2=Agrif_Update_Average) 133 82 #else 134 CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Copy, update2=Agrif_Update_Average)135 CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Average, update2=Agrif_Update_Copy)83 CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Copy, update2=Agrif_Update_Average) 84 CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Average, update2=Agrif_Update_Copy) 136 85 #endif 137 86 138 END SUBROUTINE agrif_declare_var_dom 139 140 141 SUBROUTINE Agrif_InitValues_cont 87 END SUBROUTINE agrif_declare_var_dom 88 89 SUBROUTINE Agrif_InitValues_cont 142 90 !!---------------------------------------------------------------------- 143 91 !! *** ROUTINE Agrif_InitValues_cont *** 144 !! 145 !! ** Purpose :: Declaration of variables to be interpolated 146 !!---------------------------------------------------------------------- 147 USE agrif_oce_update 148 USE agrif_oce_interp 149 USE agrif_oce_sponge 150 USE Agrif_Util 151 USE oce 152 USE dom_oce 153 USE zdf_oce 154 USE nemogcm 155 ! 156 USE lib_mpp 157 USE in_out_manager 158 ! 159 IMPLICIT NONE 160 ! 161 LOGICAL :: check_namelist 162 CHARACTER(len=15) :: cl_check1, cl_check2, cl_check3, cl_check4 163 !!---------------------------------------------------------------------- 164 165 ! 1. Declaration of the type of variable which have to be interpolated 166 !--------------------------------------------------------------------- 167 CALL agrif_declare_var 168 169 ! 2. First interpolations of potentially non zero fields 170 !------------------------------------------------------- 171 Agrif_SpecialValue = 0._wp 172 Agrif_UseSpecialValue = .TRUE. 173 CALL Agrif_Bc_variable(tsn_id,calledweight=1.,procname=interptsn) 174 CALL Agrif_Sponge 175 tabspongedone_tsn = .FALSE. 176 CALL Agrif_Bc_variable(tsn_sponge_id,calledweight=1.,procname=interptsn_sponge) 177 ! reset tsa to zero 178 tsa(:,:,:,:) = 0. 179 180 Agrif_UseSpecialValue = ln_spc_dyn 181 CALL Agrif_Bc_variable(un_interp_id,calledweight=1.,procname=interpun) 182 CALL Agrif_Bc_variable(vn_interp_id,calledweight=1.,procname=interpvn) 183 tabspongedone_u = .FALSE. 184 tabspongedone_v = .FALSE. 185 CALL Agrif_Bc_variable(un_sponge_id,calledweight=1.,procname=interpun_sponge) 186 tabspongedone_u = .FALSE. 187 tabspongedone_v = .FALSE. 188 CALL Agrif_Bc_variable(vn_sponge_id,calledweight=1.,procname=interpvn_sponge) 189 190 Agrif_UseSpecialValue = .TRUE. 191 CALL Agrif_Bc_variable(sshn_id,calledweight=1., procname=interpsshn ) 192 hbdy_w(:,:) = 0.e0 ; hbdy_e(:,:) = 0.e0 ; hbdy_n(:,:) = 0.e0 ; hbdy_s(:,:) = 0.e0 193 ssha(:,:) = 0.e0 194 195 IF ( ln_dynspg_ts ) THEN 92 !!---------------------------------------------------------------------- 93 USE agrif_oce 94 USE agrif_oce_interp 95 USE agrif_oce_sponge 96 USE dom_oce 97 USE oce 98 USE lib_mpp 99 USE lbclnk 100 ! 101 IMPLICIT NONE 102 ! 103 INTEGER :: ji, jj 104 LOGICAL :: check_namelist 105 CHARACTER(len=15) :: cl_check1, cl_check2, cl_check3, cl_check4 106 #if defined key_vertical 107 REAL(wp), DIMENSION(jpi,jpj) :: zk ! workspace 108 #endif 109 !!---------------------------------------------------------------------- 110 111 ! 1. Declaration of the type of variable which have to be interpolated 112 !--------------------------------------------------------------------- 113 CALL agrif_declare_var 114 115 ! 2. First interpolations of potentially non zero fields 116 !------------------------------------------------------- 117 118 #if defined key_vertical 119 ! Build consistent parent bathymetry and number of levels 120 ! on the child grid 121 Agrif_UseSpecialValue = .FALSE. 122 ht0_parent(:,:) = 0._wp 123 mbkt_parent(:,:) = 0 124 ! 125 CALL Agrif_Bc_variable(ht0_id ,calledweight=1.,procname=interpht0 ) 126 CALL Agrif_Bc_variable(mbkt_id,calledweight=1.,procname=interpmbkt) 127 ! 128 ! Assume step wise change of bathymetry near interface 129 ! TODO: Switch to linear interpolation of bathymetry in the s-coordinate case 130 ! and no refinement 131 DO jj = 1, jpjm1 132 DO ji = 1, jpim1 133 mbku_parent(ji,jj) = MIN( mbkt_parent(ji+1,jj ) , mbkt_parent(ji,jj) ) 134 mbkv_parent(ji,jj) = MIN( mbkt_parent(ji ,jj+1) , mbkt_parent(ji,jj) ) 135 END DO 136 END DO 137 IF ( ln_sco.AND.Agrif_Parent(ln_sco) ) THEN 138 DO jj = 1, jpjm1 139 DO ji = 1, jpim1 140 hu0_parent(ji,jj) = 0.5_wp * ( ht0_parent(ji,jj)+ht0_parent(ji+1,jj) ) 141 hv0_parent(ji,jj) = 0.5_wp * ( ht0_parent(ji,jj)+ht0_parent(ji,jj+1) ) 142 END DO 143 END DO 144 ELSE 145 DO jj = 1, jpjm1 146 DO ji = 1, jpim1 147 hu0_parent(ji,jj) = MIN( ht0_parent(ji,jj), ht0_parent(ji+1,jj)) 148 hv0_parent(ji,jj) = MIN( ht0_parent(ji,jj), ht0_parent(ji,jj+1)) 149 END DO 150 END DO 151 152 ENDIF 153 ! 154 CALL lbc_lnk( 'Agrif_InitValues_cont', hu0_parent, 'U', 1. ) 155 CALL lbc_lnk( 'Agrif_InitValues_cont', hv0_parent, 'V', 1. ) 156 zk(:,:) = REAL( mbku_parent(:,:), wp ) ; CALL lbc_lnk( 'Agrif_InitValues_cont', zk, 'U', 1. ) 157 mbku_parent(:,:) = MAX( NINT( zk(:,:) ), 1 ) 158 zk(:,:) = REAL( mbkv_parent(:,:), wp ) ; CALL lbc_lnk( 'Agrif_InitValues_cont', zk, 'V', 1. ) 159 mbkv_parent(:,:) = MAX( NINT( zk(:,:) ), 1 ) 160 #endif 161 162 Agrif_SpecialValue = 0._wp 163 Agrif_UseSpecialValue = .TRUE. 164 CALL Agrif_Bc_variable(tsn_id,calledweight=1.,procname=interptsn) 165 CALL Agrif_Sponge 166 tabspongedone_tsn = .FALSE. 167 CALL Agrif_Bc_variable(tsn_sponge_id,calledweight=1.,procname=interptsn_sponge) 168 ! reset tsa to zero 169 tsa(:,:,:,:) = 0._wp 170 196 171 Agrif_UseSpecialValue = ln_spc_dyn 197 CALL Agrif_Bc_variable(unb_id,calledweight=1.,procname=interpunb) 198 CALL Agrif_Bc_variable(vnb_id,calledweight=1.,procname=interpvnb) 199 CALL Agrif_Bc_variable(ub2b_interp_id,calledweight=1.,procname=interpub2b) 200 CALL Agrif_Bc_variable(vb2b_interp_id,calledweight=1.,procname=interpvb2b) 201 ubdy_w(:,:) = 0.e0 ; vbdy_w(:,:) = 0.e0 202 ubdy_e(:,:) = 0.e0 ; vbdy_e(:,:) = 0.e0 203 ubdy_n(:,:) = 0.e0 ; vbdy_n(:,:) = 0.e0 204 ubdy_s(:,:) = 0.e0 ; vbdy_s(:,:) = 0.e0 205 ENDIF 206 207 Agrif_UseSpecialValue = .FALSE. 208 ! reset velocities to zero 209 ua(:,:,:) = 0. 210 va(:,:,:) = 0. 211 212 ! 3. Some controls 213 !----------------- 214 check_namelist = .TRUE. 215 216 IF( check_namelist ) THEN 217 218 ! Check time steps 219 IF( NINT(Agrif_Rhot()) * NINT(rdt) .NE. Agrif_Parent(rdt) ) THEN 220 WRITE(cl_check1,*) NINT(Agrif_Parent(rdt)) 221 WRITE(cl_check2,*) NINT(rdt) 222 WRITE(cl_check3,*) NINT(Agrif_Parent(rdt)/Agrif_Rhot()) 223 CALL ctl_stop( 'Incompatible time step between ocean grids', & 224 & 'parent grid value : '//cl_check1 , & 225 & 'child grid value : '//cl_check2 , & 226 & 'value on child grid should be changed to : '//cl_check3 ) 227 ENDIF 228 229 ! Check run length 230 IF( Agrif_IRhot() * (Agrif_Parent(nitend)- & 231 Agrif_Parent(nit000)+1) .NE. (nitend-nit000+1) ) THEN 232 WRITE(cl_check1,*) (Agrif_Parent(nit000)-1)*Agrif_IRhot() + 1 233 WRITE(cl_check2,*) Agrif_Parent(nitend) *Agrif_IRhot() 234 CALL ctl_warn( 'Incompatible run length between grids' , & 235 & 'nit000 on fine grid will be changed to : '//cl_check1, & 236 & 'nitend on fine grid will be changed to : '//cl_check2 ) 237 nit000 = (Agrif_Parent(nit000)-1)*Agrif_IRhot() + 1 238 nitend = Agrif_Parent(nitend) *Agrif_IRhot() 239 ENDIF 240 241 ! Check free surface scheme 242 IF ( ( Agrif_Parent(ln_dynspg_ts ).AND.ln_dynspg_exp ).OR.& 243 & ( Agrif_Parent(ln_dynspg_exp).AND.ln_dynspg_ts ) ) THEN 244 WRITE(cl_check1,*) Agrif_Parent( ln_dynspg_ts ) 245 WRITE(cl_check2,*) ln_dynspg_ts 246 WRITE(cl_check3,*) Agrif_Parent( ln_dynspg_exp ) 247 WRITE(cl_check4,*) ln_dynspg_exp 248 CALL ctl_stop( 'Incompatible free surface scheme between grids' , & 249 & 'parent grid ln_dynspg_ts :'//cl_check1 , & 250 & 'child grid ln_dynspg_ts :'//cl_check2 , & 251 & 'parent grid ln_dynspg_exp :'//cl_check3 , & 252 & 'child grid ln_dynspg_exp :'//cl_check4 , & 253 & 'those logicals should be identical' ) 254 STOP 255 ENDIF 256 257 ! Check if identical linear free surface option 258 IF ( ( Agrif_Parent(ln_linssh ).AND.(.NOT.ln_linssh )).OR.& 259 & ( (.NOT.Agrif_Parent(ln_linssh)).AND.ln_linssh ) ) THEN 260 WRITE(cl_check1,*) Agrif_Parent(ln_linssh ) 261 WRITE(cl_check2,*) ln_linssh 262 CALL ctl_stop( 'Incompatible linearized fs option between grids', & 263 & 'parent grid ln_linssh :'//cl_check1 , & 264 & 'child grid ln_linssh :'//cl_check2 , & 265 & 'those logicals should be identical' ) 266 STOP 172 CALL Agrif_Bc_variable(un_interp_id,calledweight=1.,procname=interpun) 173 CALL Agrif_Bc_variable(vn_interp_id,calledweight=1.,procname=interpvn) 174 tabspongedone_u = .FALSE. 175 tabspongedone_v = .FALSE. 176 CALL Agrif_Bc_variable(un_sponge_id,calledweight=1.,procname=interpun_sponge) 177 tabspongedone_u = .FALSE. 178 tabspongedone_v = .FALSE. 179 CALL Agrif_Bc_variable(vn_sponge_id,calledweight=1.,procname=interpvn_sponge) 180 ua(:,:,:) = 0._wp 181 va(:,:,:) = 0._wp 182 183 Agrif_UseSpecialValue = .TRUE. 184 CALL Agrif_Bc_variable(sshn_id,calledweight=1., procname=interpsshn ) 185 hbdy(:,:) = 0._wp 186 ssha(:,:) = 0._wp 187 188 IF ( ln_dynspg_ts ) THEN 189 Agrif_UseSpecialValue = ln_spc_dyn 190 CALL Agrif_Bc_variable(unb_id,calledweight=1.,procname=interpunb) 191 CALL Agrif_Bc_variable(vnb_id,calledweight=1.,procname=interpvnb) 192 CALL Agrif_Bc_variable(ub2b_interp_id,calledweight=1.,procname=interpub2b) 193 CALL Agrif_Bc_variable(vb2b_interp_id,calledweight=1.,procname=interpvb2b) 194 ubdy(:,:) = 0._wp 195 vbdy(:,:) = 0._wp 196 ENDIF 197 198 Agrif_UseSpecialValue = .FALSE. 199 200 ! 3. Some controls 201 !----------------- 202 check_namelist = .TRUE. 203 204 IF( check_namelist ) THEN 205 206 ! Check time steps 207 IF( NINT(Agrif_Rhot()) * NINT(rdt) .NE. Agrif_Parent(rdt) ) THEN 208 WRITE(cl_check1,*) NINT(Agrif_Parent(rdt)) 209 WRITE(cl_check2,*) NINT(rdt) 210 WRITE(cl_check3,*) NINT(Agrif_Parent(rdt)/Agrif_Rhot()) 211 CALL ctl_stop( 'Incompatible time step between ocean grids', & 212 & 'parent grid value : '//cl_check1 , & 213 & 'child grid value : '//cl_check2 , & 214 & 'value on child grid should be changed to : '//cl_check3 ) 215 ENDIF 216 217 ! Check run length 218 IF( Agrif_IRhot() * (Agrif_Parent(nitend)- & 219 Agrif_Parent(nit000)+1) .NE. (nitend-nit000+1) ) THEN 220 WRITE(cl_check1,*) (Agrif_Parent(nit000)-1)*Agrif_IRhot() + 1 221 WRITE(cl_check2,*) Agrif_Parent(nitend) *Agrif_IRhot() 222 CALL ctl_warn( 'Incompatible run length between grids' , & 223 & 'nit000 on fine grid will be changed to : '//cl_check1, & 224 & 'nitend on fine grid will be changed to : '//cl_check2 ) 225 nit000 = (Agrif_Parent(nit000)-1)*Agrif_IRhot() + 1 226 nitend = Agrif_Parent(nitend) *Agrif_IRhot() 227 ENDIF 228 229 ! Check free surface scheme 230 IF ( ( Agrif_Parent(ln_dynspg_ts ).AND.ln_dynspg_exp ).OR.& 231 & ( Agrif_Parent(ln_dynspg_exp).AND.ln_dynspg_ts ) ) THEN 232 WRITE(cl_check1,*) Agrif_Parent( ln_dynspg_ts ) 233 WRITE(cl_check2,*) ln_dynspg_ts 234 WRITE(cl_check3,*) Agrif_Parent( ln_dynspg_exp ) 235 WRITE(cl_check4,*) ln_dynspg_exp 236 CALL ctl_stop( 'Incompatible free surface scheme between grids' , & 237 & 'parent grid ln_dynspg_ts :'//cl_check1 , & 238 & 'child grid ln_dynspg_ts :'//cl_check2 , & 239 & 'parent grid ln_dynspg_exp :'//cl_check3 , & 240 & 'child grid ln_dynspg_exp :'//cl_check4 , & 241 & 'those logicals should be identical' ) 242 STOP 243 ENDIF 244 245 ! Check if identical linear free surface option 246 IF ( ( Agrif_Parent(ln_linssh ).AND.(.NOT.ln_linssh )).OR.& 247 & ( (.NOT.Agrif_Parent(ln_linssh)).AND.ln_linssh ) ) THEN 248 WRITE(cl_check1,*) Agrif_Parent(ln_linssh ) 249 WRITE(cl_check2,*) ln_linssh 250 CALL ctl_stop( 'Incompatible linearized fs option between grids', & 251 & 'parent grid ln_linssh :'//cl_check1 , & 252 & 'child grid ln_linssh :'//cl_check2 , & 253 & 'those logicals should be identical' ) 254 STOP 255 ENDIF 256 267 257 ENDIF 268 258 269 259 ! check if masks and bathymetries match 270 260 IF(ln_chk_bathy) THEN 261 Agrif_UseSpecialValue = .FALSE. 271 262 ! 263 IF(lwp) WRITE(numout,*) ' ' 272 264 IF(lwp) WRITE(numout,*) 'AGRIF: Check Bathymetry and masks near bdys. Level: ', Agrif_Level() 273 265 ! 274 266 kindic_agr = 0 275 ! check if umask agree with parent along western and eastern boundaries: 276 CALL Agrif_Bc_variable(umsk_id,calledweight=1.,procname=interpumsk) 277 ! check if vmask agree with parent along northern and southern boundaries: 278 CALL Agrif_Bc_variable(vmsk_id,calledweight=1.,procname=interpvmsk) 279 ! check if tmask and vertical scale factors agree with parent over first two coarse grid points: 267 # if ! defined key_vertical 268 ! 269 ! check if tmask and vertical scale factors agree with parent in sponge area: 280 270 CALL Agrif_Bc_variable(e3t_id,calledweight=1.,procname=interpe3t) 281 271 ! 272 # else 273 ! 274 ! In case of vertical interpolation, check only that total depths agree between child and parent: 275 DO ji = 1, jpi 276 DO jj = 1, jpj 277 IF ((mbkt_parent(ji,jj)/=0).AND.(ABS(ht0_parent(ji,jj)-ht_0(ji,jj))>1.e-3)) kindic_agr = kindic_agr + 1 278 IF ((mbku_parent(ji,jj)/=0).AND.(ABS(hu0_parent(ji,jj)-hu_0(ji,jj))>1.e-3)) kindic_agr = kindic_agr + 1 279 IF ((mbkv_parent(ji,jj)/=0).AND.(ABS(hv0_parent(ji,jj)-hv_0(ji,jj))>1.e-3)) kindic_agr = kindic_agr + 1 280 END DO 281 END DO 282 # endif 282 283 CALL mpp_sum( 'agrif_user', kindic_agr ) 283 284 IF( kindic_agr /= 0 ) THEN 284 CALL ctl_stop(' Child Bathymetry is notcorrect near boundaries.')285 CALL ctl_stop('==> Child Bathymetry is NOT correct near boundaries.') 285 286 ELSE 286 IF(lwp) WRITE(numout,*) 'Child Bathymetry is ok near boundaries.' 287 END IF 288 ENDIF 289 ! 290 ENDIF 291 ! 292 END SUBROUTINE Agrif_InitValues_cont 293 294 SUBROUTINE agrif_declare_var 295 !!---------------------------------------------------------------------- 296 !! *** ROUTINE agrif_declarE_var *** 297 !! 298 !! ** Purpose :: Declaration of variables to be interpolated 299 !!---------------------------------------------------------------------- 300 USE agrif_util 301 USE agrif_oce 302 USE par_oce ! ocean parameters 303 USE zdf_oce ! vertical physics 304 USE oce 305 ! 306 IMPLICIT NONE 307 ! 308 INTEGER :: ind1, ind2, ind3 309 !!---------------------------------------------------------------------- 310 311 ! 1. Declaration of the type of variable which have to be interpolated 312 !--------------------------------------------------------------------- 313 ind1 = nbghostcells 314 ind2 = 1 + nbghostcells 315 ind3 = 2 + nbghostcells 287 IF(lwp) WRITE(numout,*) '==> Child Bathymetry is ok near boundaries.' 288 IF(lwp) WRITE(numout,*) ' ' 289 END IF 290 ! 291 ENDIF 292 316 293 # if defined key_vertical 317 CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsn_id) 318 CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsn_sponge_id) 319 320 CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_interp_id) 321 CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_interp_id) 322 CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_update_id) 323 CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_update_id) 324 CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_sponge_id) 325 CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_sponge_id) 294 ! Additional constrain that should be removed someday: 295 IF ( Agrif_Parent(jpk).GT.jpk ) THEN 296 CALL ctl_stop( ' With key_vertical, child grids must have jpk greater or equal to the parent value' ) 297 ENDIF 298 # endif 299 ! 300 END SUBROUTINE Agrif_InitValues_cont 301 302 SUBROUTINE agrif_declare_var 303 !!---------------------------------------------------------------------- 304 !! *** ROUTINE agrif_declare_var *** 305 !!---------------------------------------------------------------------- 306 USE agrif_util 307 USE agrif_oce 308 USE par_oce 309 USE zdf_oce 310 USE oce 311 ! 312 IMPLICIT NONE 313 ! 314 INTEGER :: ind1, ind2, ind3 315 !!---------------------------------------------------------------------- 316 317 ! 1. Declaration of the type of variable which have to be interpolated 318 !--------------------------------------------------------------------- 319 ind1 = nbghostcells 320 ind2 = 1 + nbghostcells 321 ind3 = 2 + nbghostcells 322 # if defined key_vertical 323 CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsn_id) 324 CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsn_sponge_id) 325 326 CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_interp_id) 327 CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_interp_id) 328 CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_update_id) 329 CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_update_id) 330 CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_sponge_id) 331 CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_sponge_id) 326 332 # else 327 CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts/),tsn_id) 328 CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts/),tsn_sponge_id) 329 330 CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),un_interp_id) 331 CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),vn_interp_id) 332 CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),un_update_id) 333 CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),vn_update_id) 334 CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),un_sponge_id) 335 CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),vn_sponge_id) 336 # endif 337 338 CALL agrif_declare_variable((/2,2,0/),(/ind3,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),e3t_id) 339 CALL agrif_declare_variable((/1,2,0/),(/ind2,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),umsk_id) 340 CALL agrif_declare_variable((/2,1,0/),(/ind3,ind2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vmsk_id) 341 342 CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,3/),scales_t_id) 343 344 CALL agrif_declare_variable((/1,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),unb_id) 345 CALL agrif_declare_variable((/2,1/),(/ind3,ind2/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),vnb_id) 346 CALL agrif_declare_variable((/1,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),ub2b_interp_id) 347 CALL agrif_declare_variable((/2,1/),(/ind3,ind2/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),vb2b_interp_id) 348 CALL agrif_declare_variable((/1,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),ub2b_update_id) 349 CALL agrif_declare_variable((/2,1/),(/ind3,ind2/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),vb2b_update_id) 350 351 CALL agrif_declare_variable((/2,2/),(/ind3,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),sshn_id) 352 353 IF( ln_zdftke.OR.ln_zdfgls ) THEN 354 ! CALL agrif_declare_variable((/2,2,0/),(/ind3,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/), en_id) 355 ! CALL agrif_declare_variable((/2,2,0/),(/ind3,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),avt_id) 333 CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts/),tsn_id) 334 CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts/),tsn_sponge_id) 335 336 CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),un_interp_id) 337 CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),vn_interp_id) 338 CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),un_update_id) 339 CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),vn_update_id) 340 CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),un_sponge_id) 341 CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),vn_sponge_id) 342 # endif 343 344 CALL agrif_declare_variable((/2,2,0/),(/ind3,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),e3t_id) 345 356 346 # if defined key_vertical 357 CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),avm_id) 347 CALL agrif_declare_variable((/2,2/),(/ind3,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),mbkt_id) 348 CALL agrif_declare_variable((/2,2/),(/ind3,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),ht0_id) 349 # endif 350 351 CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,3/),scales_t_id) 352 353 CALL agrif_declare_variable((/1,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),unb_id) 354 CALL agrif_declare_variable((/2,1/),(/ind3,ind2/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),vnb_id) 355 CALL agrif_declare_variable((/1,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),ub2b_interp_id) 356 CALL agrif_declare_variable((/2,1/),(/ind3,ind2/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),vb2b_interp_id) 357 CALL agrif_declare_variable((/1,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),ub2b_update_id) 358 CALL agrif_declare_variable((/2,1/),(/ind3,ind2/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),vb2b_update_id) 359 360 CALL agrif_declare_variable((/2,2/),(/ind3,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),sshn_id) 361 362 IF( ln_zdftke.OR.ln_zdfgls ) THEN 363 ! CALL agrif_declare_variable((/2,2,0/),(/ind3,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/), en_id) 364 ! CALL agrif_declare_variable((/2,2,0/),(/ind3,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),avt_id) 365 # if defined key_vertical 366 CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),avm_id) 358 367 # else 359 CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),avm_id) 360 # endif 361 ENDIF 362 363 ! 2. Type of interpolation 364 !------------------------- 365 CALL Agrif_Set_bcinterp(tsn_id,interp=AGRIF_linear) 366 367 CALL Agrif_Set_bcinterp(un_interp_id,interp1=Agrif_linear,interp2=AGRIF_ppm) 368 CALL Agrif_Set_bcinterp(vn_interp_id,interp1=AGRIF_ppm,interp2=Agrif_linear) 369 370 CALL Agrif_Set_bcinterp(tsn_sponge_id,interp=AGRIF_linear) 371 372 CALL Agrif_Set_bcinterp(sshn_id,interp=AGRIF_linear) 373 CALL Agrif_Set_bcinterp(unb_id,interp1=Agrif_linear,interp2=AGRIF_ppm) 374 CALL Agrif_Set_bcinterp(vnb_id,interp1=AGRIF_ppm,interp2=Agrif_linear) 375 CALL Agrif_Set_bcinterp(ub2b_interp_id,interp1=Agrif_linear,interp2=AGRIF_ppm) 376 CALL Agrif_Set_bcinterp(vb2b_interp_id,interp1=AGRIF_ppm,interp2=Agrif_linear) 377 378 379 CALL Agrif_Set_bcinterp(un_sponge_id,interp1=Agrif_linear,interp2=AGRIF_ppm) 380 CALL Agrif_Set_bcinterp(vn_sponge_id,interp1=AGRIF_ppm,interp2=Agrif_linear) 381 382 CALL Agrif_Set_bcinterp(e3t_id,interp=AGRIF_constant) 383 CALL Agrif_Set_bcinterp(umsk_id,interp=AGRIF_constant) 384 CALL Agrif_Set_bcinterp(vmsk_id,interp=AGRIF_constant) 385 386 IF( ln_zdftke.OR.ln_zdfgls ) CALL Agrif_Set_bcinterp( avm_id, interp=AGRIF_linear ) 387 388 ! 3. Location of interpolation 389 !----------------------------- 390 CALL Agrif_Set_bc( tsn_id, (/0,ind1/) ) 391 CALL Agrif_Set_bc( un_interp_id, (/0,ind1/) ) 392 CALL Agrif_Set_bc( vn_interp_id, (/0,ind1/) ) 393 394 CALL Agrif_Set_bc( tsn_sponge_id, (/-nn_sponge_len*Agrif_irhox()-1,0/) ) ! if west and rhox=3 and sponge=2 and ghost=1: columns 2 to 9 395 CALL Agrif_Set_bc( un_sponge_id, (/-nn_sponge_len*Agrif_irhox()-1,0/) ) 396 CALL Agrif_Set_bc( vn_sponge_id, (/-nn_sponge_len*Agrif_irhox()-1,0/) ) 397 398 CALL Agrif_Set_bc( sshn_id, (/0,ind1-1/) ) 399 CALL Agrif_Set_bc( unb_id, (/0,ind1-1/) ) 400 CALL Agrif_Set_bc( vnb_id, (/0,ind1-1/) ) 401 CALL Agrif_Set_bc( ub2b_interp_id, (/0,ind1-1/) ) 402 CALL Agrif_Set_bc( vb2b_interp_id, (/0,ind1-1/) ) 403 404 CALL Agrif_Set_bc( e3t_id, (/-nn_sponge_len*Agrif_irhox(),ind1-1/) ) ! if west and rhox=3 and ghost=1: column 2 to 6 405 CALL Agrif_Set_bc( umsk_id, (/0,0/) ) 406 CALL Agrif_Set_bc( vmsk_id, (/0,0/) ) 407 408 409 IF( ln_zdftke.OR.ln_zdfgls ) CALL Agrif_Set_bc( avm_id, (/0,ind1/) ) 410 411 ! 4. Update type 412 !--------------- 413 CALL Agrif_Set_Updatetype(scales_t_id, update = AGRIF_Update_Average) 368 CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),avm_id) 369 # endif 370 ENDIF 371 372 ! 2. Type of interpolation 373 !------------------------- 374 CALL Agrif_Set_bcinterp(tsn_id,interp=AGRIF_linear) 375 376 CALL Agrif_Set_bcinterp(un_interp_id,interp1=Agrif_linear,interp2=AGRIF_ppm) 377 CALL Agrif_Set_bcinterp(vn_interp_id,interp1=AGRIF_ppm,interp2=Agrif_linear) 378 379 CALL Agrif_Set_bcinterp(tsn_sponge_id,interp=AGRIF_linear) 380 381 CALL Agrif_Set_bcinterp(sshn_id,interp=AGRIF_linear) 382 CALL Agrif_Set_bcinterp(unb_id,interp1=Agrif_linear,interp2=AGRIF_ppm) 383 CALL Agrif_Set_bcinterp(vnb_id,interp1=AGRIF_ppm,interp2=Agrif_linear) 384 CALL Agrif_Set_bcinterp(ub2b_interp_id,interp1=Agrif_linear,interp2=AGRIF_ppm) 385 CALL Agrif_Set_bcinterp(vb2b_interp_id,interp1=AGRIF_ppm,interp2=Agrif_linear) 386 ! 387 ! > Divergence conserving alternative: 388 ! CALL Agrif_Set_bcinterp(sshn_id,interp=AGRIF_constant) 389 ! CALL Agrif_Set_bcinterp(unb_id,interp1=Agrif_linear,interp2=AGRIF_constant) 390 ! CALL Agrif_Set_bcinterp(vnb_id,interp1=AGRIF_constant,interp2=Agrif_linear) 391 ! CALL Agrif_Set_bcinterp(ub2b_interp_id,interp1=Agrif_linear,interp2=AGRIF_constant) 392 ! CALL Agrif_Set_bcinterp(vb2b_interp_id,interp1=AGRIF_constant,interp2=Agrif_linear) 393 !< 394 395 CALL Agrif_Set_bcinterp(un_sponge_id,interp1=Agrif_linear,interp2=AGRIF_ppm) 396 CALL Agrif_Set_bcinterp(vn_sponge_id,interp1=AGRIF_ppm,interp2=Agrif_linear) 397 398 CALL Agrif_Set_bcinterp(e3t_id,interp=AGRIF_constant) 399 400 # if defined key_vertical 401 CALL Agrif_Set_bcinterp(mbkt_id,interp=AGRIF_constant) 402 CALL Agrif_Set_bcinterp(ht0_id ,interp=AGRIF_constant) 403 # endif 404 405 IF( ln_zdftke.OR.ln_zdfgls ) CALL Agrif_Set_bcinterp( avm_id, interp=AGRIF_linear ) 406 407 ! 3. Location of interpolation 408 !----------------------------- 409 CALL Agrif_Set_bc( tsn_id, (/0,ind1-1/) ) ! if west, rhox=3 and nbghost=3: columns 2 to 4 410 CALL Agrif_Set_bc( un_interp_id, (/0,ind1-1/) ) 411 CALL Agrif_Set_bc( vn_interp_id, (/0,ind1-1/) ) 412 413 CALL Agrif_Set_bc( tsn_sponge_id, (/-nn_sponge_len*Agrif_irhox()-1,0/) ) ! if west, rhox=3, nn_sponge_len=2 414 CALL Agrif_Set_bc( un_sponge_id, (/-nn_sponge_len*Agrif_irhox()-1,0/) ) ! and nbghost=3: 415 CALL Agrif_Set_bc( vn_sponge_id, (/-nn_sponge_len*Agrif_irhox()-1,0/) ) ! columns 4 to 11 416 417 CALL Agrif_Set_bc( sshn_id, (/0,ind1-1/) ) 418 CALL Agrif_Set_bc( unb_id, (/0,ind1-1/) ) 419 CALL Agrif_Set_bc( vnb_id, (/0,ind1-1/) ) 420 CALL Agrif_Set_bc( ub2b_interp_id, (/0,ind1-1/) ) 421 CALL Agrif_Set_bc( vb2b_interp_id, (/0,ind1-1/) ) 422 423 ! CALL Agrif_Set_bc( e3t_id, (/-nn_sponge_len*Agrif_irhox(),ind1-1/) ) 424 ! JC: check near the boundary only until matching in sponge has been sorted out: 425 CALL Agrif_Set_bc( e3t_id, (/0,ind1-1/) ) 426 427 # if defined key_vertical 428 ! extend the interpolation zone by 1 more point than necessary: 429 CALL Agrif_Set_bc( mbkt_id, (/-nn_sponge_len*Agrif_irhox()-2,ind1/) ) 430 CALL Agrif_Set_bc( ht0_id, (/-nn_sponge_len*Agrif_irhox()-2,ind1/) ) 431 # endif 432 433 IF( ln_zdftke.OR.ln_zdfgls ) CALL Agrif_Set_bc( avm_id, (/0,ind1/) ) 434 435 ! 4. Update type 436 !--------------- 437 CALL Agrif_Set_Updatetype(scales_t_id, update = AGRIF_Update_Average) 414 438 415 439 # if defined UPD_HIGH 416 CALL Agrif_Set_Updatetype(tsn_id, update = Agrif_Update_Full_Weighting)417 CALL Agrif_Set_Updatetype(un_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Full_Weighting)418 CALL Agrif_Set_Updatetype(vn_update_id,update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average)419 420 CALL Agrif_Set_Updatetype(ub2b_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Full_Weighting)421 CALL Agrif_Set_Updatetype(vb2b_update_id,update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average)422 CALL Agrif_Set_Updatetype(sshn_id,update = Agrif_Update_Full_Weighting)423 CALL Agrif_Set_Updatetype(e3t_id, update = Agrif_Update_Full_Weighting)424 425 IF( ln_zdftke.OR.ln_zdfgls ) THEN426 ! CALL Agrif_Set_Updatetype( en_id, update = AGRIF_Update_Full_Weighting)427 ! CALL Agrif_Set_Updatetype(avt_id, update = AGRIF_Update_Full_Weighting)428 ! CALL Agrif_Set_Updatetype(avm_id, update = AGRIF_Update_Full_Weighting)429 ENDIF440 CALL Agrif_Set_Updatetype(tsn_id, update = Agrif_Update_Full_Weighting) 441 CALL Agrif_Set_Updatetype(un_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Full_Weighting) 442 CALL Agrif_Set_Updatetype(vn_update_id,update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average) 443 444 CALL Agrif_Set_Updatetype(ub2b_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Full_Weighting) 445 CALL Agrif_Set_Updatetype(vb2b_update_id,update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average) 446 CALL Agrif_Set_Updatetype(sshn_id,update = Agrif_Update_Full_Weighting) 447 CALL Agrif_Set_Updatetype(e3t_id, update = Agrif_Update_Full_Weighting) 448 449 IF( ln_zdftke.OR.ln_zdfgls ) THEN 450 ! CALL Agrif_Set_Updatetype( en_id, update = AGRIF_Update_Full_Weighting) 451 ! CALL Agrif_Set_Updatetype(avt_id, update = AGRIF_Update_Full_Weighting) 452 ! CALL Agrif_Set_Updatetype(avm_id, update = AGRIF_Update_Full_Weighting) 453 ENDIF 430 454 431 455 #else 432 CALL Agrif_Set_Updatetype(tsn_id, update = AGRIF_Update_Average)433 CALL Agrif_Set_Updatetype(un_update_id,update1 = Agrif_Update_Copy, update2 = Agrif_Update_Average)434 CALL Agrif_Set_Updatetype(vn_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy)435 436 CALL Agrif_Set_Updatetype(ub2b_update_id,update1 = Agrif_Update_Copy, update2 = Agrif_Update_Average)437 CALL Agrif_Set_Updatetype(vb2b_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy)438 CALL Agrif_Set_Updatetype(sshn_id,update = AGRIF_Update_Average)439 CALL Agrif_Set_Updatetype(e3t_id, update = AGRIF_Update_Average)440 441 IF( ln_zdftke.OR.ln_zdfgls ) THEN442 ! CALL Agrif_Set_Updatetype( en_id, update = AGRIF_Update_Average)443 ! CALL Agrif_Set_Updatetype(avt_id, update = AGRIF_Update_Average)444 ! CALL Agrif_Set_Updatetype(avm_id, update = AGRIF_Update_Average)445 ENDIF456 CALL Agrif_Set_Updatetype(tsn_id, update = AGRIF_Update_Average) 457 CALL Agrif_Set_Updatetype(un_update_id,update1 = Agrif_Update_Copy, update2 = Agrif_Update_Average) 458 CALL Agrif_Set_Updatetype(vn_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy) 459 460 CALL Agrif_Set_Updatetype(ub2b_update_id,update1 = Agrif_Update_Copy, update2 = Agrif_Update_Average) 461 CALL Agrif_Set_Updatetype(vb2b_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy) 462 CALL Agrif_Set_Updatetype(sshn_id,update = AGRIF_Update_Average) 463 CALL Agrif_Set_Updatetype(e3t_id, update = AGRIF_Update_Average) 464 465 IF( ln_zdftke.OR.ln_zdfgls ) THEN 466 ! CALL Agrif_Set_Updatetype( en_id, update = AGRIF_Update_Average) 467 ! CALL Agrif_Set_Updatetype(avt_id, update = AGRIF_Update_Average) 468 ! CALL Agrif_Set_Updatetype(avm_id, update = AGRIF_Update_Average) 469 ENDIF 446 470 447 471 #endif 448 !449 END SUBROUTINE agrif_declare_var472 ! 473 END SUBROUTINE agrif_declare_var 450 474 451 475 #if defined key_si3 … … 453 477 !!---------------------------------------------------------------------- 454 478 !! *** ROUTINE Agrif_InitValues_cont_ice *** 479 !!---------------------------------------------------------------------- 480 USE Agrif_Util 481 USE sbc_oce, ONLY : nn_fsbc ! clem: necessary otherwise Agrif_Parent(nn_fsbc) = nn_fsbc 482 USE ice 483 USE agrif_ice 484 USE in_out_manager 485 USE agrif_ice_interp 486 USE lib_mpp 487 ! 488 IMPLICIT NONE 489 !!---------------------------------------------------------------------- 490 ! 491 ! Declaration of the type of variable which have to be interpolated (parent=>child) 492 !---------------------------------------------------------------------------------- 493 CALL agrif_declare_var_ice 494 495 ! Controls 496 497 ! clem: For some reason, nn_fsbc(child)/=1 does not work properly (signal can be largely degraded by the agrif zoom) 498 ! the run must satisfy CFL=Uice/(dx/dt) < 0.6/nn_fsbc(child) 499 ! therefore, if nn_fsbc(child)>1 one must reduce the time-step in proportion to nn_fsbc(child), which is not acceptable 500 ! If a solution is found, the following stop could be removed because the rest of the code take nn_fsbc(child) into account 501 IF( nn_fsbc > 1 ) CALL ctl_stop('nn_fsbc(child) must be set to 1 otherwise agrif and sea-ice may not work properly') 502 503 ! stop if rhot * nn_fsbc(parent) /= N * nn_fsbc(child) with N being integer 504 IF( MOD( Agrif_irhot() * Agrif_Parent(nn_fsbc), nn_fsbc ) /= 0 ) THEN 505 CALL ctl_stop('rhot * nn_fsbc(parent) /= N * nn_fsbc(child), therefore nn_fsbc(child) should be set to 1 or nn_fsbc(parent)') 506 ENDIF 507 ! First Interpolations (using "after" ice subtime step => nbstep_ice=1) 508 !---------------------------------------------------------------------- 509 nbstep_ice = ( Agrif_irhot() * Agrif_Parent(nn_fsbc) / nn_fsbc ) ! clem: to have calledweight=1 in interp (otherwise the western border of the zoom is wrong) 510 CALL agrif_interp_ice('U') ! interpolation of ice velocities 511 CALL agrif_interp_ice('V') ! interpolation of ice velocities 512 CALL agrif_interp_ice('T') ! interpolation of ice tracers 513 nbstep_ice = 0 514 ! 515 END SUBROUTINE Agrif_InitValues_cont_ice 516 517 SUBROUTINE agrif_declare_var_ice 518 !!---------------------------------------------------------------------- 519 !! *** ROUTINE agrif_declare_var_ice *** 520 !!---------------------------------------------------------------------- 521 USE Agrif_Util 522 USE ice 523 USE par_oce, ONLY : nbghostcells 524 ! 525 IMPLICIT NONE 526 ! 527 INTEGER :: ind1, ind2, ind3 528 !!---------------------------------------------------------------------- 529 ! 530 ! 1. Declaration of the type of variable which have to be interpolated (parent=>child) 531 ! agrif_declare_variable(position,1st point index,--,--,dimensions,name) 532 ! ex.: position=> 1,1 = not-centered (in i and j) 533 ! 2,2 = centered ( - ) 534 ! index => 1,1 = one ghost line 535 ! 2,2 = two ghost lines 536 !------------------------------------------------------------------------------------- 537 ind1 = nbghostcells 538 ind2 = 1 + nbghostcells 539 ind3 = 2 + nbghostcells 540 CALL agrif_declare_variable((/2,2,0/),(/ind3,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpl*(8+nlay_s+nlay_i)/),tra_ice_id) 541 CALL agrif_declare_variable((/1,2/) ,(/ind2,ind3/) ,(/'x','y'/) ,(/1,1/) ,(/nlci,nlcj/) ,u_ice_id ) 542 CALL agrif_declare_variable((/2,1/) ,(/ind3,ind2/) ,(/'x','y'/) ,(/1,1/) ,(/nlci,nlcj/) ,v_ice_id ) 543 544 ! 2. Set interpolations (normal & tangent to the grid cell for velocities) 545 !----------------------------------- 546 CALL Agrif_Set_bcinterp(tra_ice_id, interp = AGRIF_linear) 547 CALL Agrif_Set_bcinterp(u_ice_id , interp1 = Agrif_linear,interp2 = AGRIF_ppm ) 548 CALL Agrif_Set_bcinterp(v_ice_id , interp1 = AGRIF_ppm ,interp2 = Agrif_linear) 549 550 ! 3. Set location of interpolations 551 !---------------------------------- 552 CALL Agrif_Set_bc(tra_ice_id,(/0,ind1/)) 553 CALL Agrif_Set_bc(u_ice_id ,(/0,ind1/)) 554 CALL Agrif_Set_bc(v_ice_id ,(/0,ind1/)) 555 556 ! 4. Set update type in case 2 ways (child=>parent) (normal & tangent to the grid cell for velocities) 557 !-------------------------------------------------- 558 # if defined UPD_HIGH 559 CALL Agrif_Set_Updatetype(tra_ice_id, update = Agrif_Update_Full_Weighting) 560 CALL Agrif_Set_Updatetype(u_ice_id , update1 = Agrif_Update_Average , update2 = Agrif_Update_Full_Weighting) 561 CALL Agrif_Set_Updatetype(v_ice_id , update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average ) 562 #else 563 CALL Agrif_Set_Updatetype(tra_ice_id, update = AGRIF_Update_Average) 564 CALL Agrif_Set_Updatetype(u_ice_id , update1 = Agrif_Update_Copy , update2 = Agrif_Update_Average) 565 CALL Agrif_Set_Updatetype(v_ice_id , update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy ) 566 #endif 567 568 END SUBROUTINE agrif_declare_var_ice 569 #endif 570 571 572 # if defined key_top 573 SUBROUTINE Agrif_InitValues_cont_top 574 !!---------------------------------------------------------------------- 575 !! *** ROUTINE Agrif_InitValues_cont_top *** 576 !!---------------------------------------------------------------------- 577 USE Agrif_Util 578 USE oce 579 USE dom_oce 580 USE nemogcm 581 USE par_trc 582 USE lib_mpp 583 USE trc 584 USE in_out_manager 585 USE agrif_oce_sponge 586 USE agrif_top_update 587 USE agrif_top_interp 588 USE agrif_top_sponge 455 589 !! 456 !! ** Purpose :: Initialisation of variables to be interpolated for ice 457 !!---------------------------------------------------------------------- 458 USE Agrif_Util 459 USE sbc_oce, ONLY : nn_fsbc ! clem: necessary otherwise Agrif_Parent(nn_fsbc) = nn_fsbc 460 USE ice 461 USE agrif_ice 462 USE in_out_manager 463 USE agrif_ice_interp 464 USE lib_mpp 465 ! 466 IMPLICIT NONE 467 !!---------------------------------------------------------------------- 468 ! 469 ! Declaration of the type of variable which have to be interpolated (parent=>child) 470 !---------------------------------------------------------------------------------- 471 CALL agrif_declare_var_ice 472 473 ! Controls 474 475 ! clem: For some reason, nn_fsbc(child)/=1 does not work properly (signal can be largely degraded by the agrif zoom) 476 ! the run must satisfy CFL=Uice/(dx/dt) < 0.6/nn_fsbc(child) 477 ! therefore, if nn_fsbc(child)>1 one must reduce the time-step in proportion to nn_fsbc(child), which is not acceptable 478 ! If a solution is found, the following stop could be removed because the rest of the code take nn_fsbc(child) into account 479 IF( nn_fsbc > 1 ) CALL ctl_stop('nn_fsbc(child) must be set to 1 otherwise agrif and sea-ice may not work properly') 480 481 ! stop if rhot * nn_fsbc(parent) /= N * nn_fsbc(child) with N being integer 482 IF( MOD( Agrif_irhot() * Agrif_Parent(nn_fsbc), nn_fsbc ) /= 0 ) THEN 483 CALL ctl_stop('rhot * nn_fsbc(parent) /= N * nn_fsbc(child), therefore nn_fsbc(child) should be set to 1 or nn_fsbc(parent)') 484 ENDIF 485 ! First Interpolations (using "after" ice subtime step => nbstep_ice=1) 486 !---------------------------------------------------------------------- 487 nbstep_ice = ( Agrif_irhot() * Agrif_Parent(nn_fsbc) / nn_fsbc ) ! clem: to have calledweight=1 in interp (otherwise the western border of the zoom is wrong) 488 CALL agrif_interp_ice('U') ! interpolation of ice velocities 489 CALL agrif_interp_ice('V') ! interpolation of ice velocities 490 CALL agrif_interp_ice('T') ! interpolation of ice tracers 491 nbstep_ice = 0 492 493 ! 494 END SUBROUTINE Agrif_InitValues_cont_ice 495 496 SUBROUTINE agrif_declare_var_ice 497 !!---------------------------------------------------------------------- 498 !! *** ROUTINE agrif_declare_var_ice *** 499 !! 500 !! ** Purpose :: Declaration of variables to be interpolated for ice 501 !!---------------------------------------------------------------------- 502 USE Agrif_Util 503 USE ice 504 USE par_oce, ONLY : nbghostcells 505 ! 506 IMPLICIT NONE 507 ! 508 INTEGER :: ind1, ind2, ind3 509 !!---------------------------------------------------------------------- 510 ! 511 ! 1. Declaration of the type of variable which have to be interpolated (parent=>child) 512 ! agrif_declare_variable(position,1st point index,--,--,dimensions,name) 513 ! ex.: position=> 1,1 = not-centered (in i and j) 514 ! 2,2 = centered ( - ) 515 ! index => 1,1 = one ghost line 516 ! 2,2 = two ghost lines 517 !------------------------------------------------------------------------------------- 518 ind1 = nbghostcells 519 ind2 = 1 + nbghostcells 520 ind3 = 2 + nbghostcells 521 CALL agrif_declare_variable((/2,2,0/),(/ind3,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpl*(8+nlay_s+nlay_i)/),tra_ice_id) 522 CALL agrif_declare_variable((/1,2/) ,(/ind2,ind3/) ,(/'x','y'/) ,(/1,1/) ,(/nlci,nlcj/) ,u_ice_id ) 523 CALL agrif_declare_variable((/2,1/) ,(/ind3,ind2/) ,(/'x','y'/) ,(/1,1/) ,(/nlci,nlcj/) ,v_ice_id ) 524 525 ! 2. Set interpolations (normal & tangent to the grid cell for velocities) 526 !----------------------------------- 527 CALL Agrif_Set_bcinterp(tra_ice_id, interp = AGRIF_linear) 528 CALL Agrif_Set_bcinterp(u_ice_id , interp1 = Agrif_linear,interp2 = AGRIF_ppm ) 529 CALL Agrif_Set_bcinterp(v_ice_id , interp1 = AGRIF_ppm ,interp2 = Agrif_linear) 530 531 ! 3. Set location of interpolations 532 !---------------------------------- 533 CALL Agrif_Set_bc(tra_ice_id,(/0,ind1/)) 534 CALL Agrif_Set_bc(u_ice_id ,(/0,ind1/)) 535 CALL Agrif_Set_bc(v_ice_id ,(/0,ind1/)) 536 537 ! 4. Set update type in case 2 ways (child=>parent) (normal & tangent to the grid cell for velocities) 538 !-------------------------------------------------- 539 # if defined UPD_HIGH 540 CALL Agrif_Set_Updatetype(tra_ice_id, update = Agrif_Update_Full_Weighting) 541 CALL Agrif_Set_Updatetype(u_ice_id , update1 = Agrif_Update_Average , update2 = Agrif_Update_Full_Weighting) 542 CALL Agrif_Set_Updatetype(v_ice_id , update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average ) 543 #else 544 CALL Agrif_Set_Updatetype(tra_ice_id, update = AGRIF_Update_Average) 545 CALL Agrif_Set_Updatetype(u_ice_id , update1 = Agrif_Update_Copy , update2 = Agrif_Update_Average) 546 CALL Agrif_Set_Updatetype(v_ice_id , update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy ) 547 #endif 548 549 END SUBROUTINE agrif_declare_var_ice 550 #endif 551 552 553 # if defined key_top 554 SUBROUTINE Agrif_InitValues_cont_top 555 !!---------------------------------------------------------------------- 556 !! *** ROUTINE Agrif_InitValues_cont_top *** 557 !! 558 !! ** Purpose :: Declaration of variables to be interpolated 559 !!---------------------------------------------------------------------- 560 USE Agrif_Util 561 USE oce 562 USE dom_oce 563 USE nemogcm 564 USE par_trc 565 USE lib_mpp 566 USE trc 567 USE in_out_manager 568 USE agrif_oce_sponge 569 USE agrif_top_update 570 USE agrif_top_interp 571 USE agrif_top_sponge 572 !! 573 IMPLICIT NONE 574 ! 575 CHARACTER(len=10) :: cl_check1, cl_check2, cl_check3 576 LOGICAL :: check_namelist 577 !!---------------------------------------------------------------------- 578 579 580 ! 1. Declaration of the type of variable which have to be interpolated 581 !--------------------------------------------------------------------- 582 CALL agrif_declare_var_top 583 584 ! 2. First interpolations of potentially non zero fields 585 !------------------------------------------------------- 586 Agrif_SpecialValue=0. 587 Agrif_UseSpecialValue = .TRUE. 588 CALL Agrif_Bc_variable(trn_id,calledweight=1.,procname=interptrn) 589 Agrif_UseSpecialValue = .FALSE. 590 CALL Agrif_Sponge 591 tabspongedone_trn = .FALSE. 592 CALL Agrif_Bc_variable(trn_sponge_id,calledweight=1.,procname=interptrn_sponge) 593 ! reset tsa to zero 594 tra(:,:,:,:) = 0. 595 596 597 ! 3. Some controls 598 !----------------- 599 check_namelist = .TRUE. 600 601 IF( check_namelist ) THEN 602 ! Check time steps 590 IMPLICIT NONE 591 ! 592 CHARACTER(len=10) :: cl_check1, cl_check2, cl_check3 593 LOGICAL :: check_namelist 594 !!---------------------------------------------------------------------- 595 596 ! 1. Declaration of the type of variable which have to be interpolated 597 !--------------------------------------------------------------------- 598 CALL agrif_declare_var_top 599 600 ! 2. First interpolations of potentially non zero fields 601 !------------------------------------------------------- 602 Agrif_SpecialValue=0._wp 603 Agrif_UseSpecialValue = .TRUE. 604 CALL Agrif_Bc_variable(trn_id,calledweight=1.,procname=interptrn) 605 Agrif_UseSpecialValue = .FALSE. 606 CALL Agrif_Sponge 607 tabspongedone_trn = .FALSE. 608 CALL Agrif_Bc_variable(trn_sponge_id,calledweight=1.,procname=interptrn_sponge) 609 ! reset tsa to zero 610 tra(:,:,:,:) = 0._wp 611 612 ! 3. Some controls 613 !----------------- 614 check_namelist = .TRUE. 615 616 IF( check_namelist ) THEN 617 ! Check time steps 603 618 IF( NINT(Agrif_Rhot()) * NINT(rdt) .NE. Agrif_Parent(rdt) ) THEN 604 619 WRITE(cl_check1,*) Agrif_Parent(rdt) … … 630 645 ENDIF 631 646 ! 632 END SUBROUTINE Agrif_InitValues_cont_top633 634 635 SUBROUTINE agrif_declare_var_top647 END SUBROUTINE Agrif_InitValues_cont_top 648 649 650 SUBROUTINE agrif_declare_var_top 636 651 !!---------------------------------------------------------------------- 637 652 !! *** ROUTINE agrif_declare_var_top *** 653 !!---------------------------------------------------------------------- 654 USE agrif_util 655 USE agrif_oce 656 USE dom_oce 657 USE trc 638 658 !! 639 !! ** Purpose :: Declaration of TOP variables to be interpolated 640 !!---------------------------------------------------------------------- 641 USE agrif_util 642 USE agrif_oce 643 USE dom_oce 644 USE trc 645 !! 646 IMPLICIT NONE 647 ! 648 INTEGER :: ind1, ind2, ind3 649 !!---------------------------------------------------------------------- 650 651 ! 1. Declaration of the type of variable which have to be interpolated 652 !--------------------------------------------------------------------- 653 ind1 = nbghostcells 654 ind2 = 1 + nbghostcells 655 ind3 = 2 + nbghostcells 659 IMPLICIT NONE 660 ! 661 INTEGER :: ind1, ind2, ind3 662 !!---------------------------------------------------------------------- 663 664 ! 1. Declaration of the type of variable which have to be interpolated 665 !--------------------------------------------------------------------- 666 ind1 = nbghostcells 667 ind2 = 1 + nbghostcells 668 ind3 = 2 + nbghostcells 656 669 # if defined key_vertical 657 CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra+1/),trn_id)658 CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra+1/),trn_sponge_id)670 CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra+1/),trn_id) 671 CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra+1/),trn_sponge_id) 659 672 # else 660 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra/),trn_id)661 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra/),trn_sponge_id)662 # endif 663 664 ! 2. Type of interpolation665 !-------------------------666 CALL Agrif_Set_bcinterp(trn_id,interp=AGRIF_linear)667 CALL Agrif_Set_bcinterp(trn_sponge_id,interp=AGRIF_linear)668 669 ! 3. Location of interpolation670 !-----------------------------671 CALL Agrif_Set_bc(trn_id,(/0,ind1/))672 CALL Agrif_Set_bc(trn_sponge_id,(/-nn_sponge_len*Agrif_irhox()-1,0/))673 674 ! 4. Update type675 !---------------673 CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra/),trn_id) 674 CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra/),trn_sponge_id) 675 # endif 676 677 ! 2. Type of interpolation 678 !------------------------- 679 CALL Agrif_Set_bcinterp(trn_id,interp=AGRIF_linear) 680 CALL Agrif_Set_bcinterp(trn_sponge_id,interp=AGRIF_linear) 681 682 ! 3. Location of interpolation 683 !----------------------------- 684 CALL Agrif_Set_bc(trn_id,(/0,ind1-1/)) 685 CALL Agrif_Set_bc(trn_sponge_id,(/-nn_sponge_len*Agrif_irhox()-1,0/)) 686 687 ! 4. Update type 688 !--------------- 676 689 # if defined UPD_HIGH 677 CALL Agrif_Set_Updatetype(trn_id, update = Agrif_Update_Full_Weighting)690 CALL Agrif_Set_Updatetype(trn_id, update = Agrif_Update_Full_Weighting) 678 691 #else 679 CALL Agrif_Set_Updatetype(trn_id, update = AGRIF_Update_Average)692 CALL Agrif_Set_Updatetype(trn_id, update = AGRIF_Update_Average) 680 693 #endif 681 694 ! 682 END SUBROUTINE agrif_declare_var_top683 # endif 684 685 SUBROUTINE Agrif_detect( kg, ksizex )695 END SUBROUTINE agrif_declare_var_top 696 # endif 697 698 SUBROUTINE Agrif_detect( kg, ksizex ) 686 699 !!---------------------------------------------------------------------- 687 700 !! *** ROUTINE Agrif_detect *** 688 701 !!---------------------------------------------------------------------- 689 INTEGER, DIMENSION(2) :: ksizex 690 INTEGER, DIMENSION(ksizex(1),ksizex(2)) :: kg 691 !!---------------------------------------------------------------------- 692 ! 693 RETURN 694 ! 695 END SUBROUTINE Agrif_detect 696 697 698 SUBROUTINE agrif_nemo_init 702 INTEGER, DIMENSION(2) :: ksizex 703 INTEGER, DIMENSION(ksizex(1),ksizex(2)) :: kg 704 !!---------------------------------------------------------------------- 705 ! 706 RETURN 707 ! 708 END SUBROUTINE Agrif_detect 709 710 SUBROUTINE agrif_nemo_init 699 711 !!---------------------------------------------------------------------- 700 712 !! *** ROUTINE agrif_init *** 701 713 !!---------------------------------------------------------------------- 702 USE agrif_oce703 USE agrif_ice704 USE in_out_manager705 USE lib_mpp706 !!707 IMPLICIT NONE708 !709 INTEGER :: ios ! Local integer output status for namelist read710 INTEGER :: iminspon711 NAMELIST/namagrif/ rn_sponge_tra, rn_sponge_dyn,ln_spc_dyn, ln_chk_bathy714 USE agrif_oce 715 USE agrif_ice 716 USE in_out_manager 717 USE lib_mpp 718 !! 719 IMPLICIT NONE 720 ! 721 INTEGER :: ios ! Local integer output status for namelist read 722 NAMELIST/namagrif/ ln_agrif_2way, rn_sponge_tra, rn_sponge_dyn, rn_trelax_tra, rn_trelax_dyn, & 723 & ln_spc_dyn, ln_chk_bathy 712 724 !!-------------------------------------------------------------------------------------- 713 !714 REWIND( numnam_ref ) ! Namelist namagrif in reference namelist : AGRIF zoom715 READ ( numnam_ref, namagrif, IOSTAT = ios, ERR = 901)725 ! 726 REWIND( numnam_ref ) ! Namelist namagrif in reference namelist : AGRIF zoom 727 READ ( numnam_ref, namagrif, IOSTAT = ios, ERR = 901) 716 728 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namagrif in reference namelist' ) 717 REWIND( numnam_cfg ) ! Namelist namagrif in configuration namelist : AGRIF zoom718 READ ( numnam_cfg, namagrif, IOSTAT = ios, ERR = 902 )729 REWIND( numnam_cfg ) ! Namelist namagrif in configuration namelist : AGRIF zoom 730 READ ( numnam_cfg, namagrif, IOSTAT = ios, ERR = 902 ) 719 731 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namagrif in configuration namelist' ) 720 IF(lwm) WRITE ( numond, namagrif ) 721 ! 722 IF(lwp) THEN ! control print 723 WRITE(numout,*) 724 WRITE(numout,*) 'agrif_nemo_init : AGRIF parameters' 725 WRITE(numout,*) '~~~~~~~~~~~~~~~' 726 WRITE(numout,*) ' Namelist namagrif : set AGRIF parameters' 727 WRITE(numout,*) ' sponge coefficient for tracers rn_sponge_tra = ', rn_sponge_tra, ' s' 728 WRITE(numout,*) ' sponge coefficient for dynamics rn_sponge_tra = ', rn_sponge_dyn, ' s' 729 WRITE(numout,*) ' use special values for dynamics ln_spc_dyn = ', ln_spc_dyn 730 WRITE(numout,*) ' check bathymetry ln_chk_bathy = ', ln_chk_bathy 731 ENDIF 732 ! 733 ! convert DOCTOR namelist name into OLD names 734 visc_tra = rn_sponge_tra 735 visc_dyn = rn_sponge_dyn 736 ! 737 ! Check sponge length: 738 iminspon = MIN(FLOOR(REAL(jpiglo-4)/REAL(2*Agrif_irhox())), FLOOR(REAL(jpjglo-4)/REAL(2*Agrif_irhox())) ) 739 IF (lk_mpp) iminspon = MIN(iminspon,FLOOR(REAL(jpi-2)/REAL(Agrif_irhox())), FLOOR(REAL(jpj-2)/REAL(Agrif_irhox())) ) 740 IF (nn_sponge_len > iminspon) CALL ctl_stop('agrif sponge length is too large') 741 ! 742 IF( agrif_oce_alloc() > 0 ) CALL ctl_warn('agrif agrif_oce_alloc: allocation of arrays failed') 743 ! 744 END SUBROUTINE agrif_nemo_init 732 IF(lwm) WRITE ( numond, namagrif ) 733 ! 734 IF(lwp) THEN ! control print 735 WRITE(numout,*) 736 WRITE(numout,*) 'agrif_nemo_init : AGRIF parameters' 737 WRITE(numout,*) '~~~~~~~~~~~~~~~' 738 WRITE(numout,*) ' Namelist namagrif : set AGRIF parameters' 739 WRITE(numout,*) ' Two way nesting activated ln_agrif_2way = ', ln_agrif_2way 740 WRITE(numout,*) ' sponge coefficient for tracers rn_sponge_tra = ', rn_sponge_tra, ' m^2/s' 741 WRITE(numout,*) ' sponge coefficient for dynamics rn_sponge_tra = ', rn_sponge_dyn, ' m^2/s' 742 WRITE(numout,*) ' time relaxation for tracers rn_trelax_tra = ', rn_trelax_tra, ' ad.' 743 WRITE(numout,*) ' time relaxation for dynamics rn_trelax_dyn = ', rn_trelax_dyn, ' ad.' 744 WRITE(numout,*) ' use special values for dynamics ln_spc_dyn = ', ln_spc_dyn 745 WRITE(numout,*) ' check bathymetry ln_chk_bathy = ', ln_chk_bathy 746 ENDIF 747 ! 748 ! 749 IF( agrif_oce_alloc() > 0 ) CALL ctl_warn('agrif agrif_oce_alloc: allocation of arrays failed') 750 ! 751 END SUBROUTINE agrif_nemo_init 745 752 746 753 # if defined key_mpp_mpi 747 754 748 SUBROUTINE Agrif_InvLoc( indloc, nprocloc, i, indglob )755 SUBROUTINE Agrif_InvLoc( indloc, nprocloc, i, indglob ) 749 756 !!---------------------------------------------------------------------- 750 757 !! *** ROUTINE Agrif_InvLoc *** 751 758 !!---------------------------------------------------------------------- 752 USE dom_oce 753 !! 754 IMPLICIT NONE 755 ! 756 INTEGER :: indglob, indloc, nprocloc, i 757 !!---------------------------------------------------------------------- 758 ! 759 SELECT CASE( i ) 760 CASE(1) ; indglob = indloc + nimppt(nprocloc+1) - 1 761 CASE(2) ; indglob = indloc + njmppt(nprocloc+1) - 1 762 CASE DEFAULT 763 indglob = indloc 764 END SELECT 765 ! 766 END SUBROUTINE Agrif_InvLoc 767 768 769 SUBROUTINE Agrif_get_proc_info( imin, imax, jmin, jmax ) 759 USE dom_oce 760 !! 761 IMPLICIT NONE 762 ! 763 INTEGER :: indglob, indloc, nprocloc, i 764 !!---------------------------------------------------------------------- 765 ! 766 SELECT CASE( i ) 767 CASE(1) ; indglob = indloc + nimppt(nprocloc+1) - 1 768 CASE(2) ; indglob = indloc + njmppt(nprocloc+1) - 1 769 CASE DEFAULT 770 indglob = indloc 771 END SELECT 772 ! 773 END SUBROUTINE Agrif_InvLoc 774 775 SUBROUTINE Agrif_get_proc_info( imin, imax, jmin, jmax ) 770 776 !!---------------------------------------------------------------------- 771 777 !! *** ROUTINE Agrif_get_proc_info *** 772 778 !!---------------------------------------------------------------------- 773 USE par_oce 774 !! 775 IMPLICIT NONE 776 ! 777 INTEGER, INTENT(out) :: imin, imax 778 INTEGER, INTENT(out) :: jmin, jmax 779 !!---------------------------------------------------------------------- 780 ! 781 imin = nimppt(Agrif_Procrank+1) ! ????? 782 jmin = njmppt(Agrif_Procrank+1) ! ????? 783 imax = imin + jpi - 1 784 jmax = jmin + jpj - 1 785 ! 786 END SUBROUTINE Agrif_get_proc_info 787 788 789 SUBROUTINE Agrif_estimate_parallel_cost(imin, imax,jmin, jmax, nbprocs, grid_cost) 779 USE par_oce 780 !! 781 IMPLICIT NONE 782 ! 783 INTEGER, INTENT(out) :: imin, imax 784 INTEGER, INTENT(out) :: jmin, jmax 785 !!---------------------------------------------------------------------- 786 ! 787 imin = nimppt(Agrif_Procrank+1) ! ????? 788 jmin = njmppt(Agrif_Procrank+1) ! ????? 789 imax = imin + jpi - 1 790 jmax = jmin + jpj - 1 791 ! 792 END SUBROUTINE Agrif_get_proc_info 793 794 SUBROUTINE Agrif_estimate_parallel_cost(imin, imax,jmin, jmax, nbprocs, grid_cost) 790 795 !!---------------------------------------------------------------------- 791 796 !! *** ROUTINE Agrif_estimate_parallel_cost *** 792 797 !!---------------------------------------------------------------------- 793 USE par_oce794 !!795 IMPLICIT NONE796 !797 INTEGER, INTENT(in) :: imin, imax798 INTEGER, INTENT(in) :: jmin, jmax799 INTEGER, INTENT(in) :: nbprocs800 REAL(wp), INTENT(out) :: grid_cost801 !!---------------------------------------------------------------------- 802 !803 grid_cost = REAL(imax-imin+1,wp)*REAL(jmax-jmin+1,wp) / REAL(nbprocs,wp)804 !805 END SUBROUTINE Agrif_estimate_parallel_cost798 USE par_oce 799 !! 800 IMPLICIT NONE 801 ! 802 INTEGER, INTENT(in) :: imin, imax 803 INTEGER, INTENT(in) :: jmin, jmax 804 INTEGER, INTENT(in) :: nbprocs 805 REAL(wp), INTENT(out) :: grid_cost 806 !!---------------------------------------------------------------------- 807 ! 808 grid_cost = REAL(imax-imin+1,wp)*REAL(jmax-jmin+1,wp) / REAL(nbprocs,wp) 809 ! 810 END SUBROUTINE Agrif_estimate_parallel_cost 806 811 807 812 # endif 808 813 809 814 #else 810 SUBROUTINE Subcalledbyagrif815 SUBROUTINE Subcalledbyagrif 811 816 !!---------------------------------------------------------------------- 812 817 !! *** ROUTINE Subcalledbyagrif *** 813 818 !!---------------------------------------------------------------------- 814 WRITE(*,*) 'Impossible to be here'815 END SUBROUTINE Subcalledbyagrif819 WRITE(*,*) 'Impossible to be here' 820 END SUBROUTINE Subcalledbyagrif 816 821 #endif -
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DYN/dynspg_ts.F90
r12184 r12191 505 505 ! ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d 506 506 IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e ) 507 ! 507 ! 508 508 ! ! resulting flux at mid-step (not over the full domain) 509 509 zhU(1:jpim1,1:jpj ) = e2u(1:jpim1,1:jpj ) * ua_e(1:jpim1,1:jpj ) * zhup2_e(1:jpim1,1:jpj ) ! not jpi-column … … 512 512 #if defined key_agrif 513 513 ! Set fluxes during predictor step to ensure volume conservation 514 IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN 515 IF((nbondi == -1).OR.(nbondi == 2)) THEN 516 DO jj = 1, jpj 517 zhU(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj) 518 zhV(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj) 519 END DO 520 ENDIF 521 IF((nbondi == 1).OR.(nbondi == 2)) THEN 522 DO jj=1,jpj 523 zhU(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj) 524 zhV(nlci-nbghostcells :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells :nlci-1,jj) 525 END DO 526 ENDIF 527 IF((nbondj == -1).OR.(nbondj == 2)) THEN 528 DO ji=1,jpi 529 zhV(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1) 530 zhU(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1) 531 END DO 532 ENDIF 533 IF((nbondj == 1).OR.(nbondj == 2)) THEN 534 DO ji=1,jpi 535 zhV(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2) 536 zhU(ji,nlcj-nbghostcells :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells :nlcj-1) 537 END DO 538 ENDIF 539 ENDIF 514 IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) CALL agrif_dyn_ts_flux( jn, zhU, zhV ) 540 515 #endif 541 516 IF( ln_wd_il ) CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rdtbt) !!gm wad_lmt_bt use of lbc_lnk on zhU, zhV
Note: See TracChangeset
for help on using the changeset viewer.