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.

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

Context Navigation

source: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.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: 25.0 KB

Line
1	MODULE p4zopt
2	!!======================================================================
3	!! * MODULE p4zopt *
4	!! TOP - PISCES : Compute the light availability in the water column
5	!!======================================================================
6	!! History : 1.0 ! 2004 (O. Aumont) Original code
7	!! 2.0 ! 2007-12 (C. Ethe, G. Madec) F90
8	!! 3.2 ! 2009-04 (C. Ethe, G. Madec) optimisation
9	!! 3.4 ! 2011-06 (O. Aumont, C. Ethe) Improve light availability of nano & diat
10	!!----------------------------------------------------------------------
11	!! p4z_opt : light availability in the water column
12	!!----------------------------------------------------------------------
13	USE trc ! tracer variables
14	USE oce_trc ! tracer-ocean share variables
15	USE sms_pisces ! Source Minus Sink of PISCES
16	USE iom ! I/O manager
17	USE fldread ! time interpolation
18	USE prtctl_trc ! print control for debugging
19
20
21	IMPLICIT NONE
22	PRIVATE
23
24	PUBLIC p4z_opt ! called in p4zbio.F90 module
25	PUBLIC p4z_opt_init ! called in trcsms_pisces.F90 module
26	PUBLIC p4z_opt_alloc
27
28	!! * Shared module variables
29
30	LOGICAL :: ln_varpar !: boolean for variable PAR fraction
31	REAL(wp) :: parlux !: Fraction of shortwave as PAR
32	REAL(wp) :: xparsw !: parlux/3
33	REAL(wp) :: xsi0r !: 1. /rn_si0
34
35	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_par ! structure of input par
36	INTEGER , PARAMETER :: nbtimes = 365 !: maximum number of times record in a file
37	INTEGER :: ntimes_par ! number of time steps in a file
38	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: par_varsw !: PAR fraction of shortwave
39	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ekb, ekg, ekr !: wavelength (Red-Green-Blue)
40
41	INTEGER :: nksrp ! levels below which the light cannot penetrate ( depth larger than 391 m)
42
43	REAL(wp), DIMENSION(3,61) :: xkrgb !: tabulated attenuation coefficients for RGB absorption
44
45	!!----------------------------------------------------------------------
46	!! NEMO/TOP 3.3 , NEMO Consortium (2010)
47	!! $Id: p4zopt.F90 3160 2011-11-20 14:27:18Z cetlod $
48	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
49	!!----------------------------------------------------------------------
50	CONTAINS
51
52	SUBROUTINE p4z_opt( kt, knt )
53	!!---------------------------------------------------------------------
54	!! * ROUTINE p4z_opt *
55	!!
56	!! ** Purpose : Compute the light availability in the water column
57	!! depending on the depth and the chlorophyll concentration
58	!!
59	!! ** Method : - ???
60	!!---------------------------------------------------------------------
61	!
62	INTEGER, INTENT(in) :: kt, knt ! ocean time step
63	!
64	INTEGER :: ji, jj, jk
65	INTEGER :: irgb
66	REAL(wp) :: zchl
67	REAL(wp) :: zc0 , zc1 , zc2, zc3, z1_dep
68	REAL(wp), POINTER, DIMENSION(:,: ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4
69	REAL(wp), POINTER, DIMENSION(:,: ) :: zetmp5
70	REAL(wp), POINTER, DIMENSION(:,: ) :: zqsr100, zqsr_corr
71	REAL(wp), POINTER, DIMENSION(:,:,:) :: zpar, ze0, ze1, ze2, ze3, zchl3d
72	!!---------------------------------------------------------------------
73	!
74	IF( nn_timing == 1 ) CALL timing_start('p4z_opt')
75	!
76	! Allocate temporary workspace
77	CALL wrk_alloc( jpi, jpj, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 )
78	CALL wrk_alloc( jpi, jpj, zqsr100, zqsr_corr )
79	IF( ln_p5z ) CALL wrk_alloc( jpi, jpj, zetmp5 )
80	CALL wrk_alloc( jpi, jpj, jpk, zpar , ze0, ze1, ze2, ze3, zchl3d )
81
82	IF( knt == 1 .AND. ln_varpar ) CALL p4z_opt_sbc( kt )
83
84	! Initialisation of variables used to compute PAR
85	! -----------------------------------------------
86	!$OMP PARALLEL
87	!$OMP DO schedule(static) private(jk,jj,ji)
88	DO jk = 1, jpk
89	DO jj = 1, jpj
90	DO ji = 1, jpi
91	ze1(ji,jj,jk) = 0._wp
92	ze2(ji,jj,jk) = 0._wp
93	ze3(ji,jj,jk) = 0._wp
94	END DO
95	END DO
96	END DO
97	!$OMP END DO NOWAIT
98	!
99	! !* attenuation coef. function of Chlorophyll and wavelength (Red-Green-Blue)
100	! --------------------------------------------------------
101	!$OMP DO schedule(static) private(jk,jj,ji)
102	DO jk = 1, jpk
103	DO jj = 1, jpj
104	DO ji = 1, jpi
105	zchl3d(ji,jj,jk) = trb(ji,jj,jk,jpnch) + trb(ji,jj,jk,jpdch)
106	END DO
107	END DO
108	END DO
109	!$OMP END PARALLEL
110	IF( ln_p5z ) THEN
111	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
112	DO jk = 1, jpk
113	DO jj = 1, jpj
114	DO ji = 1, jpi
115	zchl3d(ji,jj,jk) = zchl3d(ji,jj,jk) + trb(ji,jj,jk,jppch)
116	END DO
117	END DO
118	END DO
119	END IF
120	!
121	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,zchl,irgb)
122	DO jk = 1, jpkm1
123	DO jj = 1, jpj
124	DO ji = 1, jpi
125	zchl = ( zchl3d(ji,jj,jk) + rtrn ) * 1.e6
126	zchl = MIN( 10. , MAX( 0.05, zchl ) )
127	irgb = NINT( 41 + 20.* LOG10( zchl ) + rtrn )
128	!
129	ekb(ji,jj,jk) = xkrgb(1,irgb) * e3t_n(ji,jj,jk)
130	ekg(ji,jj,jk) = xkrgb(2,irgb) * e3t_n(ji,jj,jk)
131	ekr(ji,jj,jk) = xkrgb(3,irgb) * e3t_n(ji,jj,jk)
132	END DO
133	END DO
134	END DO
135	! !* Photosynthetically Available Radiation (PAR)
136	! ! --------------------------------------
137	IF( l_trcdm2dc ) THEN ! diurnal cycle
138	!
139	!$OMP PARALLEL DO schedule(static) private(jj,ji)
140	DO jj = 1, jpj
141	DO ji = 1, jpi
142	zqsr_corr(ji,jj) = qsr_mean(ji,jj) / ( 1. - fr_i(ji,jj) + rtrn )
143	END DO
144	END DO
145	!
146	CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 )
147	!
148	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
149	DO jk = 1, nksrp
150	DO jj = 1, jpj
151	DO ji = 1, jpi
152	etot_ndcy(ji,jj,jk) = ze1(ji,jj,jk) + ze2(ji,jj,jk) + ze3(ji,jj,jk)
153	enano (ji,jj,jk) = 2.1 * ze1(ji,jj,jk) + 0.42 * ze2(ji,jj,jk) + 0.4 * ze3(ji,jj,jk)
154	ediat (ji,jj,jk) = 1.6 * ze1(ji,jj,jk) + 0.69 * ze2(ji,jj,jk) + 0.7 * ze3(ji,jj,jk)
155	END DO
156	END DO
157	END DO
158	IF( ln_p5z ) THEN
159	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
160	DO jk = 1, nksrp
161	DO jj = 1, jpj
162	DO ji = 1, jpi
163	epico (ji,jj,jk) = 2.1 * ze1(ji,jj,jk) + 0.42 * ze2(ji,jj,jk) + 0.4 * ze3(ji,jj,jk)
164	END DO
165	END DO
166	END DO
167	ENDIF
168	!
169	!$OMP PARALLEL DO schedule(static) private(jj,ji)
170	DO jj = 1, jpj
171	DO ji = 1, jpi
172	zqsr_corr(ji,jj) = qsr(ji,jj) / ( 1. - fr_i(ji,jj) + rtrn )
173	END DO
174	END DO
175	!
176	CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 )
177	!
178	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
179	DO jk = 1, nksrp
180	DO jj = 1, jpj
181	DO ji = 1, jpi
182	etot(ji,jj,jk) = ze1(ji,jj,jk) + ze2(ji,jj,jk) + ze3(ji,jj,jk)
183	END DO
184	END DO
185	END DO
186	!
187	ELSE
188	!
189	!$OMP PARALLEL DO schedule(static) private(jj,ji)
190	DO jj = 1, jpj
191	DO ji = 1, jpi
192	zqsr_corr(ji,jj) = qsr(ji,jj) / ( 1. - fr_i(ji,jj) + rtrn )
193	END DO
194	END DO
195	!
196	CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 )
197	!
198	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
199	DO jk = 1, nksrp
200	DO jj = 1, jpj
201	DO ji = 1, jpi
202	etot (ji,jj,jk) = ze1(ji,jj,jk) + ze2(ji,jj,jk) + ze3(ji,jj,jk)
203	enano(ji,jj,jk) = 2.1 * ze1(ji,jj,jk) + 0.42 * ze2(ji,jj,jk) + 0.4 * ze3(ji,jj,jk)
204	ediat(ji,jj,jk) = 1.6 * ze1(ji,jj,jk) + 0.69 * ze2(ji,jj,jk) + 0.7 * ze3(ji,jj,jk)
205	END DO
206	END DO
207	END DO
208	IF( ln_p5z ) THEN
209	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
210	DO jk = 1, nksrp
211	DO jj = 1, jpj
212	DO ji = 1, jpi
213	epico(ji,jj,jk) = 2.1 * ze1(ji,jj,jk) + 0.42 * ze2(ji,jj,jk) + 0.4 * ze3(ji,jj,jk)
214	END DO
215	END DO
216	END DO
217	ENDIF
218	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
219	DO jk = 1, jpk
220	DO jj = 1, jpj
221	DO ji = 1, jpi
222	etot_ndcy(ji,jj,jk) = etot(ji,jj,jk)
223	END DO
224	END DO
225	END DO
226	ENDIF
227
228
229	IF( ln_qsr_bio ) THEN !* heat flux accros w-level (used in the dynamics)
230	! ! ------------------------
231	CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3, pe0=ze0 )
232	!
233	!$OMP PARALLEL
234	!$OMP DO schedule(static) private(jj,ji)
235	DO jj = 1, jpj
236	DO ji = 1, jpi
237	etot3(ji,jj,1) = qsr(ji,jj) * tmask(ji,jj,1)
238	END DO
239	END DO
240	!$OMP DO schedule(static) private(jk,jj,ji)
241	DO jk = 2, nksrp + 1
242	DO jj = 1, jpj
243	DO ji = 1, jpi
244	etot3(ji,jj,jk) = ( ze0(ji,jj,jk) + ze1(ji,jj,jk) + ze2(ji,jj,jk) + ze3(ji,jj,jk) ) * tmask(ji,jj,jk)
245	END DO
246	END DO
247	END DO
248	!$OMP END PARALLEL
249	! ! ------------------------
250	ENDIF
251	! !* Euphotic depth and level
252	! ------------------------
253	!$OMP PARALLEL
254	!$OMP DO schedule(static) private(jj,ji)
255	DO jj = 1, jpj
256	DO ji = 1, jpi
257	neln(ji,jj) = 1
258	heup (ji,jj) = gdepw_n(ji,jj,2)
259	heup_01(ji,jj) = gdepw_n(ji,jj,2)
260	END DO
261	END DO
262
263	DO jk = 2, nksrp
264	!$OMP DO schedule(static) private(jj,ji)
265	DO jj = 1, jpj
266	DO ji = 1, jpi
267	IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= zqsr100(ji,jj) ) THEN
268	neln(ji,jj) = jk+1 ! Euphotic level : 1rst T-level strictly below Euphotic layer
269	! ! nb: ensure the compatibility with nmld_trc definition in trd_mld_trc_zint
270	heup(ji,jj) = gdepw_n(ji,jj,jk+1) ! Euphotic layer depth
271	ENDIF
272	IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.50 ) THEN
273	heup_01(ji,jj) = gdepw_n(ji,jj,jk+1) ! Euphotic layer depth (light level definition)
274	ENDIF
275	END DO
276	END DO
277	END DO
278	!
279	!$OMP DO schedule(static) private(jj,ji)
280	DO jj = 1, jpj
281	DO ji = 1, jpi
282	heup (ji,jj) = MIN( 300., heup (ji,jj) )
283	heup_01(ji,jj) = MIN( 300., heup_01(ji,jj) )
284	! !* mean light over the mixed layer
285	zdepmoy(ji,jj) = 0.e0 ! -------------------------------
286	zetmp1 (ji,jj) = 0.e0
287	zetmp2 (ji,jj) = 0.e0
288	zetmp3 (ji,jj) = 0.e0
289	zetmp4 (ji,jj) = 0.e0
290	END DO
291	END DO
292
293	DO jk = 1, nksrp
294	!$OMP DO schedule(static) private(jj,ji)
295	DO jj = 1, jpj
296	DO ji = 1, jpi
297	IF( gdepw_n(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
298	zetmp1 (ji,jj) = zetmp1 (ji,jj) + etot (ji,jj,jk) * e3t_n(ji,jj,jk) ! remineralisation
299	zetmp2 (ji,jj) = zetmp2 (ji,jj) + etot_ndcy(ji,jj,jk) * e3t_n(ji,jj,jk) ! production
300	zetmp3 (ji,jj) = zetmp3 (ji,jj) + enano (ji,jj,jk) * e3t_n(ji,jj,jk) ! production
301	zetmp4 (ji,jj) = zetmp4 (ji,jj) + ediat (ji,jj,jk) * e3t_n(ji,jj,jk) ! production
302	zdepmoy(ji,jj) = zdepmoy(ji,jj) + e3t_n(ji,jj,jk)
303	ENDIF
304	END DO
305	END DO
306	END DO
307	!
308	!$OMP DO schedule(static) private(jk,jj,ji)
309	DO jk = 1, jpk
310	DO jj = 1, jpj
311	DO ji = 1, jpi
312	emoy(ji,jj,jk) = etot(ji,jj,jk) ! remineralisation
313	zpar(ji,jj,jk) = etot_ndcy(ji,jj,jk) ! diagnostic : PAR with no diurnal cycle
314	END DO
315	END DO
316	END DO
317	!
318	!$OMP DO schedule(static) private(jk,jj,ji,z1_dep)
319	DO jk = 1, nksrp
320	DO jj = 1, jpj
321	DO ji = 1, jpi
322	IF( gdepw_n(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
323	z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
324	emoy (ji,jj,jk) = zetmp1(ji,jj) * z1_dep
325	zpar (ji,jj,jk) = zetmp2(ji,jj) * z1_dep
326	enano(ji,jj,jk) = zetmp3(ji,jj) * z1_dep
327	ediat(ji,jj,jk) = zetmp4(ji,jj) * z1_dep
328	ENDIF
329	END DO
330	END DO
331	END DO
332	!$OMP END PARALLEL
333	!
334	IF( ln_p5z ) THEN
335	!$OMP PARALLEL
336	!$OMP DO schedule(static) private(jj,ji)
337	DO jj = 1, jpj
338	DO ji = 1, jpi
339	zetmp5 (ji,jj) = 0.e0
340	END DO
341	END DO
342	DO jk = 1, nksrp
343	!$OMP DO schedule(static) private(jj,ji,z1_dep)
344	DO jj = 1, jpj
345	DO ji = 1, jpi
346	IF( gdepw_n(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
347	z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
348	zetmp5(ji,jj) = zetmp5 (ji,jj) + epico(ji,jj,jk) * e3t_n(ji,jj,jk) ! production
349	epico(ji,jj,jk) = zetmp5(ji,jj) * z1_dep
350	ENDIF
351	END DO
352	END DO
353	END DO
354	!$OMP END PARALLEL
355	ENDIF
356	IF( lk_iomput ) THEN
357	IF( knt == nrdttrc ) THEN
358	IF( iom_use( "Heup" ) ) CALL iom_put( "Heup" , heup(:,: ) * tmask(:,:,1) ) ! euphotic layer deptht
359	IF( iom_use( "PARDM" ) ) CALL iom_put( "PARDM", zpar(:,:,:) * tmask(:,:,:) ) ! Photosynthetically Available Radiation
360	IF( iom_use( "PAR" ) ) CALL iom_put( "PAR" , emoy(:,:,:) * tmask(:,:,:) ) ! Photosynthetically Available Radiation
361	ENDIF
362	ENDIF
363	!
364	CALL wrk_dealloc( jpi, jpj, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 )
365	CALL wrk_dealloc( jpi, jpj, zqsr100, zqsr_corr )
366	IF( ln_p5z ) CALL wrk_dealloc( jpi, jpj, zetmp5 )
367	CALL wrk_dealloc( jpi, jpj, jpk, zpar , ze0, ze1, ze2, ze3, zchl3d )
368	!
369	IF( nn_timing == 1 ) CALL timing_stop('p4z_opt')
370	!
371	END SUBROUTINE p4z_opt
372
373	SUBROUTINE p4z_opt_par( kt, pqsr, pe1, pe2, pe3, pe0, pqsr100 )
374	!!----------------------------------------------------------------------
375	!! * routine p4z_opt_par *
376	!!
377	!! ** purpose : compute PAR of each wavelength (Red-Green-Blue)
378	!! for a given shortwave radiation
379	!!
380	!!----------------------------------------------------------------------
381	!! * arguments
382	INTEGER, INTENT(in) :: kt ! ocean time-step
383	REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: pqsr ! shortwave
384	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pe1 , pe2 , pe3 ! PAR ( R-G-B)
385	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL :: pe0
386	REAL(wp), DIMENSION(jpi,jpj) , INTENT(out) , OPTIONAL :: pqsr100
387	!! * local variables
388	INTEGER :: ji, jj, jk ! dummy loop indices
389	REAL(wp), DIMENSION(jpi,jpj) :: zqsr ! shortwave
390	!!----------------------------------------------------------------------
391
392	! Real shortwave
393	IF( ln_varpar ) THEN
394	!$OMP PARALLEL DO schedule(static) private(jj,ji)
395	DO jj = 1, jpj
396	DO ji = 1, jpi
397	zqsr(ji,jj) = par_varsw(ji,jj) * pqsr(ji,jj)
398	END DO
399	END DO
400	ELSE
401	!$OMP PARALLEL DO schedule(static) private(jj,ji)
402	DO jj = 1, jpj
403	DO ji = 1, jpi
404	zqsr(ji,jj) = xparsw * pqsr(ji,jj)
405	END DO
406	END DO
407	ENDIF
408
409	! Light at the euphotic depth
410	IF( PRESENT( pqsr100 ) ) THEN
411	!$OMP PARALLEL DO schedule(static) private(jj,ji)
412	DO jj = 1, jpj
413	DO ji = 1, jpi
414	pqsr100(ji,jj) = 0.01 * 3. * zqsr(ji,jj)
415	END DO
416	END DO
417	ENDIF
418
419	IF( PRESENT( pe0 ) ) THEN ! W-level
420	!
421	!$OMP PARALLEL
422	!$OMP DO schedule(static) private(jj,ji)
423	DO jj = 1, jpj
424	DO ji = 1, jpi
425	pe0(ji,jj,1) = pqsr(ji,jj) - 3. * zqsr(ji,jj) ! ( 1 - 3 * alpha ) * q
426	pe1(ji,jj,1) = zqsr(ji,jj)
427	pe2(ji,jj,1) = zqsr(ji,jj)
428	pe3(ji,jj,1) = zqsr(ji,jj)
429	END DO
430	END DO
431	!
432	DO jk = 2, nksrp + 1
433	!$OMP DO schedule(static) private(jj,ji)
434	DO jj = 1, jpj
435	DO ji = 1, jpi
436	pe0(ji,jj,jk) = pe0(ji,jj,jk-1) * EXP( -e3t_n(ji,jj,jk-1) * xsi0r )
437	pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -ekb(ji,jj,jk-1 ) )
438	pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -ekg(ji,jj,jk-1 ) )
439	pe3(ji,jj,jk) = pe3(ji,jj,jk-1) * EXP( -ekr(ji,jj,jk-1 ) )
440	END DO
441	!
442	END DO
443	!
444	END DO
445	!$OMP END PARALLEL
446	!
447	ELSE ! T- level
448	!
449	!$OMP PARALLEL
450	!$OMP DO schedule(static) private(jj,ji)
451	DO jj = 1, jpj
452	DO ji = 1, jpi
453	pe1(ji,jj,1) = zqsr(ji,jj) * EXP( -0.5 * ekb(ji,jj,1) )
454	pe2(ji,jj,1) = zqsr(ji,jj) * EXP( -0.5 * ekg(ji,jj,1) )
455	pe3(ji,jj,1) = zqsr(ji,jj) * EXP( -0.5 * ekr(ji,jj,1) )
456	END DO
457	END DO
458	!
459	DO jk = 2, nksrp
460	!$OMP DO schedule(static) private(jj,ji)
461	DO jj = 1, jpj
462	DO ji = 1, jpi
463	pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -0.5 * ( ekb(ji,jj,jk-1) + ekb(ji,jj,jk) ) )
464	pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -0.5 * ( ekg(ji,jj,jk-1) + ekg(ji,jj,jk) ) )
465	pe3(ji,jj,jk) = pe3(ji,jj,jk-1) * EXP( -0.5 * ( ekr(ji,jj,jk-1) + ekr(ji,jj,jk) ) )
466	END DO
467	END DO
468	END DO
469	!$OMP END PARALLEL
470	!
471	ENDIF
472	!
473	END SUBROUTINE p4z_opt_par
474
475
476	SUBROUTINE p4z_opt_sbc( kt )
477	!!----------------------------------------------------------------------
478	!! * routine p4z_opt_sbc *
479	!!
480	!! ** purpose : read and interpolate the variable PAR fraction
481	!! of shortwave radiation
482	!!
483	!! ** method : read the files and interpolate the appropriate variables
484	!!
485	!! ** input : external netcdf files
486	!!
487	!!----------------------------------------------------------------------
488	!! * arguments
489	INTEGER , INTENT(in) :: kt ! ocean time step
490
491	!! * local declarations
492	INTEGER :: ji,jj
493	REAL(wp) :: zcoef
494	!!---------------------------------------------------------------------
495	!
496	IF( nn_timing == 1 ) CALL timing_start('p4z_optsbc')
497	!
498	! Compute par_varsw at nit000 or only if there is more than 1 time record in par coefficient file
499	IF( ln_varpar ) THEN
500	IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_par > 1 ) ) THEN
501	CALL fld_read( kt, 1, sf_par )
502	par_varsw(:,:) = ( sf_par(1)%fnow(:,:,1) ) / 3.0
503	ENDIF
504	ENDIF
505	!
506	IF( nn_timing == 1 ) CALL timing_stop('p4z_optsbc')
507	!
508	END SUBROUTINE p4z_opt_sbc
509
510	SUBROUTINE p4z_opt_init
511	!!----------------------------------------------------------------------
512	!! * ROUTINE p4z_opt_init *
513	!!
514	!! ** Purpose : Initialization of tabulated attenuation coef
515	!! and of the percentage of PAR in Shortwave
516	!!
517	!! ** Input : external ascii and netcdf files
518	!!----------------------------------------------------------------------
519	!
520	INTEGER :: numpar
521	INTEGER :: ierr
522	INTEGER :: ios ! Local integer output status for namelist read
523	INTEGER :: ji, jj, jk ! dummy loop indices
524	REAL(wp), DIMENSION(nbtimes) :: zsteps ! times records
525	!
526	CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files
527	TYPE(FLD_N) :: sn_par ! informations about the fields to be read
528	!
529	NAMELIST/nampisopt/cn_dir, sn_par, ln_varpar, parlux
530
531	!!----------------------------------------------------------------------
532
533	IF( nn_timing == 1 ) CALL timing_start('p4z_opt_init')
534
535	REWIND( numnatp_ref ) ! Namelist nampisopt in reference namelist : Pisces attenuation coef. and PAR
536	READ ( numnatp_ref, nampisopt, IOSTAT = ios, ERR = 901)
537	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisopt in reference namelist', lwp )
538
539	REWIND( numnatp_cfg ) ! Namelist nampisopt in configuration namelist : Pisces attenuation coef. and PAR
540	READ ( numnatp_cfg, nampisopt, IOSTAT = ios, ERR = 902 )
541	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisopt in configuration namelist', lwp )
542	IF(lwm) WRITE ( numonp, nampisopt )
543
544	IF(lwp) THEN
545	WRITE(numout,*) ' '
546	WRITE(numout,*) ' namelist : nampisopt '
547	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ '
548	WRITE(numout,*) ' PAR as a variable fraction of SW ln_varpar = ', ln_varpar
549	WRITE(numout,*) ' Default value for the PAR fraction parlux = ', parlux
550	ENDIF
551	!
552	xparsw = parlux / 3.0
553	xsi0r = 1.e0 / rn_si0
554	!
555	! Variable PAR at the surface of the ocean
556	! ----------------------------------------
557	IF( ln_varpar ) THEN
558	IF(lwp) WRITE(numout,*) ' initialize variable par fraction '
559	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
560	!
561	ALLOCATE( par_varsw(jpi,jpj) )
562	!
563	ALLOCATE( sf_par(1), STAT=ierr ) !* allocate and fill sf_sst (forcing structure) with sn_sst
564	IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'p4z_opt_init: unable to allocate sf_par structure' )
565	!
566	CALL fld_fill( sf_par, (/ sn_par /), cn_dir, 'p4z_opt_init', 'Variable PAR fraction ', 'nampisopt' )
567	ALLOCATE( sf_par(1)%fnow(jpi,jpj,1) )
568	IF( sn_par%ln_tint ) ALLOCATE( sf_par(1)%fdta(jpi,jpj,1,2) )
569
570	CALL iom_open ( TRIM( sn_par%clname ) , numpar )
571	CALL iom_gettime( numpar, zsteps, kntime=ntimes_par) ! get number of record in file
572	ENDIF
573	!
574	CALL trc_oce_rgb( xkrgb ) ! tabulated attenuation coefficients
575	nksrp = trc_oce_ext_lev( r_si2, 0.33e2 ) ! max level of light extinction (Blue Chl=0.01)
576	!
577	IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksrp, ' ref depth = ', gdepw_1d(nksrp+1), ' m'
578	!
579	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
580	DO jk = 1, jpk
581	DO jj = 1, jpj
582	DO ji = 1, jpi
583	ekr (ji,jj,jk) = 0._wp
584	ekb (ji,jj,jk) = 0._wp
585	ekg (ji,jj,jk) = 0._wp
586	etot (ji,jj,jk) = 0._wp
587	etot_ndcy(ji,jj,jk) = 0._wp
588	enano (ji,jj,jk) = 0._wp
589	ediat (ji,jj,jk) = 0._wp
590	END DO
591	END DO
592	END DO
593	IF( ln_qsr_bio ) THEN
594	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
595	DO jk = 1, jpk
596	DO jj = 1, jpj
597	DO ji = 1, jpi
598	etot3 (ji,jj,jk) = 0._wp
599	END DO
600	END DO
601	END DO
602	END IF
603
604	IF( ln_p5z ) THEN
605	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
606	DO jk = 1, jpk
607	DO jj = 1, jpj
608	DO ji = 1, jpi
609	epico (ji,jj,jk) = 0._wp
610	END DO
611	END DO
612	END DO
613	END IF
614	!
615	IF( nn_timing == 1 ) CALL timing_stop('p4z_opt_init')
616	!
617	END SUBROUTINE p4z_opt_init
618
619
620	INTEGER FUNCTION p4z_opt_alloc()
621	!!----------------------------------------------------------------------
622	!! * ROUTINE p4z_opt_alloc *
623	!!----------------------------------------------------------------------
624	!
625	ALLOCATE( ekb(jpi,jpj,jpk), ekr(jpi,jpj,jpk), &
626	ekg(jpi,jpj,jpk), STAT= p4z_opt_alloc )
627	!
628	IF( p4z_opt_alloc /= 0 ) CALL ctl_warn('p4z_opt_alloc : failed to allocate arrays.')
629	!
630	END FUNCTION p4z_opt_alloc
631
632	!!======================================================================
633	END MODULE p4zopt

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

Download in other formats: