New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

agrif_oce.F90 in NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/src/NST – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/src/NST/agrif_oce.F90 @ 11205

Last change on this file since 11205 was 11205, checked in by jchanut, 5 years ago
#2199 1) Make sponge independent of sub-domain size. Partially masked open boundary segments are not taken into account anymore. To do so, sponge coefficients should be read in a file for realistic applications (then nesting tools need to be modified accordingly). 2) Replace East-West-North-South barotropic data arrays by a global 2d array. Then apply barotropic open boundary conditions thanks to mi0/mi1, mj0/mj1 indexes. 3) Call AGRIF bdy update one more time in dynspg_ts during extrapolation phase. This removes a dozen lines of code in dynspg_ts routine.
Property svn:keywords set to `Id`
File size: 11.6 KB

Line
1	MODULE agrif_oce
2	!!======================================================================
3	!! * MODULE agrif_oce *
4	!! AGRIF : define in memory AGRIF variables
5	!!----------------------------------------------------------------------
6	!! History : 2.0 ! 2007-12 (R. Benshila) Original code
7	!!----------------------------------------------------------------------
8	#if defined key_agrif
9	!!----------------------------------------------------------------------
10	!! 'key_agrif' AGRIF zoom
11	!!----------------------------------------------------------------------
12	USE par_oce ! ocean parameters
13	USE dom_oce ! domain parameters
14
15	IMPLICIT NONE
16	PRIVATE
17
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
22	! !!* 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)
25	REAL(wp), PUBLIC :: rn_sponge_tra = 2800. !: sponge coeff. for tracers
26	REAL(wp), PUBLIC :: rn_sponge_dyn = 2800. !: sponge coeff. for dynamics
27	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
33	LOGICAL , PUBLIC :: spongedoneT = .FALSE. !: tracer sponge layer indicator
34	LOGICAL , PUBLIC :: spongedoneU = .FALSE. !: dynamics sponge layer indicator
35	LOGICAL , PUBLIC :: lk_agrif_fstep = .TRUE. !: if true: first step
36	LOGICAL , PUBLIC :: lk_agrif_debug = .FALSE. !: if true: print debugging info
37
38	LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_tsn
39	# if defined key_top
40	LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_trn
41	# endif
42	LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_u
43	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
46
47	! Barotropic arrays used to store open boundary data during time-splitting loop:
48	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ubdy, vbdy, hbdy
49
50
51	INTEGER, PUBLIC :: tsn_id ! AGRIF profile for tracers interpolation and update
52	INTEGER, PUBLIC :: un_interp_id, vn_interp_id ! AGRIF profiles for interpolations
53	INTEGER, PUBLIC :: un_update_id, vn_update_id ! AGRIF profiles for udpates
54	INTEGER, PUBLIC :: tsn_sponge_id, un_sponge_id, vn_sponge_id ! AGRIF profiles for sponge layers
55	# if defined key_top
56	INTEGER, PUBLIC :: trn_id, trn_sponge_id
57	# endif
58	INTEGER, PUBLIC :: unb_id, vnb_id, ub2b_interp_id, vb2b_interp_id
59	INTEGER, PUBLIC :: ub2b_update_id, vb2b_update_id
60	INTEGER, PUBLIC :: e3t_id, e1u_id, e2v_id, sshn_id
61	INTEGER, PUBLIC :: scales_t_id
62	INTEGER, PUBLIC :: avt_id, avm_id, en_id ! TKE related identificators
63	INTEGER, PUBLIC :: umsk_id, vmsk_id
64	INTEGER, PUBLIC :: kindic_agr
65
66	!!----------------------------------------------------------------------
67	!! NEMO/NST 4.0 , NEMO Consortium (2018)
68	!! $Id$
69	!! Software governed by the CeCILL license (see ./LICENSE)
70	!!----------------------------------------------------------------------
71	CONTAINS
72
73	INTEGER FUNCTION agrif_oce_alloc()
74	!!----------------------------------------------------------------------
75	!! * FUNCTION agrif_oce_alloc *
76	!!----------------------------------------------------------------------
77	INTEGER, DIMENSION(2) :: ierr
78	!!----------------------------------------------------------------------
79	ierr(:) = 0
80	!
81	ALLOCATE( fsaht_spu(jpi,jpj), fsaht_spv(jpi,jpj), &
82	& fsahm_spt(jpi,jpj), fsahm_spf(jpi,jpj), &
83	& tabspongedone_tsn(jpi,jpj), &
84	# if defined key_top
85	& tabspongedone_trn(jpi,jpj), &
86	# endif
87	& tabspongedone_u (jpi,jpj), &
88	& tabspongedone_v (jpi,jpj), STAT = ierr(1) )
89
90	ALLOCATE( ubdy(jpi,jpj), vbdy(jpi,jpj), hbdy(jpi,jpj), STAT = ierr(2) )
91
92	agrif_oce_alloc = MAXVAL(ierr)
93	!
94	END FUNCTION agrif_oce_alloc
95
96	#if defined key_vertical
97	SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout)
98	!!----------------------------------------------------------------------
99	!! * FUNCTION reconstructandremap *
100	!!----------------------------------------------------------------------
101	IMPLICIT NONE
102	INTEGER N, Nout
103	REAL(wp) tabin(N), tabout(Nout)
104	REAL(wp) hin(N), hout(Nout)
105	REAL(wp) coeffremap(N,3),zwork(N,3)
106	REAL(wp) zwork2(N+1,3)
107	INTEGER jk
108	DOUBLE PRECISION, PARAMETER :: dsmll=1.0d-8
109	REAL(wp) q,q01,q02,q001,q002,q0
110	REAL(wp) z_win(1:N+1), z_wout(1:Nout+1)
111	REAL(wp),PARAMETER :: dpthin = 1.D-3
112	INTEGER :: k1, kbox, ktop, ka, kbot
113	REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop
114
115	z_win(1)=0.; z_wout(1)= 0.
116	DO jk=1,N
117	z_win(jk+1)=z_win(jk)+hin(jk)
118	ENDDO
119
120	DO jk=1,Nout
121	z_wout(jk+1)=z_wout(jk)+hout(jk)
122	ENDDO
123
124	DO jk=2,N
125	zwork(jk,1)=1./(hin(jk-1)+hin(jk))
126	ENDDO
127
128	DO jk=2,N-1
129	q0 = 1./(hin(jk-1)+hin(jk)+hin(jk+1))
130	zwork(jk,2)=hin(jk-1)+2.*hin(jk)+hin(jk+1)
131	zwork(jk,3)=q0
132	ENDDO
133
134	DO jk= 2,N
135	zwork2(jk,1)=zwork(jk,1)*(tabin(jk)-tabin(jk-1))
136	ENDDO
137
138	coeffremap(:,1) = tabin(:)
139
140	DO jk=2,N-1
141	q001 = hin(jk)*zwork2(jk+1,1)
142	q002 = hin(jk)*zwork2(jk,1)
143	IF (q001*q002 < 0) then
144	q001 = 0.
145	q002 = 0.
146	ENDIF
147	q=zwork(jk,2)
148	q01=q*zwork2(jk+1,1)
149	q02=q*zwork2(jk,1)
150	IF (abs(q001) > abs(q02)) q001 = q02
151	IF (abs(q002) > abs(q01)) q002 = q01
152
153	q=(q001-q002)*zwork(jk,3)
154	q001=q001-q*hin(jk+1)
155	q002=q002+q*hin(jk-1)
156
157	coeffremap(jk,3)=coeffremap(jk,1)+q001
158	coeffremap(jk,2)=coeffremap(jk,1)-q002
159
160	zwork2(jk,1)=(2.q001-q002)*2
161	zwork2(jk,2)=(2.q002-q001)*2
162	ENDDO
163
164	DO jk=1,N
165	IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN
166	coeffremap(jk,3) = coeffremap(jk,1)
167	coeffremap(jk,2) = coeffremap(jk,1)
168	zwork2(jk,1) = 0.
169	zwork2(jk,2) = 0.
170	ENDIF
171	ENDDO
172
173	DO jk=2,N
174	q002=max(zwork2(jk-1,2),dsmll)
175	q001=max(zwork2(jk,1),dsmll)
176	zwork2(jk,3)=(q001coeffremap(jk-1,3)+q002coeffremap(jk,2))/(q001+q002)
177	ENDDO
178
179	zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3)
180	zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3)
181
182	DO jk=1,N
183	q01=zwork2(jk+1,3)-coeffremap(jk,1)
184	q02=coeffremap(jk,1)-zwork2(jk,3)
185	q001=2.*q01
186	q002=2.*q02
187	IF (q01*q02<0) then
188	q01=0.
189	q02=0.
190	ELSEIF (abs(q01)>abs(q002)) then
191	q01=q002
192	ELSEIF (abs(q02)>abs(q001)) then
193	q02=q001
194	ENDIF
195	coeffremap(jk,2)=coeffremap(jk,1)-q02
196	coeffremap(jk,3)=coeffremap(jk,1)+q01
197	ENDDO
198
199	zbot=0.0
200	kbot=1
201	DO jk=1,Nout
202	ztop=zbot !top is bottom of previous layer
203	ktop=kbot
204	IF (ztop.GE.z_win(ktop+1)) then
205	ktop=ktop+1
206	ENDIF
207
208	zbot=z_wout(jk+1)
209	zthk=zbot-ztop
210
211	IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(Nout+1)) THEN
212
213	kbot=ktop
214	DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.N)
215	kbot=kbot+1
216	ENDDO
217	zbox=zbot
218	DO k1= jk+1,Nout
219	IF (z_wout(k1+1)-z_wout(k1).GT.dpthin) THEN
220	exit !thick layer
221	ELSE
222	zbox=z_wout(k1+1) !include thin adjacent layers
223	IF(zbox.EQ.z_wout(Nout+1)) THEN
224	exit !at bottom
225	ENDIF
226	ENDIF
227	ENDDO
228	zthk=zbox-ztop
229
230	kbox=ktop
231	DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.N)
232	kbox=kbox+1
233	ENDDO
234
235	IF(ktop.EQ.kbox) THEN
236	IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN
237	IF(hin(kbox).GT.dpthin) THEN
238	q001 = (zbox-z_win(kbox))/hin(kbox)
239	q002 = (ztop-z_win(kbox))/hin(kbox)
240	q01=q0012+q0022+q001q002+1.-2.(q001+q002)
241	q02=q01-1.+(q001+q002)
242	q0=1.-q01-q02
243	ELSE
244	q0 = 1.0
245	q01 = 0.
246	q02 = 0.
247	ENDIF
248	tabout(jk)=q0coeffremap(kbox,1)+q01coeffremap(kbox,2)+q02*coeffremap(kbox,3)
249	ELSE
250	tabout(jk) = tabin(kbox)
251	ENDIF
252	ELSE
253	IF(ktop.LE.jk .AND. kbox.GE.jk) THEN
254	ka = jk
255	ELSEIF (kbox-ktop.GE.3) THEN
256	ka = (kbox+ktop)/2
257	ELSEIF (hin(ktop).GE.hin(kbox)) THEN
258	ka = ktop
259	ELSE
260	ka = kbox
261	ENDIF !choose ka
262
263	offset=coeffremap(ka,1)
264
265	qtop = z_win(ktop+1)-ztop !partial layer thickness
266	IF(hin(ktop).GT.dpthin) THEN
267	q=(ztop-z_win(ktop))/hin(ktop)
268	q01=q*(q-1.)
269	q02=q01+q
270	q0=1-q01-q02
271	ELSE
272	q0 = 1.
273	q01 = 0.
274	q02 = 0.
275	ENDIF
276
277	tsum =((q0coeffremap(ktop,1)+q01coeffremap(ktop,2)+q02coeffremap(ktop,3))-offset)qtop
278
279	DO k1= ktop+1,kbox-1
280	tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1)
281	ENDDO !k1
282
283	qbot = zbox-z_win(kbox) !partial layer thickness
284	IF(hin(kbox).GT.dpthin) THEN
285	q=qbot/hin(kbox)
286	q01=(q-1.)**2
287	q02=q01-1.+q
288	q0=1-q01-q02
289	ELSE
290	q0 = 1.0
291	q01 = 0.
292	q02 = 0.
293	ENDIF
294
295	tsum = tsum +((q0coeffremap(kbox,1)+q01coeffremap(kbox,2)+q02coeffremap(kbox,3))-offset)qbot
296
297	rpsum=1.0d0/zthk
298	tabout(jk)=offset+tsum*rpsum
299
300	ENDIF !single or multiple layers
301	ELSE
302	IF (jk==1) THEN
303	write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1)
304	ENDIF
305	tabout(jk) = tabout(jk-1)
306
307	ENDIF !normal:thin layer
308	ENDDO !jk
309
310	return
311	end subroutine reconstructandremap
312	#endif
313
314	#endif
315	!!======================================================================
316	END MODULE agrif_oce

Note: See TracBrowser for help on using the repository browser.

Download in other formats: