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.

p4zsink.F90 in trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago
update trunk with OpenMP parallelization
File size: 21.6 KB

Line
1	MODULE p4zsink
2	!!======================================================================
3	!! * MODULE p4zsink *
4	!! TOP : PISCES vertical flux of particulate matter due to gravitational sinking
5	!!======================================================================
6	!! History : 1.0 ! 2004 (O. Aumont) Original code
7	!! 2.0 ! 2007-12 (C. Ethe, G. Madec) F90
8	!! 3.4 ! 2011-06 (O. Aumont, C. Ethe) Change aggregation formula
9	!! 3.5 ! 2012-07 (O. Aumont) Introduce potential time-splitting
10	!!----------------------------------------------------------------------
11	!! p4z_sink : Compute vertical flux of particulate matter due to gravitational sinking
12	!! p4z_sink_init : Unitialisation of sinking speed parameters
13	!! p4z_sink_alloc : Allocate sinking speed variables
14	!!----------------------------------------------------------------------
15	USE oce_trc ! shared variables between ocean and passive tracers
16	USE trc ! passive tracers common variables
17	USE sms_pisces ! PISCES Source Minus Sink variables
18	USE prtctl_trc ! print control for debugging
19	USE iom ! I/O manager
20	USE lib_mpp
21
22	IMPLICIT NONE
23	PRIVATE
24
25	PUBLIC p4z_sink ! called in p4zbio.F90
26	PUBLIC p4z_sink_init ! called in trcsms_pisces.F90
27	PUBLIC p4z_sink_alloc
28
29	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinking, sinking2 !: POC sinking fluxes
30	! ! (different meanings depending on the parameterization)
31	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkingn, sinking2n !: POC sinking fluxes
32	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkingp, sinking2p !: POC sinking fluxes
33	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkcal, sinksil !: CaCO3 and BSi sinking fluxes
34	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkfer !: Small BFe sinking fluxes
35	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkfer2 !: Big iron sinking fluxes
36	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkfep !: Fep sinking fluxes
37
38	INTEGER :: ik100
39
40	!!----------------------------------------------------------------------
41	!! NEMO/TOP 3.3 , NEMO Consortium (2010)
42	!! $Id: p4zsink.F90 3160 2011-11-20 14:27:18Z cetlod $
43	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
44	!!----------------------------------------------------------------------
45	CONTAINS
46
47	!!----------------------------------------------------------------------
48	!! 'standard sinking parameterisation' ???
49	!!----------------------------------------------------------------------
50
51	SUBROUTINE p4z_sink ( kt, knt )
52	!!---------------------------------------------------------------------
53	!! * ROUTINE p4z_sink *
54	!!
55	!! ** Purpose : Compute vertical flux of particulate matter due to
56	!! gravitational sinking
57	!!
58	!! ** Method : - ???
59	!!---------------------------------------------------------------------
60	INTEGER, INTENT(in) :: kt, knt
61	INTEGER :: ji, jj, jk, jit
62	INTEGER :: iiter1, iiter2
63	REAL(wp) :: zagg1, zagg2, zagg3, zagg4
64	REAL(wp) :: zagg , zaggfe, zaggdoc, zaggdoc2, zaggdoc3
65	REAL(wp) :: zfact, zwsmax, zmax
66	CHARACTER (len=25) :: charout
67	REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d
68	REAL(wp), POINTER, DIMENSION(:,: ) :: zw2d
69	!!---------------------------------------------------------------------
70	!
71	IF( nn_timing == 1 ) CALL timing_start('p4z_sink')
72
73
74	! Initialization of some global variables
75	! ---------------------------------------
76	!$OMP PARALLEL
77	!$OMP DO schedule(static) private(jk, jj, ji)
78	DO jk = 1, jpk
79	DO jj = 1, jpj
80	DO ji = 1,jpi
81	prodpoc(ji,jj,jk) = 0.
82	conspoc(ji,jj,jk) = 0.
83	prodgoc(ji,jj,jk) = 0.
84	consgoc(ji,jj,jk) = 0.
85	END DO
86	END DO
87	END DO
88
89	!
90	! Sinking speeds of detritus is increased with depth as shown
91	! by data and from the coagulation theory
92	! -----------------------------------------------------------
93	!$OMP DO schedule(static) private(jk, jj, ji, zmax, zfact)
94	DO jk = 1, jpkm1
95	DO jj = 1, jpj
96	DO ji = 1,jpi
97	zmax = MAX( heup_01(ji,jj), hmld(ji,jj) )
98	zfact = MAX( 0., gdepw_n(ji,jj,jk+1) - zmax ) / wsbio2scale
99	wsbio4(ji,jj,jk) = wsbio2 + MAX(0., ( wsbio2max - wsbio2 )) * zfact
100	END DO
101	END DO
102	END DO
103
104	! limit the values of the sinking speeds to avoid numerical instabilities
105	!$OMP DO schedule(static) private(jk, jj, ji)
106	DO jk = 1, jpk
107	DO jj = 1, jpj
108	DO ji = 1, jpi
109	wsbio3(ji,jj,jk) = wsbio
110	END DO
111	END DO
112	END DO
113	!$OMP END PARALLEL
114
115	!
116	! OA This is (I hope) a temporary solution for the problem that may
117	! OA arise in specific situation where the CFL criterion is broken
118	! OA for vertical sedimentation of particles. To avoid this, a time
119	! OA splitting algorithm has been coded. A specific maximum
120	! OA iteration number is provided and may be specified in the namelist
121	! OA This is to avoid very large iteration number when explicit free
122	! OA surface is used (for instance). When niter?max is set to 1,
123	! OA this computation is skipped. The crude old threshold method is
124	! OA then applied. This also happens when niter exceeds nitermax.
125	IF( MAX( niter1max, niter2max ) == 1 ) THEN
126	iiter1 = 1
127	iiter2 = 1
128	ELSE
129	iiter1 = 1
130	iiter2 = 1
131	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,zwsmax) REDUCTION(MAX:iiter1, iiter2)
132	DO jk = 1, jpkm1
133	DO jj = 1, jpj
134	DO ji = 1, jpi
135	IF( tmask(ji,jj,jk) == 1) THEN
136	zwsmax = 0.5 * e3t_n(ji,jj,jk) / xstep
137	iiter1 = MAX( iiter1, INT( wsbio3(ji,jj,jk) / zwsmax ) )
138	iiter2 = MAX( iiter2, INT( wsbio4(ji,jj,jk) / zwsmax ) )
139	ENDIF
140	END DO
141	END DO
142	END DO
143	IF( lk_mpp ) THEN
144	CALL mpp_max( iiter1 )
145	CALL mpp_max( iiter2 )
146	ENDIF
147	iiter1 = MIN( iiter1, niter1max )
148	iiter2 = MIN( iiter2, niter2max )
149	ENDIF
150
151	!$OMP PARALLEL
152	!$OMP DO schedule(static) private(jk, jj, ji, zwsmax)
153	DO jk = 1,jpkm1
154	DO jj = 1, jpj
155	DO ji = 1, jpi
156	IF( tmask(ji,jj,jk) == 1 ) THEN
157	zwsmax = 0.5 * e3t_n(ji,jj,jk) / xstep
158	wsbio3(ji,jj,jk) = MIN( wsbio3(ji,jj,jk), zwsmax * REAL( iiter1, wp ) )
159	wsbio4(ji,jj,jk) = MIN( wsbio4(ji,jj,jk), zwsmax * REAL( iiter2, wp ) )
160	ENDIF
161	END DO
162	END DO
163	END DO
164
165	! Initializa to zero all the sinking arrays
166	! -----------------------------------------
167	!$OMP DO schedule(static) private(jk, jj, ji)
168	DO jk = 1, jpk
169	DO jj = 1, jpj
170	DO ji = 1, jpi
171	sinking (ji,jj,jk) = 0.e0
172	sinking2(ji,jj,jk) = 0.e0
173	sinkcal (ji,jj,jk) = 0.e0
174	sinkfer (ji,jj,jk) = 0.e0
175	sinksil (ji,jj,jk) = 0.e0
176	sinkfer2(ji,jj,jk) = 0.e0
177	wscal (ji,jj,jk) = wsbio4(ji,jj,jk)
178	END DO
179	END DO
180	END DO
181	!$OMP END PARALLEL
182
183	! Compute the sedimentation term using p4zsink2 for all the sinking particles
184	! -----------------------------------------------------
185	DO jit = 1, iiter1
186	CALL p4z_sink2( wsbio3, sinking , jppoc, iiter1 )
187	CALL p4z_sink2( wsbio3, sinkfer , jpsfe, iiter1 )
188	END DO
189
190	DO jit = 1, iiter2
191	CALL p4z_sink2( wsbio4, sinking2, jpgoc, iiter2 )
192	CALL p4z_sink2( wsbio4, sinkfer2, jpbfe, iiter2 )
193	CALL p4z_sink2( wsbio4, sinksil , jpgsi, iiter2 )
194	CALL p4z_sink2( wscal , sinkcal , jpcal, iiter2 )
195	END DO
196
197	IF( ln_p5z ) THEN
198	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
199	DO jk = 1, jpk
200	DO jj = 1, jpj
201	DO ji = 1, jpi
202	sinkingn (ji,jj,jk) = 0.e0
203	sinking2n(ji,jj,jk) = 0.e0
204	sinkingp (ji,jj,jk) = 0.e0
205	sinking2p(ji,jj,jk) = 0.e0
206	END DO
207	END DO
208	END DO
209
210	! Compute the sedimentation term using p4zsink2 for all the sinking particles
211	! -----------------------------------------------------
212	DO jit = 1, iiter1
213	CALL p4z_sink2( wsbio3, sinkingn , jppon, iiter1 )
214	CALL p4z_sink2( wsbio3, sinkingp , jppop, iiter1 )
215	END DO
216
217	DO jit = 1, iiter2
218	CALL p4z_sink2( wsbio4, sinking2n, jpgon, iiter2 )
219	CALL p4z_sink2( wsbio4, sinking2p, jpgop, iiter2 )
220	END DO
221	ENDIF
222
223	IF( ln_ligand ) THEN
224	!$OMP PARALLEL
225	!$OMP DO schedule(static) private(jk, jj, ji)
226	DO jk = 1, jpk
227	DO jj = 1, jpj
228	DO ji = 1, jpi
229	wsfep (ji,jj,jk) = wfep
230	END DO
231	END DO
232	END DO
233	!$OMP DO schedule(static) private(jk, jj, ji, zwsmax)
234	DO jk = 1,jpkm1
235	DO jj = 1, jpj
236	DO ji = 1, jpi
237	IF( tmask(ji,jj,jk) == 1 ) THEN
238	zwsmax = 0.5 * e3t_n(ji,jj,jk) / xstep
239	wsfep(ji,jj,jk) = MIN( wsfep(ji,jj,jk), zwsmax * REAL( iiter1, wp ) )
240	ENDIF
241	END DO
242	END DO
243	END DO
244	!$OMP END DO NOWAIT
245	!
246	!$OMP DO schedule(static) private(jk, jj, ji)
247	DO jk = 1, jpk
248	DO jj = 1, jpj
249	DO ji = 1, jpi
250	sinkfep(ji,jj,jk) = 0.e0
251	END DO
252	END DO
253	END DO
254	!$OMP END PARALLEL
255	DO jit = 1, iiter1
256	CALL p4z_sink2( wsfep, sinkfep , jpfep, iiter1 )
257	END DO
258	ENDIF
259
260	! Total carbon export per year
261	IF( iom_use( "tcexp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc ) ) &
262	& t_oce_co2_exp = glob_sum( ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * e1e2t(:,:) * tmask(:,:,1) )
263	!
264	IF( lk_iomput ) THEN
265	IF( knt == nrdttrc ) THEN
266	CALL wrk_alloc( jpi, jpj, zw2d )
267	CALL wrk_alloc( jpi, jpj, jpk, zw3d )
268	zfact = 1.e+3 * rfact2r ! conversion from mol/l/kt to mol/m3/s
269	!
270	IF( iom_use( "EPC100" ) ) THEN
271	!$OMP PARALLEL DO schedule(static) private(jj, ji)
272	DO jj = 1, jpj
273	DO ji = 1, jpi
274	zw2d(ji,jj) = ( sinking(ji,jj,ik100) + sinking2(ji,jj,ik100) ) * zfact * tmask(ji,jj,1) ! Export of carbon at 100m
275	END DO
276	END DO
277	CALL iom_put( "EPC100" , zw2d )
278	ENDIF
279	IF( iom_use( "EPFE100" ) ) THEN
280	!$OMP PARALLEL DO schedule(static) private(jj, ji)
281	DO jj = 1, jpj
282	DO ji = 1, jpi
283	zw2d(ji,jj) = ( sinkfer(ji,jj,ik100) + sinkfer2(ji,jj,ik100) ) * zfact * tmask(ji,jj,1) ! Export of iron at 100m
284	END DO
285	END DO
286	CALL iom_put( "EPFE100" , zw2d )
287	ENDIF
288	IF( iom_use( "EPCAL100" ) ) THEN
289	!$OMP PARALLEL DO schedule(static) private(jj, ji)
290	DO jj = 1, jpj
291	DO ji = 1, jpi
292	zw2d(ji,jj) = sinkcal(ji,jj,ik100) * zfact * tmask(ji,jj,1) ! Export of calcite at 100m
293	END DO
294	END DO
295	CALL iom_put( "EPCAL100" , zw2d )
296	ENDIF
297	IF( iom_use( "EPSI100" ) ) THEN
298	!$OMP PARALLEL DO schedule(static) private(jj, ji)
299	DO jj = 1, jpj
300	DO ji = 1, jpi
301	zw2d(ji,jj) = sinksil(ji,jj,ik100) * zfact * tmask(ji,jj,1) ! Export of bigenic silica at 100m
302	END DO
303	END DO
304	CALL iom_put( "EPSI100" , zw2d )
305	ENDIF
306	IF( iom_use( "EXPC" ) ) THEN
307	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
308	DO jk = 1, jpk
309	DO jj = 1, jpj
310	DO ji = 1, jpi
311	zw3d(ji,jj,jk) = ( sinking(ji,jj,jk) + sinking2(ji,jj,jk) ) * zfact * tmask(ji,jj,jk) ! Export of carbon in the water column
312	END DO
313	END DO
314	END DO
315	CALL iom_put( "EXPC" , zw3d )
316	ENDIF
317	IF( iom_use( "EXPFE" ) ) THEN
318	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
319	DO jk = 1, jpk
320	DO jj = 1, jpj
321	DO ji = 1, jpi
322	zw3d(ji,jj,jk) = ( sinkfer(ji,jj,jk) + sinkfer2(ji,jj,jk) ) * zfact * tmask(ji,jj,jk) ! Export of iron
323	END DO
324	END DO
325	END DO
326	CALL iom_put( "EXPFE" , zw3d )
327	ENDIF
328	IF( iom_use( "EXPCAL" ) ) THEN
329	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
330	DO jk = 1, jpk
331	DO jj = 1, jpj
332	DO ji = 1, jpi
333	zw3d(ji,jj,jk) = sinkcal(ji,jj,jk) * zfact * tmask(ji,jj,jk) ! Export of calcite
334	END DO
335	END DO
336	END DO
337	CALL iom_put( "EXPCAL" , zw3d )
338	ENDIF
339	IF( iom_use( "EXPSI" ) ) THEN
340	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
341	DO jk = 1, jpk
342	DO jj = 1, jpj
343	DO ji = 1, jpi
344	zw3d(ji,jj,jk) = sinksil(ji,jj,jk) * zfact * tmask(ji,jj,jk) ! Export of bigenic silica
345	END DO
346	END DO
347	END DO
348	CALL iom_put( "EXPSI" , zw3d )
349	ENDIF
350	IF( iom_use( "tcexp" ) ) CALL iom_put( "tcexp" , t_oce_co2_exp * zfact ) ! molC/s
351	!
352	CALL wrk_dealloc( jpi, jpj, zw2d )
353	CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
354	ENDIF
355	ENDIF
356	!
357	IF(ln_ctl) THEN ! print mean trends (used for debugging)
358	WRITE(charout, FMT="('sink')")
359	CALL prt_ctl_trc_info(charout)
360	CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
361	ENDIF
362	!
363	IF( nn_timing == 1 ) CALL timing_stop('p4z_sink')
364	!
365	END SUBROUTINE p4z_sink
366
367	SUBROUTINE p4z_sink_init
368	!!----------------------------------------------------------------------
369	!! * ROUTINE p4z_sink_init *
370	!!----------------------------------------------------------------------
371	INTEGER :: jk
372
373	ik100 = 10 ! last level where depth less than 100 m
374	DO jk = jpkm1, 1, -1
375	IF( gdept_1d(jk) > 100. ) ik100 = jk - 1
376	END DO
377	IF (lwp) WRITE(numout,*)
378	IF (lwp) WRITE(numout,*) ' Level corresponding to 100m depth ', ik100 + 1
379	IF (lwp) WRITE(numout,*)
380	!
381	t_oce_co2_exp = 0._wp
382	!
383	END SUBROUTINE p4z_sink_init
384
385	SUBROUTINE p4z_sink2( pwsink, psinkflx, jp_tra, kiter )
386	!!---------------------------------------------------------------------
387	!! * ROUTINE p4z_sink2 *
388	!!
389	!! ** Purpose : Compute the sedimentation terms for the various sinking
390	!! particles. The scheme used to compute the trends is based
391	!! on MUSCL.
392	!!
393	!! ** Method : - this ROUTINE compute not exactly the advection but the
394	!! transport term, i.e. div(u*tra).
395	!!---------------------------------------------------------------------
396	!
397	INTEGER , INTENT(in ) :: jp_tra ! tracer index index
398	INTEGER , INTENT(in ) :: kiter ! number of iterations for time-splitting
399	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pwsink ! sinking speed
400	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: psinkflx ! sinking fluxe
401	!!
402	INTEGER :: ji, jj, jk, jn
403	REAL(wp) :: zigma,zew,zign, zflx, zstep
404	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztraz, zakz, zwsink2, ztrb
405	!!---------------------------------------------------------------------
406	!
407	IF( nn_timing == 1 ) CALL timing_start('p4z_sink2')
408	!
409	! Allocate temporary workspace
410	CALL wrk_alloc( jpi, jpj, jpk, ztraz, zakz, zwsink2, ztrb )
411
412	zstep = rfact2 / REAL( kiter, wp ) / 2.
413
414	!$OMP PARALLEL
415	!$OMP DO schedule(static) private(jk, jj, ji)
416	DO jk = 1, jpk
417	DO jj = 1, jpj
418	DO ji = 1, jpi
419	ztraz(ji,jj,jk) = 0.e0
420	zakz (ji,jj,jk) = 0.e0
421	ztrb (ji,jj,jk) = trb(ji,jj,jk,jp_tra)
422	END DO
423	END DO
424	END DO
425	!$OMP END DO NOWAIT
426	!$OMP DO schedule(static) private(jk, jj, ji)
427	DO jk = 1, jpkm1
428	DO jj = 1, jpj
429	DO ji = 1, jpi
430	zwsink2(ji,jj,jk+1) = -pwsink(ji,jj,jk) / rday * tmask(ji,jj,jk+1)
431	END DO
432	END DO
433	END DO
434
435	!$OMP DO schedule(static) private(jj, ji)
436	DO jj = 1, jpj
437	DO ji = 1, jpi
438	zwsink2(ji,jj,1) = 0.e0
439	END DO
440	END DO
441	!$OMP END DO NOWAIT
442
443	! Vertical advective flux
444	DO jn = 1, 2
445	! first guess of the slopes interior values
446	!$OMP DO schedule(static) private(jk,jj,ji)
447	DO jk = 2, jpkm1
448	DO jj = 1, jpj
449	DO ji = 1, jpi
450	ztraz(ji,jj,jk) = ( trb(ji,jj,jk-1,jp_tra) - trb(ji,jj,jk,jp_tra) ) * tmask(ji,jj,jk)
451	END DO
452	END DO
453	END DO
454	!$OMP END DO NOWAIT
455	!$OMP DO schedule(static) private(jj, ji)
456	DO jj = 1, jpj
457	DO ji = 1, jpi
458	ztraz(ji,jj,1 ) = 0.0
459	ztraz(ji,jj,jpk) = 0.0
460	END DO
461	END DO
462
463	! slopes
464	!$OMP DO schedule(static) private(jk, jj, ji, zign)
465	DO jk = 2, jpkm1
466	DO jj = 1,jpj
467	DO ji = 1, jpi
468	zign = 0.25 + SIGN( 0.25, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) )
469	zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign
470	END DO
471	END DO
472	END DO
473
474	! Slopes limitation
475	!$OMP DO schedule(static) private(jk, jj, ji)
476	DO jk = 2, jpkm1
477	DO jj = 1, jpj
478	DO ji = 1, jpi
479	zakz(ji,jj,jk) = SIGN( 1., zakz(ji,jj,jk) ) * &
480	& MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) )
481	END DO
482	END DO
483	END DO
484
485	! vertical advective flux
486	!$OMP DO schedule(static) private(jk, jj, ji, zigma, zew)
487	DO jk = 1, jpkm1
488	DO jj = 1, jpj
489	DO ji = 1, jpi
490	zigma = zwsink2(ji,jj,jk+1) * zstep / e3w_n(ji,jj,jk+1)
491	zew = zwsink2(ji,jj,jk+1)
492	psinkflx(ji,jj,jk+1) = -zew * ( trb(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
493	END DO
494	END DO
495	END DO
496	!
497	! Boundary conditions
498	!$OMP DO schedule(static) private(jj, ji)
499	DO jj = 1, jpj
500	DO ji = 1, jpi
501	psinkflx(ji,jj,1 ) = 0.e0
502	psinkflx(ji,jj,jpk) = 0.e0
503	END DO
504	END DO
505
506	!$OMP DO schedule(static) private(jk, jj, ji, zflx)
507	DO jk=1,jpkm1
508	DO jj = 1,jpj
509	DO ji = 1, jpi
510	zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
511	trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + zflx
512	END DO
513	END DO
514	END DO
515
516	ENDDO
517
518	!$OMP DO schedule(static) private(jk, jj, ji, zflx)
519	DO jk = 1,jpkm1
520	DO jj = 1,jpj
521	DO ji = 1, jpi
522	zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
523	ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx
524	END DO
525	END DO
526	END DO
527
528	!$OMP DO schedule(static) private(jk, jj, ji)
529	DO jk = 1, jpk
530	DO jj = 1, jpj
531	DO ji = 1, jpi
532	trb(ji,jj,jk,jp_tra) = ztrb(ji,jj,jk)
533	psinkflx(ji,jj,jk) = 2. * psinkflx(ji,jj,jk)
534	END DO
535	END DO
536	END DO
537	!$OMP END PARALLEL
538	!
539	CALL wrk_dealloc( jpi, jpj, jpk, ztraz, zakz, zwsink2, ztrb )
540	!
541	IF( nn_timing == 1 ) CALL timing_stop('p4z_sink2')
542	!
543	END SUBROUTINE p4z_sink2
544
545
546	INTEGER FUNCTION p4z_sink_alloc()
547	!!----------------------------------------------------------------------
548	!! * ROUTINE p4z_sink_alloc *
549	!!----------------------------------------------------------------------
550	INTEGER :: ierr(3)
551
552	ierr(:) = 0
553	!
554	ALLOCATE( sinking(jpi,jpj,jpk) , sinking2(jpi,jpj,jpk) , &
555	& sinkcal(jpi,jpj,jpk) , sinksil (jpi,jpj,jpk) , &
556	& sinkfer2(jpi,jpj,jpk) , &
557	& sinkfer(jpi,jpj,jpk) , STAT=ierr(1) )
558	!
559	IF( ln_ligand ) ALLOCATE( sinkfep(jpi,jpj,jpk) , STAT=ierr(2) )
560
561	IF( ln_p5z ) ALLOCATE( sinkingn(jpi,jpj,jpk), sinking2n(jpi,jpj,jpk) , &
562	& sinkingp(jpi,jpj,jpk), sinking2p(jpi,jpj,jpk) , STAT=ierr(3) )
563	!
564	p4z_sink_alloc = MAXVAL( ierr )
565	IF( p4z_sink_alloc /= 0 ) CALL ctl_warn('p4z_sink_alloc : failed to allocate arrays.')
566	!
567	END FUNCTION p4z_sink_alloc
568
569	!!======================================================================
570	END MODULE p4zsink

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

Download in other formats: