- Timestamp:
- 2019-12-11T16:56:06+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
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 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.