Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

agrif_oce.F90 @ 9006

Last change on this file since 9006 was 8998, checked in by timgraham, 6 years ago
First commit of Jerome's modified versions of agrif_opa routines
Property svn:keywords set to `Id`
File size: 11.9 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 = .TRUE. !: 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
48	! time-splitting loop:
49	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: ubdy_w, vbdy_w, hbdy_w
50	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: ubdy_e, vbdy_e, hbdy_e
51	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: ubdy_n, vbdy_n, hbdy_n
52	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: ubdy_s, vbdy_s, hbdy_s
53
54	INTEGER :: tsn_id ! AGRIF profile for tracers interpolation and update
55	INTEGER :: un_interp_id, vn_interp_id ! AGRIF profiles for interpolations
56	INTEGER :: un_update_id, vn_update_id ! AGRIF profiles for udpates
57	INTEGER :: tsn_sponge_id, un_sponge_id, vn_sponge_id ! AGRIF profiles for sponge layers
58	# if defined key_top
59	INTEGER :: trn_id, trn_sponge_id
60	# endif
61	INTEGER :: unb_id, vnb_id, ub2b_interp_id, vb2b_interp_id
62	INTEGER :: ub2b_update_id, vb2b_update_id
63	INTEGER :: e3t_id, e1u_id, e2v_id, sshn_id
64	INTEGER :: scales_t_id
65	# if defined key_zdftke \|\| defined key_zdfgls
66	INTEGER :: avt_id, avm_id, en_id
67	# endif
68	INTEGER :: umsk_id, vmsk_id
69	INTEGER :: kindic_agr
70
71	!!----------------------------------------------------------------------
72	!! NEMO/NST 3.3.1 , NEMO Consortium (2011)
73	!! $Id$
74	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
75	!!----------------------------------------------------------------------
76	CONTAINS
77
78	INTEGER FUNCTION agrif_oce_alloc()
79	!!----------------------------------------------------------------------
80	!! * FUNCTION agrif_oce_alloc *
81	!!----------------------------------------------------------------------
82	INTEGER, DIMENSION(2) :: ierr
83	!!----------------------------------------------------------------------
84	ierr(:) = 0
85	!
86	ALLOCATE( fsaht_spu(jpi,jpj), fsaht_spv(jpi,jpj), &
87	& fsahm_spt(jpi,jpj), fsahm_spf(jpi,jpj), &
88	& tabspongedone_tsn(jpi,jpj), &
89	# if defined key_top
90	& tabspongedone_trn(jpi,jpj), &
91	# endif
92	& tabspongedone_u (jpi,jpj), &
93	& tabspongedone_v (jpi,jpj), STAT = ierr(1) )
94
95	ALLOCATE( ubdy_w(jpj), vbdy_w(jpj), hbdy_w(jpj), &
96	& ubdy_e(jpj), vbdy_e(jpj), hbdy_e(jpj), &
97	& ubdy_n(jpi), vbdy_n(jpi), hbdy_n(jpi), &
98	& ubdy_s(jpi), vbdy_s(jpi), hbdy_s(jpi), STAT = ierr(2) )
99
100	agrif_oce_alloc = MAXVAL(ierr)
101	!
102	END FUNCTION agrif_oce_alloc
103
104	#if defined key_vertical
105	SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout)
106	!!----------------------------------------------------------------------
107	!! * FUNCTION reconstructandremap *
108	!!----------------------------------------------------------------------
109	IMPLICIT NONE
110	INTEGER N, Nout
111	REAL(wp) tabin(N), tabout(Nout)
112	REAL(wp) hin(N), hout(Nout)
113	REAL(wp) coeffremap(N,3),zwork(N,3)
114	REAL(wp) zwork2(N+1,3)
115	INTEGER jk
116	DOUBLE PRECISION, PARAMETER :: dsmll=1.0d-8
117	REAL(wp) q,q01,q02,q001,q002,q0
118	REAL(wp) z_win(1:N+1), z_wout(1:Nout+1)
119	REAL(wp),PARAMETER :: dpthin = 1.D-3
120	INTEGER :: k1, kbox, ktop, ka, kbot
121	REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop
122
123	z_win(1)=0.; z_wout(1)= 0.
124	DO jk=1,N
125	z_win(jk+1)=z_win(jk)+hin(jk)
126	ENDDO
127
128	DO jk=1,Nout
129	z_wout(jk+1)=z_wout(jk)+hout(jk)
130	ENDDO
131
132	DO jk=2,N
133	zwork(jk,1)=1./(hin(jk-1)+hin(jk))
134	ENDDO
135
136	DO jk=2,N-1
137	q0 = 1./(hin(jk-1)+hin(jk)+hin(jk+1))
138	zwork(jk,2)=hin(jk-1)+2.*hin(jk)+hin(jk+1)
139	zwork(jk,3)=q0
140	ENDDO
141
142	DO jk= 2,N
143	zwork2(jk,1)=zwork(jk,1)*(tabin(jk)-tabin(jk-1))
144	ENDDO
145
146	coeffremap(:,1) = tabin(:)
147
148	DO jk=2,N-1
149	q001 = hin(jk)*zwork2(jk+1,1)
150	q002 = hin(jk)*zwork2(jk,1)
151	IF (q001*q002 < 0) then
152	q001 = 0.
153	q002 = 0.
154	ENDIF
155	q=zwork(jk,2)
156	q01=q*zwork2(jk+1,1)
157	q02=q*zwork2(jk,1)
158	IF (abs(q001) > abs(q02)) q001 = q02
159	IF (abs(q002) > abs(q01)) q002 = q01
160
161	q=(q001-q002)*zwork(jk,3)
162	q001=q001-q*hin(jk+1)
163	q002=q002+q*hin(jk-1)
164
165	coeffremap(jk,3)=coeffremap(jk,1)+q001
166	coeffremap(jk,2)=coeffremap(jk,1)-q002
167
168	zwork2(jk,1)=(2.q001-q002)*2
169	zwork2(jk,2)=(2.q002-q001)*2
170	ENDDO
171
172	DO jk=1,N
173	IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN
174	coeffremap(jk,3) = coeffremap(jk,1)
175	coeffremap(jk,2) = coeffremap(jk,1)
176	zwork2(jk,1) = 0.
177	zwork2(jk,2) = 0.
178	ENDIF
179	ENDDO
180
181	DO jk=2,N
182	q002=max(zwork2(jk-1,2),dsmll)
183	q001=max(zwork2(jk,1),dsmll)
184	zwork2(jk,3)=(q001coeffremap(jk-1,3)+q002coeffremap(jk,2))/(q001+q002)
185	ENDDO
186
187	zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3)
188	zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3)
189
190	DO jk=1,N
191	q01=zwork2(jk+1,3)-coeffremap(jk,1)
192	q02=coeffremap(jk,1)-zwork2(jk,3)
193	q001=2.*q01
194	q002=2.*q02
195	IF (q01*q02<0) then
196	q01=0.
197	q02=0.
198	ELSEIF (abs(q01)>abs(q002)) then
199	q01=q002
200	ELSEIF (abs(q02)>abs(q001)) then
201	q02=q001
202	ENDIF
203	coeffremap(jk,2)=coeffremap(jk,1)-q02
204	coeffremap(jk,3)=coeffremap(jk,1)+q01
205	ENDDO
206
207	zbot=0.0
208	kbot=1
209	DO jk=1,Nout
210	ztop=zbot !top is bottom of previous layer
211	ktop=kbot
212	IF (ztop.GE.z_win(ktop+1)) then
213	ktop=ktop+1
214	ENDIF
215
216	zbot=z_wout(jk+1)
217	zthk=zbot-ztop
218
219	IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(Nout+1)) THEN
220
221	kbot=ktop
222	DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.N)
223	kbot=kbot+1
224	ENDDO
225	zbox=zbot
226	DO k1= jk+1,Nout
227	IF (z_wout(k1+1)-z_wout(k1).GT.dpthin) THEN
228	exit !thick layer
229	ELSE
230	zbox=z_wout(k1+1) !include thin adjacent layers
231	IF(zbox.EQ.z_wout(Nout+1)) THEN
232	exit !at bottom
233	ENDIF
234	ENDIF
235	ENDDO
236	zthk=zbox-ztop
237
238	kbox=ktop
239	DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.N)
240	kbox=kbox+1
241	ENDDO
242
243	IF(ktop.EQ.kbox) THEN
244	IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN
245	IF(hin(kbox).GT.dpthin) THEN
246	q001 = (zbox-z_win(kbox))/hin(kbox)
247	q002 = (ztop-z_win(kbox))/hin(kbox)
248	q01=q0012+q0022+q001q002+1.-2.(q001+q002)
249	q02=q01-1.+(q001+q002)
250	q0=1.-q01-q02
251	ELSE
252	q0 = 1.0
253	q01 = 0.
254	q02 = 0.
255	ENDIF
256	tabout(jk)=q0coeffremap(kbox,1)+q01coeffremap(kbox,2)+q02*coeffremap(kbox,3)
257	ELSE
258	tabout(jk) = tabin(kbox)
259	ENDIF
260	ELSE
261	IF(ktop.LE.jk .AND. kbox.GE.jk) THEN
262	ka = jk
263	ELSEIF (kbox-ktop.GE.3) THEN
264	ka = (kbox+ktop)/2
265	ELSEIF (hin(ktop).GE.hin(kbox)) THEN
266	ka = ktop
267	ELSE
268	ka = kbox
269	ENDIF !choose ka
270
271	offset=coeffremap(ka,1)
272
273	qtop = z_win(ktop+1)-ztop !partial layer thickness
274	IF(hin(ktop).GT.dpthin) THEN
275	q=(ztop-z_win(ktop))/hin(ktop)
276	q01=q*(q-1.)
277	q02=q01+q
278	q0=1-q01-q02
279	ELSE
280	q0 = 1.
281	q01 = 0.
282	q02 = 0.
283	ENDIF
284
285	tsum =((q0coeffremap(ktop,1)+q01coeffremap(ktop,2)+q02coeffremap(ktop,3))-offset)qtop
286
287	DO k1= ktop+1,kbox-1
288	tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1)
289	ENDDO !k1
290
291	qbot = zbox-z_win(kbox) !partial layer thickness
292	IF(hin(kbox).GT.dpthin) THEN
293	q=qbot/hin(kbox)
294	q01=(q-1.)**2
295	q02=q01-1.+q
296	q0=1-q01-q02
297	ELSE
298	q0 = 1.0
299	q01 = 0.
300	q02 = 0.
301	ENDIF
302
303	tsum = tsum +((q0coeffremap(kbox,1)+q01coeffremap(kbox,2)+q02coeffremap(kbox,3))-offset)qbot
304
305	rpsum=1.0d0/zthk
306	tabout(jk)=offset+tsum*rpsum
307
308	ENDIF !single or multiple layers
309	ELSE
310	IF (jk==1) THEN
311	write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1)
312	ENDIF
313	tabout(jk) = tabout(jk-1)
314
315	ENDIF !normal:thin layer
316	ENDDO !jk
317
318	return
319	end subroutine reconstructandremap
320	#endif
321
322	#endif
323	!!======================================================================
324	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/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90 @ 9006

Download in other formats: