- Timestamp:
- 2019-12-12T17:41:04+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11943_MERGE_2019/src/NST/agrif_oce.F90
r11949 r12229 17 17 18 18 PUBLIC agrif_oce_alloc ! routine called by nemo_init in nemogcm.F90 19 #if defined key_vertical 20 PUBLIC reconstructandremap ! remapping routine 21 #endif 19 22 20 ! !!* Namelist namagrif: AGRIF parameters 23 LOGICAL , PUBLIC :: ln_spc_dyn = .FALSE. !: 24 INTEGER , PUBLIC, PARAMETER :: nn_sponge_len = 2 !: Sponge width (in number of parent grid points) 21 LOGICAL , PUBLIC :: ln_agrif_2way = .TRUE. !: activate two way nesting 22 LOGICAL , PUBLIC :: ln_spc_dyn = .FALSE. !: use zeros (.false.) or not (.true.) in 23 !: bdys dynamical fields interpolation 25 24 REAL(wp), PUBLIC :: rn_sponge_tra = 2800. !: sponge coeff. for tracers 26 25 REAL(wp), PUBLIC :: rn_sponge_dyn = 2800. !: sponge coeff. for dynamics 26 REAL(wp), PUBLIC :: rn_trelax_tra = 0.01 !: time relaxation parameter for tracers 27 REAL(wp), PUBLIC :: rn_trelax_dyn = 0.01 !: time relaxation parameter for momentum 27 28 LOGICAL , PUBLIC :: ln_chk_bathy = .FALSE. !: check of parent bathymetry 28 LOGICAL , PUBLIC :: lk_agrif_clp = .FALSE. !: Force clamped bcs 29 ! !!! OLD namelist names 30 REAL(wp), PUBLIC :: visc_tra !: sponge coeff. for tracers 31 REAL(wp), PUBLIC :: visc_dyn !: sponge coeff. for dynamics 32 29 ! 30 INTEGER , PUBLIC, PARAMETER :: nn_sponge_len = 2 !: Sponge width (in number of parent grid points) 33 31 LOGICAL , PUBLIC :: spongedoneT = .FALSE. !: tracer sponge layer indicator 34 32 LOGICAL , PUBLIC :: spongedoneU = .FALSE. !: dynamics sponge layer indicator … … 42 40 LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_u 43 41 LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_v 44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fsaht_spu, fsaht_spv !: sponge diffusivities 45 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fsahm_spt, fsahm_spf !: sponge viscosities 42 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: utint_stage 43 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vtint_stage 44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fspu, fspv !: sponge arrays 45 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fspt, fspf !: " " 46 46 47 47 ! Barotropic arrays used to store open boundary data during time-splitting loop: 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ubdy_w, vbdy_w, hbdy_w 49 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ubdy_e, vbdy_e, hbdy_e 50 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ubdy_n, vbdy_n, hbdy_n 51 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ubdy_s, vbdy_s, hbdy_s 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ubdy, vbdy, hbdy 52 49 INTEGER , PUBLIC, SAVE :: Kbb_a, Kmm_a, Krhs_a !: AGRIF module-specific copies of time-level indices 53 50 51 # if defined key_vertical 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht0_parent, hu0_parent, hv0_parent 53 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt_parent, mbku_parent, mbkv_parent 54 # endif 54 55 55 56 INTEGER, PUBLIC :: tsn_id ! AGRIF profile for tracers interpolation and update … … 65 66 INTEGER, PUBLIC :: scales_t_id 66 67 INTEGER, PUBLIC :: avt_id, avm_id, en_id ! TKE related identificators 67 INTEGER, PUBLIC :: umsk_id, vmsk_id68 INTEGER, PUBLIC :: mbkt_id, ht0_id 68 69 INTEGER, PUBLIC :: kindic_agr 69 70 … … 83 84 ierr(:) = 0 84 85 ! 85 ALLOCATE( fsaht_spu(jpi,jpj), fsaht_spv(jpi,jpj), & 86 & fsahm_spt(jpi,jpj), fsahm_spf(jpi,jpj), & 87 & tabspongedone_tsn(jpi,jpj), & 86 ALLOCATE( fspu(jpi,jpj), fspv(jpi,jpj), & 87 & fspt(jpi,jpj), fspf(jpi,jpj), & 88 & tabspongedone_tsn(jpi,jpj), & 89 & utint_stage(jpi,jpj), vtint_stage(jpi,jpj), & 88 90 # if defined key_top 89 91 & tabspongedone_trn(jpi,jpj), & 90 # endif 92 # endif 93 # if defined key_vertical 94 & ht0_parent(jpi,jpj), mbkt_parent(jpi,jpj), & 95 & hu0_parent(jpi,jpj), mbku_parent(jpi,jpj), & 96 & hv0_parent(jpi,jpj), mbkv_parent(jpi,jpj), & 97 # endif 91 98 & tabspongedone_u (jpi,jpj), & 92 99 & tabspongedone_v (jpi,jpj), STAT = ierr(1) ) 93 100 94 ALLOCATE( ubdy_w(nbghostcells,jpj), vbdy_w(nbghostcells,jpj), hbdy_w(nbghostcells,jpj), & 95 & ubdy_e(nbghostcells,jpj), vbdy_e(nbghostcells,jpj), hbdy_e(nbghostcells,jpj), & 96 & ubdy_n(jpi,nbghostcells), vbdy_n(jpi,nbghostcells), hbdy_n(jpi,nbghostcells), & 97 & ubdy_s(jpi,nbghostcells), vbdy_s(jpi,nbghostcells), hbdy_s(jpi,nbghostcells), STAT = ierr(2) ) 101 ALLOCATE( ubdy(jpi,jpj), vbdy(jpi,jpj), hbdy(jpi,jpj), STAT = ierr(2) ) 98 102 99 103 agrif_oce_alloc = MAXVAL(ierr) … … 101 105 END FUNCTION agrif_oce_alloc 102 106 103 #if defined key_vertical104 SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout)105 !!----------------------------------------------------------------------106 !! *** FUNCTION reconstructandremap ***107 !!----------------------------------------------------------------------108 IMPLICIT NONE109 INTEGER N, Nout110 REAL(wp) tabin(N), tabout(Nout)111 REAL(wp) hin(N), hout(Nout)112 REAL(wp) coeffremap(N,3),zwork(N,3)113 REAL(wp) zwork2(N+1,3)114 INTEGER jk115 DOUBLE PRECISION, PARAMETER :: dsmll=1.0d-8116 REAL(wp) q,q01,q02,q001,q002,q0117 REAL(wp) z_win(1:N+1), z_wout(1:Nout+1)118 REAL(wp),PARAMETER :: dpthin = 1.D-3119 INTEGER :: k1, kbox, ktop, ka, kbot120 REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop121 122 z_win(1)=0.; z_wout(1)= 0.123 DO jk=1,N124 z_win(jk+1)=z_win(jk)+hin(jk)125 ENDDO126 127 DO jk=1,Nout128 z_wout(jk+1)=z_wout(jk)+hout(jk)129 ENDDO130 131 DO jk=2,N132 zwork(jk,1)=1./(hin(jk-1)+hin(jk))133 ENDDO134 135 DO jk=2,N-1136 q0 = 1./(hin(jk-1)+hin(jk)+hin(jk+1))137 zwork(jk,2)=hin(jk-1)+2.*hin(jk)+hin(jk+1)138 zwork(jk,3)=q0139 ENDDO140 141 DO jk= 2,N142 zwork2(jk,1)=zwork(jk,1)*(tabin(jk)-tabin(jk-1))143 ENDDO144 145 coeffremap(:,1) = tabin(:)146 147 DO jk=2,N-1148 q001 = hin(jk)*zwork2(jk+1,1)149 q002 = hin(jk)*zwork2(jk,1)150 IF (q001*q002 < 0) then151 q001 = 0.152 q002 = 0.153 ENDIF154 q=zwork(jk,2)155 q01=q*zwork2(jk+1,1)156 q02=q*zwork2(jk,1)157 IF (abs(q001) > abs(q02)) q001 = q02158 IF (abs(q002) > abs(q01)) q002 = q01159 160 q=(q001-q002)*zwork(jk,3)161 q001=q001-q*hin(jk+1)162 q002=q002+q*hin(jk-1)163 164 coeffremap(jk,3)=coeffremap(jk,1)+q001165 coeffremap(jk,2)=coeffremap(jk,1)-q002166 167 zwork2(jk,1)=(2.*q001-q002)**2168 zwork2(jk,2)=(2.*q002-q001)**2169 ENDDO170 171 DO jk=1,N172 IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN173 coeffremap(jk,3) = coeffremap(jk,1)174 coeffremap(jk,2) = coeffremap(jk,1)175 zwork2(jk,1) = 0.176 zwork2(jk,2) = 0.177 ENDIF178 ENDDO179 180 DO jk=2,N181 q002=max(zwork2(jk-1,2),dsmll)182 q001=max(zwork2(jk,1),dsmll)183 zwork2(jk,3)=(q001*coeffremap(jk-1,3)+q002*coeffremap(jk,2))/(q001+q002)184 ENDDO185 186 zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3)187 zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3)188 189 DO jk=1,N190 q01=zwork2(jk+1,3)-coeffremap(jk,1)191 q02=coeffremap(jk,1)-zwork2(jk,3)192 q001=2.*q01193 q002=2.*q02194 IF (q01*q02<0) then195 q01=0.196 q02=0.197 ELSEIF (abs(q01)>abs(q002)) then198 q01=q002199 ELSEIF (abs(q02)>abs(q001)) then200 q02=q001201 ENDIF202 coeffremap(jk,2)=coeffremap(jk,1)-q02203 coeffremap(jk,3)=coeffremap(jk,1)+q01204 ENDDO205 206 zbot=0.0207 kbot=1208 DO jk=1,Nout209 ztop=zbot !top is bottom of previous layer210 ktop=kbot211 IF (ztop.GE.z_win(ktop+1)) then212 ktop=ktop+1213 ENDIF214 215 zbot=z_wout(jk+1)216 zthk=zbot-ztop217 218 IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(Nout+1)) THEN219 220 kbot=ktop221 DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.N)222 kbot=kbot+1223 ENDDO224 zbox=zbot225 DO k1= jk+1,Nout226 IF (z_wout(k1+1)-z_wout(k1).GT.dpthin) THEN227 exit !thick layer228 ELSE229 zbox=z_wout(k1+1) !include thin adjacent layers230 IF(zbox.EQ.z_wout(Nout+1)) THEN231 exit !at bottom232 ENDIF233 ENDIF234 ENDDO235 zthk=zbox-ztop236 237 kbox=ktop238 DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.N)239 kbox=kbox+1240 ENDDO241 242 IF(ktop.EQ.kbox) THEN243 IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN244 IF(hin(kbox).GT.dpthin) THEN245 q001 = (zbox-z_win(kbox))/hin(kbox)246 q002 = (ztop-z_win(kbox))/hin(kbox)247 q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002)248 q02=q01-1.+(q001+q002)249 q0=1.-q01-q02250 ELSE251 q0 = 1.0252 q01 = 0.253 q02 = 0.254 ENDIF255 tabout(jk)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3)256 ELSE257 tabout(jk) = tabin(kbox)258 ENDIF259 ELSE260 IF(ktop.LE.jk .AND. kbox.GE.jk) THEN261 ka = jk262 ELSEIF (kbox-ktop.GE.3) THEN263 ka = (kbox+ktop)/2264 ELSEIF (hin(ktop).GE.hin(kbox)) THEN265 ka = ktop266 ELSE267 ka = kbox268 ENDIF !choose ka269 270 offset=coeffremap(ka,1)271 272 qtop = z_win(ktop+1)-ztop !partial layer thickness273 IF(hin(ktop).GT.dpthin) THEN274 q=(ztop-z_win(ktop))/hin(ktop)275 q01=q*(q-1.)276 q02=q01+q277 q0=1-q01-q02278 ELSE279 q0 = 1.280 q01 = 0.281 q02 = 0.282 ENDIF283 284 tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop285 286 DO k1= ktop+1,kbox-1287 tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1)288 ENDDO !k1289 290 qbot = zbox-z_win(kbox) !partial layer thickness291 IF(hin(kbox).GT.dpthin) THEN292 q=qbot/hin(kbox)293 q01=(q-1.)**2294 q02=q01-1.+q295 q0=1-q01-q02296 ELSE297 q0 = 1.0298 q01 = 0.299 q02 = 0.300 ENDIF301 302 tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot303 304 rpsum=1.0d0/zthk305 tabout(jk)=offset+tsum*rpsum306 307 ENDIF !single or multiple layers308 ELSE309 IF (jk==1) THEN310 write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1)311 ENDIF312 tabout(jk) = tabout(jk-1)313 314 ENDIF !normal:thin layer315 ENDDO !jk316 317 return318 end subroutine reconstructandremap319 #endif320 321 107 #endif 322 108 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.