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.

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

Context Navigation

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

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago
update trunk with OpenMP parallelization
Property svn:keywords set to `Id`
File size: 29.1 KB

Line
1	MODULE p4zsbc
2	!!======================================================================
3	!! * MODULE p4sbc *
4	!! TOP : PISCES surface boundary conditions of external inputs of nutrients
5	!!======================================================================
6	!! History : 3.5 ! 2012-07 (O. Aumont, C. Ethe) Original code
7	!!----------------------------------------------------------------------
8	!! p4z_sbc : Read and interpolate time-varying nutrients fluxes
9	!! p4z_sbc_init : Initialization of p4z_sbc
10	!!----------------------------------------------------------------------
11	USE oce_trc ! shared variables between ocean and passive tracers
12	USE trc ! passive tracers common variables
13	USE sms_pisces ! PISCES Source Minus Sink variables
14	USE iom ! I/O manager
15	USE fldread ! time interpolation
16
17	IMPLICIT NONE
18	PRIVATE
19
20	PUBLIC p4z_sbc
21	PUBLIC p4z_sbc_init
22
23	LOGICAL , PUBLIC :: ln_dust !: boolean for dust input from the atmosphere
24	LOGICAL , PUBLIC :: ln_solub !: boolean for variable solubility of atmospheric iron
25	LOGICAL , PUBLIC :: ln_river !: boolean for river input of nutrients
26	LOGICAL , PUBLIC :: ln_ndepo !: boolean for atmospheric deposition of N
27	LOGICAL , PUBLIC :: ln_ironsed !: boolean for Fe input from sediments
28	LOGICAL , PUBLIC :: ln_hydrofe !: boolean for Fe input from hydrothermal vents
29	LOGICAL , PUBLIC :: ln_ironice !: boolean for Fe input from sea ice
30	REAL(wp), PUBLIC :: sedfeinput !: Coastal release of Iron
31	REAL(wp), PUBLIC :: dustsolub !: Solubility of the dust
32	REAL(wp), PUBLIC :: mfrac !: Mineral Content of the dust
33	REAL(wp), PUBLIC :: icefeinput !: Iron concentration in sea ice
34	REAL(wp), PUBLIC :: wdust !: Sinking speed of the dust
35	REAL(wp), PUBLIC :: nitrfix !: Nitrogen fixation rate
36	REAL(wp), PUBLIC :: diazolight !: Nitrogen fixation sensitivty to light
37	REAL(wp), PUBLIC :: concfediaz !: Fe half-saturation Cste for diazotrophs
38	REAL(wp) :: hratio !: Fe:3He ratio assumed for vent iron supply
39	REAL(wp), PUBLIC :: fep_rats !: Fep/Fer ratio from sed sources
40	REAL(wp), PUBLIC :: fep_rath !: Fep/Fer ratio from hydro sources
41	REAL(wp), PUBLIC :: lgw_rath !: Weak ligand ratio from hydro sources
42
43
44	LOGICAL , PUBLIC :: ll_sbc
45
46	LOGICAL :: ll_solub
47
48	INTEGER , PARAMETER :: jpriv = 7 !: Maximum number of river input fields
49	INTEGER , PARAMETER :: jr_dic = 1 !: index of dissolved inorganic carbon
50	INTEGER , PARAMETER :: jr_doc = 2 !: index of dissolved organic carbon
51	INTEGER , PARAMETER :: jr_din = 3 !: index of dissolved inorganic nitrogen
52	INTEGER , PARAMETER :: jr_don = 4 !: index of dissolved organic nitrogen
53	INTEGER , PARAMETER :: jr_dip = 5 !: index of dissolved inorganic phosporus
54	INTEGER , PARAMETER :: jr_dop = 6 !: index of dissolved organic phosphorus
55	INTEGER , PARAMETER :: jr_dsi = 7 !: index of dissolved silicate
56
57
58	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_dust ! structure of input dust
59	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_solub ! structure of input dust
60	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_river ! structure of input riverdic
61	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_ndepo ! structure of input nitrogen deposition
62	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_ironsed ! structure of input iron from sediment
63	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_hydrofe ! structure of input iron from hydrothermal vents
64
65	INTEGER , PARAMETER :: nbtimes = 365 !: maximum number of times record in a file
66	INTEGER :: ntimes_dust, ntimes_riv, ntimes_ndep ! number of time steps in a file
67	INTEGER :: ntimes_solub, ntimes_hydro ! number of time steps in a file
68
69	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dust, solub !: dust fields
70	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rivdic, rivalk !: river input fields
71	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rivdin, rivdip !: river input fields
72	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rivdon, rivdop !: river input fields
73	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rivdoc !: river input fields
74	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rivdsi !: river input fields
75	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: nitdep !: atmospheric N deposition
76	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ironsed !: Coastal supply of iron
77	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hydrofe !: Hydrothermal vent supply of iron
78
79	REAL(wp), PUBLIC :: sumdepsi, rivalkinput, rivdicinput, nitdepinput
80	REAL(wp), PUBLIC :: rivdininput, rivdipinput, rivdsiinput
81
82	!! * Substitutions
83	# include "vectopt_loop_substitute.h90"
84	!!----------------------------------------------------------------------
85	!! NEMO/TOP 3.3 , NEMO Consortium (2010)
86	!! $Id$
87	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
88	!!----------------------------------------------------------------------
89	CONTAINS
90
91	SUBROUTINE p4z_sbc( kt )
92	!!----------------------------------------------------------------------
93	!! * routine p4z_sbc *
94	!!
95	!! ** purpose : read and interpolate the external sources of nutrients
96	!!
97	!! ** method : read the files and interpolate the appropriate variables
98	!!
99	!! ** input : external netcdf files
100	!!
101	!!----------------------------------------------------------------------
102	!! * arguments
103	INTEGER, INTENT( in ) :: kt ! ocean time step
104
105	!! * local declarations
106	INTEGER :: ji,jj
107	REAL(wp) :: zcoef, zyyss
108	!!---------------------------------------------------------------------
109	!
110	IF( nn_timing == 1 ) CALL timing_start('p4z_sbc')
111
112	!
113	! Compute dust at nit000 or only if there is more than 1 time record in dust file
114	IF( ln_dust ) THEN
115	IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_dust > 1 ) ) THEN
116	CALL fld_read( kt, 1, sf_dust )
117	IF( nn_ice_tr == -1 .AND. .NOT. ln_ironice ) THEN
118	!$OMP PARALLEL DO schedule(static) private(jj, ji)
119	DO jj = 1, jpj
120	DO ji = 1, jpi
121	dust(ji,jj) = sf_dust(1)%fnow(ji,jj,1)
122	END DO
123	END DO
124	ELSE
125	!$OMP PARALLEL DO schedule(static) private(jj, ji)
126	DO jj = 1, jpj
127	DO ji = 1, jpi
128	dust(ji,jj) = sf_dust(1)%fnow(ji,jj,1) * ( 1.0 - fr_i(ji,jj) )
129	END DO
130	END DO
131	ENDIF
132	ENDIF
133	ENDIF
134
135	IF( ll_solub ) THEN
136	IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_solub > 1 ) ) THEN
137	CALL fld_read( kt, 1, sf_solub )
138	!$OMP PARALLEL DO schedule(static) private(jj, ji)
139	DO jj = 1, jpj
140	DO ji = 1, jpi
141	solub(ji,jj) = sf_solub(1)%fnow(ji,jj,1)
142	END DO
143	END DO
144	ENDIF
145	ENDIF
146
147	! N/P and Si releases due to coastal rivers
148	! Compute river at nit000 or only if there is more than 1 time record in river file
149	! -----------------------------------------
150	IF( ln_river ) THEN
151	IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_riv > 1 ) ) THEN
152	CALL fld_read( kt, 1, sf_river )
153	IF( ln_p4z ) THEN
154	!$OMP PARALLEL DO schedule(static) private(jj, ji, zcoef)
155	DO jj = 1, jpj
156	DO ji = 1, jpi
157	zcoef = ryyss * e1e2t(ji,jj) * h_rnf(ji,jj)
158	rivalk(ji,jj) = sf_river(jr_dic)%fnow(ji,jj,1) &
159	& * 1.E3 / ( 12. * zcoef + rtrn )
160	rivdic(ji,jj) = ( sf_river(jr_dic)%fnow(ji,jj,1) + sf_river(jr_doc)%fnow(ji,jj,1) ) &
161	& * 1.E3 / ( 12. * zcoef + rtrn )
162	rivdin(ji,jj) = ( sf_river(jr_din)%fnow(ji,jj,1) + sf_river(jr_don)%fnow(ji,jj,1) ) &
163	& * 1.E3 / rno3 / ( 14. * zcoef + rtrn )
164	rivdip(ji,jj) = ( sf_river(jr_dip)%fnow(ji,jj,1) + sf_river(jr_dop)%fnow(ji,jj,1) ) &
165	& * 1.E3 / po4r / ( 31. * zcoef + rtrn )
166	rivdsi(ji,jj) = sf_river(jr_dsi)%fnow(ji,jj,1) &
167	& * 1.E3 / ( 28.1 * zcoef + rtrn )
168	END DO
169	END DO
170	ELSE ! ln_p5z
171	!$OMP PARALLEL DO schedule(static) private(jj, ji, zcoef)
172	DO jj = 1, jpj
173	DO ji = 1, jpi
174	zcoef = ryyss * e1e2t(ji,jj) * h_rnf(ji,jj)
175	rivalk(ji,jj) = sf_river(jr_dic)%fnow(ji,jj,1) &
176	& * 1.E3 / ( 12. * zcoef + rtrn )
177	rivdic(ji,jj) = ( sf_river(jr_dic)%fnow(ji,jj,1) ) &
178	& * 1.E3 / ( 12. * zcoef + rtrn ) * tmask(ji,jj,1)
179	rivdin(ji,jj) = ( sf_river(jr_din)%fnow(ji,jj,1) ) &
180	& * 1.E3 / rno3 / ( 14. * zcoef + rtrn ) * tmask(ji,jj,1)
181	rivdip(ji,jj) = ( sf_river(jr_dip)%fnow(ji,jj,1) ) &
182	& * 1.E3 / po4r / ( 31. * zcoef + rtrn ) * tmask(ji,jj,1)
183	rivdoc(ji,jj) = ( sf_river(jr_doc)%fnow(ji,jj,1) ) &
184	& * 1.E3 / ( 12. * zcoef + rtrn ) * tmask(ji,jj,1)
185	rivdon(ji,jj) = ( sf_river(jr_don)%fnow(ji,jj,1) ) &
186	& * 1.E3 / rno3 / ( 14. * zcoef + rtrn ) * tmask(ji,jj,1)
187	rivdop(ji,jj) = ( sf_river(jr_dop)%fnow(ji,jj,1) ) &
188	& * 1.E3 / po4r / ( 31. * zcoef + rtrn ) * tmask(ji,jj,1)
189	END DO
190	END DO
191	ENDIF
192	ENDIF
193	ENDIF
194
195	! Compute N deposition at nit000 or only if there is more than 1 time record in N deposition file
196	IF( ln_ndepo ) THEN
197	IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_ndep > 1 ) ) THEN
198	zcoef = rno3 * 14E6 * ryyss
199	CALL fld_read( kt, 1, sf_ndepo )
200	!$OMP PARALLEL DO schedule(static) private(jj, ji)
201	DO jj = 1, jpj
202	DO ji = 1, jpi
203	nitdep(ji,jj) = sf_ndepo(1)%fnow(ji,jj,1) / zcoef / e3t_n(ji,jj,1)
204	END DO
205	END DO
206	ENDIF
207	IF( .NOT.ln_linssh ) THEN
208	zcoef = rno3 * 14E6 * ryyss
209	!$OMP PARALLEL DO schedule(static) private(jj, ji)
210	DO jj = 1, jpj
211	DO ji = 1, jpi
212	nitdep(ji,jj) = sf_ndepo(1)%fnow(ji,jj,1) / zcoef / e3t_n(ji,jj,1)
213	END DO
214	END DO
215	ENDIF
216	ENDIF
217	!
218	IF( nn_timing == 1 ) CALL timing_stop('p4z_sbc')
219	!
220	END SUBROUTINE p4z_sbc
221
222	SUBROUTINE p4z_sbc_init
223
224	!!----------------------------------------------------------------------
225	!! * routine p4z_sbc_init *
226	!!
227	!! ** purpose : initialization of the external sources of nutrients
228	!!
229	!! ** method : read the files and compute the budget
230	!! called at the first timestep (nittrc000)
231	!!
232	!! ** input : external netcdf files
233	!!
234	!!----------------------------------------------------------------------
235	!
236	INTEGER :: ji, jj, jk, jm, ifpr
237	INTEGER :: ii0, ii1, ij0, ij1
238	INTEGER :: numdust, numsolub, numriv, numiron, numdepo, numhydro
239	INTEGER :: ierr, ierr1, ierr2, ierr3
240	INTEGER :: ios ! Local integer output status for namelist read
241	INTEGER :: ik50 ! last level where depth less than 50 m
242	INTEGER :: isrow ! index for ORCA1 starting row
243	REAL(wp) :: zexpide, zdenitide, zmaskt
244	REAL(wp) :: ztimes_dust, ztimes_riv, ztimes_ndep
245	REAL(wp), DIMENSION(nbtimes) :: zsteps ! times records
246	REAL(wp), DIMENSION(:), ALLOCATABLE :: rivinput
247	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zdust, zndepo, zriver, zcmask
248	!
249	CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files
250	TYPE(FLD_N), DIMENSION(jpriv) :: slf_river ! array of namelist informations on the fields to read
251	TYPE(FLD_N) :: sn_dust, sn_solub, sn_ndepo, sn_ironsed, sn_hydrofe ! informations about the fields to be read
252	TYPE(FLD_N) :: sn_riverdoc, sn_riverdic, sn_riverdsi ! informations about the fields to be read
253	TYPE(FLD_N) :: sn_riverdin, sn_riverdon, sn_riverdip, sn_riverdop
254	!
255	NAMELIST/nampissbc/cn_dir, sn_dust, sn_solub, sn_riverdic, sn_riverdoc, sn_riverdin, sn_riverdon, &
256	& sn_riverdip, sn_riverdop, sn_riverdsi, sn_ndepo, sn_ironsed, sn_hydrofe, &
257	& ln_dust, ln_solub, ln_river, ln_ndepo, ln_ironsed, ln_ironice, ln_hydrofe, &
258	& sedfeinput, dustsolub, icefeinput, wdust, mfrac, nitrfix, diazolight, concfediaz, &
259	& hratio, fep_rats, fep_rath, lgw_rath
260	!!----------------------------------------------------------------------
261	!
262	IF( nn_timing == 1 ) CALL timing_start('p4z_sbc_init')
263	!
264	! !* set file information
265	REWIND( numnatp_ref ) ! Namelist nampissbc in reference namelist : Pisces external sources of nutrients
266	READ ( numnatp_ref, nampissbc, IOSTAT = ios, ERR = 901)
267	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampissbc in reference namelist', lwp )
268
269	REWIND( numnatp_cfg ) ! Namelist nampissbc in configuration namelist : Pisces external sources of nutrients
270	READ ( numnatp_cfg, nampissbc, IOSTAT = ios, ERR = 902 )
271	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampissbc in configuration namelist', lwp )
272	IF(lwm) WRITE ( numonp, nampissbc )
273
274	IF ( ( nn_ice_tr >= 0 ) .AND. ln_ironice ) THEN
275	IF(lwp) THEN
276	WRITE(numout,*) ' ln_ironice incompatible with nn_ice_tr = ', nn_ice_tr
277	WRITE(numout,*) ' Specify your sea ice iron concentration in nampisice instead '
278	WRITE(numout,*) ' ln_ironice is forced to .FALSE. '
279	ln_ironice = .FALSE.
280	ENDIF
281	ENDIF
282
283	IF(lwp) THEN
284	WRITE(numout,*) ' '
285	WRITE(numout,*) ' namelist : nampissbc '
286	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ '
287	WRITE(numout,*) ' dust input from the atmosphere ln_dust = ', ln_dust
288	WRITE(numout,*) ' Variable solubility of iron input ln_solub = ', ln_solub
289	WRITE(numout,*) ' river input of nutrients ln_river = ', ln_river
290	WRITE(numout,*) ' atmospheric deposition of n ln_ndepo = ', ln_ndepo
291	WRITE(numout,*) ' Fe input from sediments ln_ironsed = ', ln_ironsed
292	WRITE(numout,*) ' Fe input from seaice ln_ironice = ', ln_ironice
293	WRITE(numout,*) ' fe input from hydrothermal vents ln_hydrofe = ', ln_hydrofe
294	WRITE(numout,*) ' coastal release of iron sedfeinput = ', sedfeinput
295	WRITE(numout,*) ' solubility of the dust dustsolub = ', dustsolub
296	WRITE(numout,*) ' Mineral Fe content of the dust mfrac = ', mfrac
297	WRITE(numout,*) ' Iron concentration in sea ice icefeinput = ', icefeinput
298	WRITE(numout,*) ' sinking speed of the dust wdust = ', wdust
299	WRITE(numout,*) ' nitrogen fixation rate nitrfix = ', nitrfix
300	WRITE(numout,*) ' nitrogen fixation sensitivty to light diazolight = ', diazolight
301	WRITE(numout,*) ' fe half-saturation cste for diazotrophs concfediaz = ', concfediaz
302	WRITE(numout,*) ' Fe to 3He ratio assumed for vent iron supply hratio = ', hratio
303	IF( ln_ligand ) THEN
304	WRITE(numout,*) ' Fep/Fer ratio from sed sources fep_rats = ', fep_rats
305	WRITE(numout,*) ' Fep/Fer ratio from sed hydro sources fep_rath = ', fep_rath
306	WRITE(numout,*) ' Weak ligand ratio from sed hydro sources lgw_rath = ', lgw_rath
307	ENDIF
308	END IF
309
310	IF( ln_dust .OR. ln_river .OR. ln_ndepo ) THEN ; ll_sbc = .TRUE.
311	ELSE ; ll_sbc = .FALSE.
312	ENDIF
313
314	IF( ln_dust .AND. ln_solub ) THEN ; ll_solub = .TRUE.
315	ELSE ; ll_solub = .FALSE.
316	ENDIF
317
318	! set the number of level over which river runoffs are applied
319	! online configuration : computed in sbcrnf
320	IF( l_offline ) THEN
321	!$OMP PARALLEL DO schedule(static) private(jj, ji)
322	DO jj = 1, jpj
323	DO ji = 1, jpi
324	nk_rnf(ji,jj) = 1
325	h_rnf (ji,jj) = gdept_n(ji,jj,1)
326	END DO
327	END DO
328	ENDIF
329
330	! dust input from the atmosphere
331	! ------------------------------
332	IF( ln_dust ) THEN
333	!
334	IF(lwp) WRITE(numout,*) ' initialize dust input from atmosphere '
335	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
336	!
337	ALLOCATE( dust(jpi,jpj) ) ! allocation
338	!
339	ALLOCATE( sf_dust(1), STAT=ierr ) !* allocate and fill sf_sst (forcing structure) with sn_sst
340	IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_dust structure' )
341	!
342	CALL fld_fill( sf_dust, (/ sn_dust /), cn_dir, 'p4z_sed_init', 'Atmospheric dust deposition', 'nampissed' )
343	ALLOCATE( sf_dust(1)%fnow(jpi,jpj,1) )
344	IF( sn_dust%ln_tint ) ALLOCATE( sf_dust(1)%fdta(jpi,jpj,1,2) )
345	!
346	IF( Agrif_Root() ) THEN ! Only on the master grid
347	! Get total input dust ; need to compute total atmospheric supply of Si in a year
348	CALL iom_open ( TRIM( sn_dust%clname ) , numdust )
349	CALL iom_gettime( numdust, zsteps, kntime=ntimes_dust) ! get number of record in file
350	ALLOCATE( zdust(jpi,jpj,ntimes_dust) )
351	DO jm = 1, ntimes_dust
352	CALL iom_get( numdust, jpdom_data, TRIM( sn_dust%clvar ), zdust(:,:,jm), jm )
353	END DO
354	CALL iom_close( numdust )
355	ztimes_dust = 1._wp / REAL(ntimes_dust, wp)
356	sumdepsi = 0.e0
357	DO jm = 1, ntimes_dust
358	sumdepsi = sumdepsi + glob_sum( zdust(:,:,jm) * e1e2t(:,:) * tmask(:,:,1) * ztimes_dust )
359	ENDDO
360	sumdepsi = sumdepsi / ( nyear_len(1) * rday ) * 12. * 8.8 * 0.075 * mfrac / 28.1
361	DEALLOCATE( zdust)
362	ENDIF
363	ELSE
364	sumdepsi = 0._wp
365	END IF
366
367	! Solubility of dust deposition of iron
368	! Only if ln_dust and ln_solubility set to true (ll_solub = .true.)
369	! -----------------------------------------------------------------
370	IF( ll_solub ) THEN
371	!
372	IF(lwp) WRITE(numout,*) ' initialize variable solubility of Fe '
373	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
374	!
375	ALLOCATE( solub(jpi,jpj) ) ! allocation
376	!
377	ALLOCATE( sf_solub(1), STAT=ierr ) !* allocate and fill sf_sst (forcing structure) with sn_sst
378	IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_solub structure' )
379	!
380	CALL fld_fill( sf_solub, (/ sn_solub /), cn_dir, 'p4z_sed_init', 'Solubility of atm. iron ', 'nampissed' )
381	ALLOCATE( sf_solub(1)%fnow(jpi,jpj,1) )
382	IF( sn_solub%ln_tint ) ALLOCATE( sf_solub(1)%fdta(jpi,jpj,1,2) )
383	! get number of record in file
384	CALL iom_open ( TRIM( sn_solub%clname ) , numsolub )
385	CALL iom_gettime( numsolub, zsteps, kntime=ntimes_solub) ! get number of record in file
386	CALL iom_close( numsolub )
387	ENDIF
388
389	! nutrient input from rivers
390	! --------------------------
391	IF( ln_river ) THEN
392	!
393	slf_river(jr_dic) = sn_riverdic ; slf_river(jr_doc) = sn_riverdoc ; slf_river(jr_din) = sn_riverdin
394	slf_river(jr_don) = sn_riverdon ; slf_river(jr_dip) = sn_riverdip ; slf_river(jr_dop) = sn_riverdop
395	slf_river(jr_dsi) = sn_riverdsi
396	!
397	ALLOCATE( rivdic(jpi,jpj), rivalk(jpi,jpj), rivdin(jpi,jpj), rivdip(jpi,jpj), rivdsi(jpi,jpj) )
398	IF( ln_p5z ) ALLOCATE( rivdon(jpi,jpj), rivdop(jpi,jpj), rivdoc(jpi,jpj) )
399	!
400	ALLOCATE( sf_river(jpriv), rivinput(jpriv), STAT=ierr1 ) !* allocate and fill sf_river (forcing structure) with sn_river_
401	rivinput(:) = 0.0
402
403	IF( ierr1 > 0 ) CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_irver structure' )
404	!
405	CALL fld_fill( sf_river, slf_river, cn_dir, 'p4z_sed_init', 'Input from river ', 'nampissed' )
406	DO ifpr = 1, jpriv
407	ALLOCATE( sf_river(ifpr)%fnow(jpi,jpj,1 ) )
408	IF( slf_river(ifpr)%ln_tint ) ALLOCATE( sf_river(ifpr)%fdta(jpi,jpj,1,2) )
409	END DO
410	IF( Agrif_Root() ) THEN ! Only on the master grid
411	! Get total input rivers ; need to compute total river supply in a year
412	DO ifpr = 1, jpriv
413	CALL iom_open ( TRIM( slf_river(ifpr)%clname ), numriv )
414	CALL iom_gettime( numriv, zsteps, kntime=ntimes_riv)
415	ALLOCATE( zriver(jpi,jpj,ntimes_riv) )
416	DO jm = 1, ntimes_riv
417	CALL iom_get( numriv, jpdom_data, TRIM( slf_river(ifpr)%clvar ), zriver(:,:,jm), jm )
418	END DO
419	CALL iom_close( numriv )
420	ztimes_riv = 1._wp / REAL(ntimes_riv, wp)
421	DO jm = 1, ntimes_riv
422	rivinput(ifpr) = rivinput(ifpr) + glob_sum( zriver(:,:,jm) * tmask(:,:,1) * ztimes_riv )
423	END DO
424	DEALLOCATE( zriver)
425	END DO
426	! N/P and Si releases due to coastal rivers
427	! -----------------------------------------
428	rivdicinput = (rivinput(jr_dic) + rivinput(jr_doc) ) * 1E3 / 12._wp
429	rivdininput = (rivinput(jr_din) + rivinput(jr_don) ) * 1E3 / rno3 / 14._wp
430	rivdipinput = (rivinput(jr_dip) + rivinput(jr_dop) ) * 1E3 / po4r / 31._wp
431	rivdsiinput = rivinput(jr_dsi) * 1E3 / 28.1_wp
432	rivalkinput = rivinput(jr_dic) * 1E3 / 12._wp
433	!
434	ENDIF
435	ELSE
436	rivdicinput = 0._wp
437	rivdininput = 0._wp
438	rivdipinput = 0._wp
439	rivdsiinput = 0._wp
440	rivalkinput = 0._wp
441	END IF
442	! nutrient input from dust
443	! ------------------------
444	IF( ln_ndepo ) THEN
445	!
446	IF(lwp) WRITE(numout,*) ' initialize the nutrient input by dust from ndeposition.orca.nc'
447	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
448	!
449	ALLOCATE( nitdep(jpi,jpj) ) ! allocation
450	!
451	ALLOCATE( sf_ndepo(1), STAT=ierr3 ) !* allocate and fill sf_sst (forcing structure) with sn_sst
452	IF( ierr3 > 0 ) CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_ndepo structure' )
453	!
454	CALL fld_fill( sf_ndepo, (/ sn_ndepo /), cn_dir, 'p4z_sed_init', 'Nutrient atmospheric depositon ', 'nampissed' )
455	ALLOCATE( sf_ndepo(1)%fnow(jpi,jpj,1) )
456	IF( sn_ndepo%ln_tint ) ALLOCATE( sf_ndepo(1)%fdta(jpi,jpj,1,2) )
457	!
458	IF( Agrif_Root() ) THEN ! Only on the master grid
459	! Get total input dust ; need to compute total atmospheric supply of N in a year
460	CALL iom_open ( TRIM( sn_ndepo%clname ), numdepo )
461	CALL iom_gettime( numdepo, zsteps, kntime=ntimes_ndep)
462	ALLOCATE( zndepo(jpi,jpj,ntimes_ndep) )
463	DO jm = 1, ntimes_ndep
464	CALL iom_get( numdepo, jpdom_data, TRIM( sn_ndepo%clvar ), zndepo(:,:,jm), jm )
465	END DO
466	CALL iom_close( numdepo )
467	ztimes_ndep = 1._wp / REAL(ntimes_ndep, wp)
468	nitdepinput = 0._wp
469	DO jm = 1, ntimes_ndep
470	nitdepinput = nitdepinput + glob_sum( zndepo(:,:,jm) * e1e2t(:,:) * tmask(:,:,1) * ztimes_ndep )
471	ENDDO
472	nitdepinput = nitdepinput / rno3 / 14E6
473	DEALLOCATE( zndepo)
474	ENDIF
475	ELSE
476	nitdepinput = 0._wp
477	ENDIF
478
479	! coastal and island masks
480	! ------------------------
481	IF( ln_ironsed ) THEN
482	!
483	IF(lwp) WRITE(numout,*) ' computation of an island mask to enhance coastal supply of iron'
484	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
485	!
486	ALLOCATE( ironsed(jpi,jpj,jpk) ) ! allocation
487	!
488	CALL iom_open ( TRIM( sn_ironsed%clname ), numiron )
489	ALLOCATE( zcmask(jpi,jpj,jpk) )
490	CALL iom_get ( numiron, jpdom_data, TRIM( sn_ironsed%clvar ), zcmask(:,:,:), 1 )
491	CALL iom_close( numiron )
492	!
493	ik50 = 5 ! last level where depth less than 50 m
494	DO jk = jpkm1, 1, -1
495	IF( gdept_1d(jk) > 50. ) ik50 = jk - 1
496	END DO
497	IF (lwp) WRITE(numout,*)
498	IF (lwp) WRITE(numout,*) ' Level corresponding to 50m depth ', ik50,' ', gdept_1d(ik50+1)
499	IF (lwp) WRITE(numout,*)
500	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zmaskt)
501	DO jk = 1, ik50
502	DO jj = 2, jpjm1
503	DO ji = fs_2, fs_jpim1
504	IF( tmask(ji,jj,jk) /= 0. ) THEN
505	zmaskt = tmask(ji+1,jj,jk) * tmask(ji-1,jj,jk) * tmask(ji,jj+1,jk) &
506	& * tmask(ji,jj-1,jk) * tmask(ji,jj,jk+1)
507	IF( zmaskt == 0. ) zcmask(ji,jj,jk ) = MAX( 0.1, zcmask(ji,jj,jk) )
508	END IF
509	END DO
510	END DO
511	END DO
512	!
513	CALL lbc_lnk( zcmask , 'T', 1. ) ! lateral boundary conditions on cmask (sign unchanged)
514	!
515	!$OMP PARALLEL
516	!$OMP DO schedule(static) private(jk, jj, ji, zexpide, zdenitide)
517	DO jk = 1, jpk
518	DO jj = 1, jpj
519	DO ji = 1, jpi
520	zexpide = MIN( 8.,( gdept_n(ji,jj,jk) / 500. )**(-1.5) )
521	zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2
522	zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 )
523	END DO
524	END DO
525	END DO
526	!$OMP END DO NOWAIT
527	! Coastal supply of iron
528	! -------------------------
529	!$OMP DO schedule(static) private(jj,ji)
530	DO jj = 1, jpj
531	DO ji = 1, jpi
532	ironsed(ji,jj,jpk) = 0._wp
533	END DO
534	END DO
535	!$OMP DO schedule(static) private(jk,jj,ji)
536	DO jk = 1, jpkm1
537	DO jj = 1, jpj
538	DO ji = 1, jpi
539	ironsed(ji,jj,jk) = sedfeinput * zcmask(ji,jj,jk) / ( e3t_0(ji,jj,jk) * rday )
540	END DO
541	END DO
542	END DO
543	!$OMP END PARALLEL
544	DEALLOCATE( zcmask)
545	ENDIF
546	!
547	! Iron from Hydrothermal vents
548	! ------------------------
549	IF( ln_hydrofe ) THEN
550	!
551	IF(lwp) WRITE(numout,*) ' Input of iron from hydrothermal vents '
552	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
553	!
554	ALLOCATE( hydrofe(jpi,jpj,jpk) ) ! allocation
555	!
556	CALL iom_open ( TRIM( sn_hydrofe%clname ), numhydro )
557	CALL iom_get ( numhydro, jpdom_data, TRIM( sn_hydrofe%clvar ), hydrofe(:,:,:), 1 )
558	CALL iom_close( numhydro )
559	!
560	DO jk = 1, jpk
561	hydrofe(:,:,jk) = ( hydrofe(:,:,jk) * hratio ) / ( e1e2t(:,:) * e3t_0(:,:,jk) * ryyss + rtrn ) / 1000._wp
562	ENDDO
563	!
564	ENDIF
565	!
566	IF( ll_sbc ) CALL p4z_sbc( nit000 )
567	!
568	IF(lwp) THEN
569	WRITE(numout,*)
570	WRITE(numout,*) ' Total input of elements from river supply'
571	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
572	WRITE(numout,) ' N Supply : ', rivdininputrno31E3/1E1214.,' TgN/yr'
573	WRITE(numout,) ' Si Supply : ', rivdsiinput1E3/1E12*28.1,' TgSi/yr'
574	WRITE(numout,) ' P Supply : ', rivdipinput1E3po4r/1E1231.,' TgP/yr'
575	WRITE(numout,) ' Alk Supply : ', rivalkinput1E3/1E12,' Teq/yr'
576	WRITE(numout,) ' DIC Supply : ', rivdicinput1E3*12./1E12,'TgC/yr'
577	WRITE(numout,*)
578	WRITE(numout,*) ' Total input of elements from atmospheric supply'
579	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
580	WRITE(numout,) ' N Supply : ', nitdepinputrno31E3/1E1214.,' TgN/yr'
581	WRITE(numout,*)
582	ENDIF
583	!
584	IF( nn_timing == 1 ) CALL timing_stop('p4z_sbc_init')
585	!
586	END SUBROUTINE p4z_sbc_init
587
588	!!======================================================================
589	END MODULE p4zsbc

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

Download in other formats: