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 | # if defined key_top |
---|
62 | INTEGER :: trn_id, trn_sponge_id |
---|
63 | # endif |
---|
64 | INTEGER :: unb_id, vnb_id, ub2b_interp_id, vb2b_interp_id |
---|
65 | INTEGER :: ub2b_update_id, vb2b_update_id |
---|
66 | INTEGER :: e3t_id, e1u_id, e2v_id, sshn_id |
---|
67 | INTEGER :: scales_t_id |
---|
68 | # if defined key_zdftke |
---|
69 | INTEGER :: avt_id, avm_id, en_id |
---|
70 | # endif |
---|
71 | INTEGER :: umsk_id, vmsk_id |
---|
72 | INTEGER :: kindic_agr |
---|
73 | |
---|
74 | !!---------------------------------------------------------------------- |
---|
75 | !! NEMO/NST 3.3.1 , NEMO Consortium (2011) |
---|
76 | !! $Id$ |
---|
77 | !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) |
---|
78 | !!---------------------------------------------------------------------- |
---|
79 | CONTAINS |
---|
80 | |
---|
81 | INTEGER FUNCTION agrif_oce_alloc() |
---|
82 | !!---------------------------------------------------------------------- |
---|
83 | !! *** FUNCTION agrif_oce_alloc *** |
---|
84 | !!---------------------------------------------------------------------- |
---|
85 | INTEGER, DIMENSION(2) :: ierr |
---|
86 | !!---------------------------------------------------------------------- |
---|
87 | ierr(:) = 0 |
---|
88 | ! |
---|
89 | ALLOCATE( fsaht_spu(jpi,jpj), fsaht_spv(jpi,jpj), & |
---|
90 | & fsahm_spt(jpi,jpj), fsahm_spf(jpi,jpj), & |
---|
91 | & tabspongedone_tsn(jpi,jpj), & |
---|
92 | # if defined key_top |
---|
93 | & tabspongedone_trn(jpi,jpj), & |
---|
94 | # endif |
---|
95 | & tabspongedone_u (jpi,jpj), & |
---|
96 | & tabspongedone_v (jpi,jpj), STAT = ierr(1) ) |
---|
97 | |
---|
98 | ALLOCATE( ubdy_w(jpj), vbdy_w(jpj), hbdy_w(jpj), & |
---|
99 | & ubdy_e(jpj), vbdy_e(jpj), hbdy_e(jpj), & |
---|
100 | & ubdy_n(jpi), vbdy_n(jpi), hbdy_n(jpi), & |
---|
101 | & ubdy_s(jpi), vbdy_s(jpi), hbdy_s(jpi), STAT = ierr(2) ) |
---|
102 | |
---|
103 | agrif_oce_alloc = MAXVAL(ierr) |
---|
104 | ! |
---|
105 | END FUNCTION agrif_oce_alloc |
---|
106 | |
---|
107 | subroutine reconstructandremap(tabin,hin,tabout,hout,N,Nout) |
---|
108 | implicit none |
---|
109 | integer N, Nout |
---|
110 | real tabin(N), tabout(Nout) |
---|
111 | real hin(N), hout(Nout) |
---|
112 | real coeffremap(N,3),zwork(N,3) |
---|
113 | real zwork2(N+1,3) |
---|
114 | integer k |
---|
115 | double precision, parameter :: dsmll=1.0d-8 |
---|
116 | real q,q01,q02,q001,q002,q0 |
---|
117 | real z_win(1:N+1), z_wout(1:Nout+1) |
---|
118 | real,parameter :: dpthin = 1.D-3 |
---|
119 | integer :: k1, kbox, ktop, ka, kbot |
---|
120 | real :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop |
---|
121 | |
---|
122 | z_win(1)=0.; z_wout(1)= 0. |
---|
123 | do k=1,N |
---|
124 | z_win(k+1)=z_win(k)+hin(k) |
---|
125 | enddo |
---|
126 | |
---|
127 | do k=1,Nout |
---|
128 | z_wout(k+1)=z_wout(k)+hout(k) |
---|
129 | enddo |
---|
130 | |
---|
131 | do k=2,N |
---|
132 | zwork(k,1)=1./(hin(k-1)+hin(k)) |
---|
133 | enddo |
---|
134 | |
---|
135 | do k=2,N-1 |
---|
136 | q0 = 1./(hin(k-1)+hin(k)+hin(k+1)) |
---|
137 | zwork(k,2)=hin(k-1)+2.*hin(k)+hin(k+1) |
---|
138 | zwork(k,3)=q0 |
---|
139 | enddo |
---|
140 | |
---|
141 | do k= 2,N |
---|
142 | zwork2(k,1)=zwork(k,1)*(tabin(k)-tabin(k-1)) |
---|
143 | enddo |
---|
144 | |
---|
145 | coeffremap(:,1) = tabin(:) |
---|
146 | |
---|
147 | do k=2,N-1 |
---|
148 | q001 = hin(k)*zwork2(k+1,1) |
---|
149 | q002 = hin(k)*zwork2(k,1) |
---|
150 | if (q001*q002 < 0) then |
---|
151 | q001 = 0. |
---|
152 | q002 = 0. |
---|
153 | endif |
---|
154 | q=zwork(k,2) |
---|
155 | q01=q*zwork2(k+1,1) |
---|
156 | q02=q*zwork2(k,1) |
---|
157 | if (abs(q001) > abs(q02)) q001 = q02 |
---|
158 | if (abs(q002) > abs(q01)) q002 = q01 |
---|
159 | |
---|
160 | q=(q001-q002)*zwork(k,3) |
---|
161 | q001=q001-q*hin(k+1) |
---|
162 | q002=q002+q*hin(k-1) |
---|
163 | |
---|
164 | coeffremap(k,3)=coeffremap(k,1)+q001 |
---|
165 | coeffremap(k,2)=coeffremap(k,1)-q002 |
---|
166 | |
---|
167 | zwork2(k,1)=(2.*q001-q002)**2 |
---|
168 | zwork2(k,2)=(2.*q002-q001)**2 |
---|
169 | enddo |
---|
170 | |
---|
171 | do k=1,N |
---|
172 | if (k.eq.1 .or. k.eq.N .or. hin(k).le.dpthin) then |
---|
173 | coeffremap(k,3) = coeffremap(k,1) |
---|
174 | coeffremap(k,2) = coeffremap(k,1) |
---|
175 | zwork2(k,1) = 0. |
---|
176 | zwork2(k,2) = 0. |
---|
177 | endif |
---|
178 | enddo |
---|
179 | |
---|
180 | do k=2,N |
---|
181 | q002=max(zwork2(k-1,2),dsmll) |
---|
182 | q001=max(zwork2(k,1),dsmll) |
---|
183 | zwork2(k,3)=(q001*coeffremap(k-1,3)+q002*coeffremap(k,2))/(q001+q002) |
---|
184 | enddo |
---|
185 | |
---|
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 k=1,N |
---|
190 | q01=zwork2(k+1,3)-coeffremap(k,1) |
---|
191 | q02=coeffremap(k,1)-zwork2(k,3) |
---|
192 | q001=2.*q01 |
---|
193 | q002=2.*q02 |
---|
194 | if (q01*q02<0) then |
---|
195 | q01=0. |
---|
196 | q02=0. |
---|
197 | elseif (abs(q01)>abs(q002)) then |
---|
198 | q01=q002 |
---|
199 | elseif (abs(q02)>abs(q001)) then |
---|
200 | q02=q001 |
---|
201 | endif |
---|
202 | coeffremap(k,2)=coeffremap(k,1)-q02 |
---|
203 | coeffremap(k,3)=coeffremap(k,1)+q01 |
---|
204 | enddo |
---|
205 | |
---|
206 | zbot=0.0 |
---|
207 | kbot=1 |
---|
208 | do k=1,Nout |
---|
209 | ztop=zbot !top is bottom of previous layer |
---|
210 | ktop=kbot |
---|
211 | if (ztop.ge.z_win(ktop+1)) then |
---|
212 | ktop=ktop+1 |
---|
213 | endif |
---|
214 | |
---|
215 | zbot=z_wout(k+1) |
---|
216 | zthk=zbot-ztop |
---|
217 | |
---|
218 | if (zthk.gt.dpthin .and. ztop.lt.z_wout(Nout+1)) then |
---|
219 | |
---|
220 | kbot=ktop |
---|
221 | do while (z_win(kbot+1).lt.zbot.and.kbot.lt.N) |
---|
222 | kbot=kbot+1 |
---|
223 | enddo |
---|
224 | zbox=zbot |
---|
225 | do k1= k+1,Nout |
---|
226 | if (z_wout(k1+1)-z_wout(k1).gt.dpthin) then |
---|
227 | exit !thick layer |
---|
228 | else |
---|
229 | zbox=z_wout(k1+1) !include thin adjacent layers |
---|
230 | if (zbox.eq.z_wout(Nout+1)) then |
---|
231 | exit !at bottom |
---|
232 | endif |
---|
233 | endif |
---|
234 | enddo |
---|
235 | zthk=zbox-ztop |
---|
236 | |
---|
237 | kbox=ktop |
---|
238 | do while (z_win(kbox+1).lt.zbox.and.kbox.lt.N) |
---|
239 | kbox=kbox+1 |
---|
240 | enddo |
---|
241 | |
---|
242 | if (ktop.eq.kbox) then |
---|
243 | |
---|
244 | |
---|
245 | if (z_wout(k) .ne.z_win(kbox) .or.z_wout(k+1).ne.z_win(kbox+1) ) then |
---|
246 | |
---|
247 | if (hin(kbox).gt.dpthin) then |
---|
248 | q001 = (zbox-z_win(kbox))/hin(kbox) |
---|
249 | q002 = (ztop-z_win(kbox))/hin(kbox) |
---|
250 | q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002) |
---|
251 | q02=q01-1.+(q001+q002) |
---|
252 | q0=1.-q01-q02 |
---|
253 | else |
---|
254 | q0 = 1.0 |
---|
255 | q01 = 0. |
---|
256 | q02 = 0. |
---|
257 | endif |
---|
258 | tabout(k)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3) |
---|
259 | |
---|
260 | else |
---|
261 | tabout(k) = tabin(kbox) |
---|
262 | |
---|
263 | endif |
---|
264 | |
---|
265 | else |
---|
266 | |
---|
267 | if (ktop.le.k .and. kbox.ge.k) then |
---|
268 | ka = k |
---|
269 | elseif (kbox-ktop.ge.3) then |
---|
270 | ka = (kbox+ktop)/2 |
---|
271 | elseif (hin(ktop).ge.hin(kbox)) then |
---|
272 | ka = ktop |
---|
273 | else |
---|
274 | ka = kbox |
---|
275 | endif !choose ka |
---|
276 | |
---|
277 | offset=coeffremap(ka,1) |
---|
278 | |
---|
279 | qtop = z_win(ktop+1)-ztop !partial layer thickness |
---|
280 | if (hin(ktop).gt.dpthin) then |
---|
281 | q=(ztop-z_win(ktop))/hin(ktop) |
---|
282 | q01=q*(q-1.) |
---|
283 | q02=q01+q |
---|
284 | q0=1-q01-q02 |
---|
285 | else |
---|
286 | q0 = 1. |
---|
287 | q01 = 0. |
---|
288 | q02 = 0. |
---|
289 | endif |
---|
290 | |
---|
291 | tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop |
---|
292 | |
---|
293 | do k1= ktop+1,kbox-1 |
---|
294 | tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1) |
---|
295 | enddo !k1 |
---|
296 | |
---|
297 | |
---|
298 | qbot = zbox-z_win(kbox) !partial layer thickness |
---|
299 | if (hin(kbox).gt.dpthin) then |
---|
300 | q=qbot/hin(kbox) |
---|
301 | q01=(q-1.)**2 |
---|
302 | q02=q01-1.+q |
---|
303 | q0=1-q01-q02 |
---|
304 | else |
---|
305 | q0 = 1.0 |
---|
306 | q01 = 0. |
---|
307 | q02 = 0. |
---|
308 | endif |
---|
309 | |
---|
310 | tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot |
---|
311 | |
---|
312 | rpsum=1.0d0/zthk |
---|
313 | tabout(k)=offset+tsum*rpsum |
---|
314 | |
---|
315 | endif !single or multiple layers |
---|
316 | else |
---|
317 | if (k==1) then |
---|
318 | write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(k+1),hout(1) |
---|
319 | endif |
---|
320 | tabout(k) = tabout(k-1) |
---|
321 | |
---|
322 | endif !normal:thin layer |
---|
323 | enddo !k |
---|
324 | |
---|
325 | return |
---|
326 | end subroutine reconstructandremap |
---|
327 | |
---|
328 | #endif |
---|
329 | !!====================================================================== |
---|
330 | END MODULE agrif_oce |
---|