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 @ 11219

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

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

Download in other formats: