Context Navigation

agrif_oce.F90 @ 6258

Last change on this file since 6258 was 6258, checked in by timgraham, 8 years ago

First inclusion of Laurent Debreu's modified code for vertical refinement.
Still a lot of outstanding issues:
1) conv preprocessor fails for limrhg.F90 at the moment (for now I've run without ice model)
2) conv preprocessor fails for STO code - removed this code from testing for now
3) conv preprocessor fails for cpl_oasis.F90 - can work round this by modifying code but the preprocessor should be fixed to deal with this.

After that code compiles and can be run for horizontal grid refinement. Not yet working for vertical refinement.

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	PUBLIC reconstructandremap
20
21	! !!* Namelist namagrif: AGRIF parameters
22	LOGICAL , PUBLIC :: ln_spc_dyn = .FALSE. !:
23	INTEGER , PUBLIC :: nn_cln_update = 3 !: update frequency
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
29	! !!! OLD namelist names
30	INTEGER , PUBLIC :: nbcline = 0 !: update counter
31	INTEGER , PUBLIC :: nbclineupdate !: update frequency
32	REAL(wp), PUBLIC :: visc_tra !: sponge coeff. for tracers
33	REAL(wp), PUBLIC :: visc_dyn !: sponge coeff. for dynamics
34
35	LOGICAL , PUBLIC :: spongedoneT = .FALSE. !: tracer sponge layer indicator
36	LOGICAL , PUBLIC :: spongedoneU = .FALSE. !: dynamics sponge layer indicator
37	LOGICAL , PUBLIC :: lk_agrif_fstep = .TRUE. !: if true: first step
38	LOGICAL , PUBLIC :: lk_agrif_doupd = .TRUE. !: if true: send update from current grid
39	LOGICAL , PUBLIC :: lk_agrif_debug = .FALSE. !: if true: print debugging info
40
41	LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_tsn
42	# if defined key_top
43	LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_trn
44	# endif
45	LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_u
46	LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_v
47	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fsaht_spu, fsaht_spv !: sponge diffusivities
48	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fsahm_spt, fsahm_spf !: sponge viscosities
49
50	! Barotropic arrays used to store open boundary data during
51	! time-splitting loop:
52	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: ubdy_w, vbdy_w, hbdy_w
53	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: ubdy_e, vbdy_e, hbdy_e
54	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: ubdy_n, vbdy_n, hbdy_n
55	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: ubdy_s, vbdy_s, hbdy_s
56
57	INTEGER :: tsn_id ! AGRIF profile for tracers interpolation and update
58	INTEGER :: un_interp_id, vn_interp_id ! AGRIF profiles for interpolations
59	INTEGER :: un_update_id, vn_update_id ! AGRIF profiles for udpates
60	INTEGER :: tsn_sponge_id, un_sponge_id, vn_sponge_id ! AGRIF profiles for sponge layers
61	! VERTICAL REFINEMENT BEGIN
62	INTEGER :: scales_t_id, scales_u_id, scales_v_id
63	! VERTICAL REFINEMENT END
64
65	# if defined key_top
66	INTEGER :: trn_id, trn_sponge_id
67	# endif
68	INTEGER :: unb_id, vnb_id, ub2b_interp_id, vb2b_interp_id
69	INTEGER :: ub2b_update_id, vb2b_update_id
70	INTEGER :: e3t_id, e1u_id, e2v_id, sshn_id
71	# if defined key_zdftke
72	INTEGER :: avt_id, avm_id, en_id
73	# endif
74	INTEGER :: umsk_id, vmsk_id
75	INTEGER :: kindic_agr
76
77	!!----------------------------------------------------------------------
78	!! NEMO/NST 3.3.1 , NEMO Consortium (2011)
79	!! $Id: agrif_oce.F90 5081 2015-02-13 09:51:27Z smasson $
80	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
81	!!----------------------------------------------------------------------
82	CONTAINS
83
84	INTEGER FUNCTION agrif_oce_alloc()
85	!!----------------------------------------------------------------------
86	!! * FUNCTION agrif_oce_alloc *
87	!!----------------------------------------------------------------------
88	INTEGER, DIMENSION(2) :: ierr
89	!!----------------------------------------------------------------------
90	ierr(:) = 0
91	!
92	ALLOCATE( fsaht_spu(jpi,jpj), fsaht_spv(jpi,jpj), &
93	& fsahm_spt(jpi,jpj), fsahm_spf(jpi,jpj), &
94	& tabspongedone_tsn(jpi,jpj), &
95	# if defined key_top
96	& tabspongedone_trn(jpi,jpj), &
97	# endif
98	& tabspongedone_u (jpi,jpj), &
99	& tabspongedone_v (jpi,jpj), STAT = ierr(1) )
100
101	ALLOCATE( ubdy_w(jpj), vbdy_w(jpj), hbdy_w(jpj), &
102	& ubdy_e(jpj), vbdy_e(jpj), hbdy_e(jpj), &
103	& ubdy_n(jpi), vbdy_n(jpi), hbdy_n(jpi), &
104	& ubdy_s(jpi), vbdy_s(jpi), hbdy_s(jpi), STAT = ierr(2) )
105
106	agrif_oce_alloc = MAXVAL(ierr)
107	!
108	END FUNCTION agrif_oce_alloc
109
110	subroutine reconstructandremap(tabin,hin,tabout,hout,N,Nout)
111	implicit none
112	integer N, Nout
113	real tabin(N), tabout(Nout)
114	real hin(N), hout(Nout)
115	real coeffremap(N,3),zwork(N,3)
116	real zwork2(N+1,3)
117	integer k
118	double precision, parameter :: dsmll=1.0d-8
119	real q,q01,q02,q001,q002,q0
120	real z_win(1:N+1), z_wout(1:Nout+1)
121	real,parameter :: dpthin = 1.D-3
122	integer :: k1, kbox, ktop, ka, kbot
123	real :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop
124
125	z_win(1)=0.; z_wout(1)= 0.
126	do k=1,N
127	z_win(k+1)=z_win(k)+hin(k)
128	enddo
129
130	do k=1,Nout
131	z_wout(k+1)=z_wout(k)+hout(k)
132	enddo
133
134	do k=2,N
135	zwork(k,1)=1./(hin(k-1)+hin(k))
136	enddo
137
138	do k=2,N-1
139	q0 = 1./(hin(k-1)+hin(k)+hin(k+1))
140	zwork(k,2)=hin(k-1)+2.*hin(k)+hin(k+1)
141	zwork(k,3)=q0
142	enddo
143
144	do k= 2,N
145	zwork2(k,1)=zwork(k,1)*(tabin(k)-tabin(k-1))
146	enddo
147
148	coeffremap(:,1) = tabin(:)
149
150	do k=2,N-1
151	q001 = hin(k)*zwork2(k+1,1)
152	q002 = hin(k)*zwork2(k,1)
153	if (q001*q002 < 0) then
154	q001 = 0.
155	q002 = 0.
156	endif
157	q=zwork(k,2)
158	q01=q*zwork2(k+1,1)
159	q02=q*zwork2(k,1)
160	if (abs(q001) > abs(q02)) q001 = q02
161	if (abs(q002) > abs(q01)) q002 = q01
162
163	q=(q001-q002)*zwork(k,3)
164	q001=q001-q*hin(k+1)
165	q002=q002+q*hin(k-1)
166
167	coeffremap(k,3)=coeffremap(k,1)+q001
168	coeffremap(k,2)=coeffremap(k,1)-q002
169
170	zwork2(k,1)=(2.q001-q002)*2
171	zwork2(k,2)=(2.q002-q001)*2
172	enddo
173
174	do k=1,N
175	if (k.eq.1 .or. k.eq.N .or. hin(k).le.dpthin) then
176	coeffremap(k,3) = coeffremap(k,1)
177	coeffremap(k,2) = coeffremap(k,1)
178	zwork2(k,1) = 0.
179	zwork2(k,2) = 0.
180	endif
181	enddo
182
183	do k=2,N
184	q002=max(zwork2(k-1,2),dsmll)
185	q001=max(zwork2(k,1),dsmll)
186	zwork2(k,3)=(q001coeffremap(k-1,3)+q002coeffremap(k,2))/(q001+q002)
187	enddo
188
189	zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3)
190	zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3)
191
192	do k=1,N
193	q01=zwork2(k+1,3)-coeffremap(k,1)
194	q02=coeffremap(k,1)-zwork2(k,3)
195	q001=2.*q01
196	q002=2.*q02
197	if (q01*q02<0) then
198	q01=0.
199	q02=0.
200	elseif (abs(q01)>abs(q002)) then
201	q01=q002
202	elseif (abs(q02)>abs(q001)) then
203	q02=q001
204	endif
205	coeffremap(k,2)=coeffremap(k,1)-q02
206	coeffremap(k,3)=coeffremap(k,1)+q01
207	enddo
208
209	zbot=0.0
210	kbot=1
211	do k=1,Nout
212	ztop=zbot !top is bottom of previous layer
213	ktop=kbot
214	if (ztop.ge.z_win(ktop+1)) then
215	ktop=ktop+1
216	endif
217
218	zbot=z_wout(k+1)
219	zthk=zbot-ztop
220
221	if (zthk.gt.dpthin .and. ztop.lt.z_wout(Nout+1)) then
222
223	kbot=ktop
224	do while (z_win(kbot+1).lt.zbot.and.kbot.lt.N)
225	kbot=kbot+1
226	enddo
227	zbox=zbot
228	do k1= k+1,Nout
229	if (z_wout(k1+1)-z_wout(k1).gt.dpthin) then
230	exit !thick layer
231	else
232	zbox=z_wout(k1+1) !include thin adjacent layers
233	if (zbox.eq.z_wout(Nout+1)) then
234	exit !at bottom
235	endif
236	endif
237	enddo
238	zthk=zbox-ztop
239
240	kbox=ktop
241	do while (z_win(kbox+1).lt.zbox.and.kbox.lt.N)
242	kbox=kbox+1
243	enddo
244
245	if (ktop.eq.kbox) then
246
247
248	if (z_wout(k) .ne.z_win(kbox) .or.z_wout(k+1).ne.z_win(kbox+1) ) then
249
250	if (hin(kbox).gt.dpthin) then
251	q001 = (zbox-z_win(kbox))/hin(kbox)
252	q002 = (ztop-z_win(kbox))/hin(kbox)
253	q01=q0012+q0022+q001q002+1.-2.(q001+q002)
254	q02=q01-1.+(q001+q002)
255	q0=1.-q01-q02
256	else
257	q0 = 1.0
258	q01 = 0.
259	q02 = 0.
260	endif
261	tabout(k)=q0coeffremap(kbox,1)+q01coeffremap(kbox,2)+q02*coeffremap(kbox,3)
262
263	else
264	tabout(k) = tabin(kbox)
265
266	endif
267
268	else
269
270	if (ktop.le.k .and. kbox.ge.k) then
271	ka = k
272	elseif (kbox-ktop.ge.3) then
273	ka = (kbox+ktop)/2
274	elseif (hin(ktop).ge.hin(kbox)) then
275	ka = ktop
276	else
277	ka = kbox
278	endif !choose ka
279
280	offset=coeffremap(ka,1)
281
282	qtop = z_win(ktop+1)-ztop !partial layer thickness
283	if (hin(ktop).gt.dpthin) then
284	q=(ztop-z_win(ktop))/hin(ktop)
285	q01=q*(q-1.)
286	q02=q01+q
287	q0=1-q01-q02
288	else
289	q0 = 1.
290	q01 = 0.
291	q02 = 0.
292	endif
293
294	tsum =((q0coeffremap(ktop,1)+q01coeffremap(ktop,2)+q02coeffremap(ktop,3))-offset)qtop
295
296	do k1= ktop+1,kbox-1
297	tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1)
298	enddo !k1
299
300
301	qbot = zbox-z_win(kbox) !partial layer thickness
302	if (hin(kbox).gt.dpthin) then
303	q=qbot/hin(kbox)
304	q01=(q-1.)**2
305	q02=q01-1.+q
306	q0=1-q01-q02
307	else
308	q0 = 1.0
309	q01 = 0.
310	q02 = 0.
311	endif
312
313	tsum = tsum +((q0coeffremap(kbox,1)+q01coeffremap(kbox,2)+q02coeffremap(kbox,3))-offset)qbot
314
315	rpsum=1.0d0/zthk
316	tabout(k)=offset+tsum*rpsum
317
318	endif !single or multiple layers
319	else
320	if (k==1) then
321	write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(k+1),hout(1)
322	endif
323	tabout(k) = tabout(k-1)
324
325	endif !normal:thin layer
326	enddo !k
327
328	return
329	end subroutine reconstructandremap
330
331	#endif
332	!!======================================================================
333	END MODULE agrif_oce

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

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

Context Navigation

source: branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90 @ 6258

Download in other formats: