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 NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/TOP/PISCES/P4Z – NEMO

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/TOP/PISCES/P4Z/p4zsbc.F90 @ 11671

Last change on this file since 11671 was 11671, checked in by acc, 5 years ago

Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Final, non-substantive changes to complete this branch. These changes remove all REWIND statements on the old namelist fortran units (now character variables for internal files). These changes have been left until last since they are easily repeated via a script and it may be preferable to use the previous revision for merge purposes and reapply these last changes separately. This branch has been fully SETTE tested.
  • Property svn:keywords set to Id
File size: 25.5 KB
Line 
1MODULE 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   REAL(wp), PUBLIC ::   sedfeinput   !: Coastal release of Iron
30   REAL(wp), PUBLIC ::   dustsolub    !: Solubility of the dust
31   REAL(wp), PUBLIC ::   mfrac        !: Mineral Content of the dust
32   REAL(wp), PUBLIC ::   icefeinput   !: Iron concentration in sea ice
33   REAL(wp), PUBLIC ::   wdust        !: Sinking speed of the dust
34   REAL(wp), PUBLIC ::   nitrfix      !: Nitrogen fixation rate   
35   REAL(wp), PUBLIC ::   diazolight   !: Nitrogen fixation sensitivty to light
36   REAL(wp), PUBLIC ::   concfediaz   !: Fe half-saturation Cste for diazotrophs
37   REAL(wp)         ::   hratio       !: Fe:3He ratio assumed for vent iron supply
38   REAL(wp)         ::   distcoast    !: Distance off the coast for Iron from sediments
39   REAL(wp), PUBLIC ::   lgw_rath     !: Weak ligand ratio from hydro sources
40
41   LOGICAL , PUBLIC ::   ll_sbc
42   LOGICAL          ::   ll_solub
43
44   INTEGER , PARAMETER  :: jpriv  = 7   !: Maximum number of river input fields
45   INTEGER , PARAMETER  :: jr_dic = 1   !: index of dissolved inorganic carbon
46   INTEGER , PARAMETER  :: jr_doc = 2   !: index of dissolved organic carbon
47   INTEGER , PARAMETER  :: jr_din = 3   !: index of dissolved inorganic nitrogen
48   INTEGER , PARAMETER  :: jr_don = 4   !: index of dissolved organic nitrogen
49   INTEGER , PARAMETER  :: jr_dip = 5   !: index of dissolved inorganic phosporus
50   INTEGER , PARAMETER  :: jr_dop = 6   !: index of dissolved organic phosphorus
51   INTEGER , PARAMETER  :: jr_dsi = 7   !: index of dissolved silicate
52
53   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_dust      ! structure of input dust
54   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_solub     ! structure of input dust
55   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_river     ! structure of input riverdic
56   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_ndepo     ! structure of input nitrogen deposition
57   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_ironsed   ! structure of input iron from sediment
58   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_hydrofe   ! structure of input iron from hydrothermal vents
59
60   INTEGER , PARAMETER ::   nbtimes = 365                          ! maximum number of times record in a file
61   INTEGER             ::   ntimes_dust, ntimes_riv, ntimes_ndep   ! number of time steps in a file
62   INTEGER             ::   ntimes_solub, ntimes_hydro             ! number of time steps in a file
63
64   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dust  , solub    !: dust fields
65   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rivdic, rivalk   !: river input fields
66   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rivdin, rivdip   !: river input fields
67   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rivdon, rivdop   !: river input fields
68   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rivdoc           !: river input fields
69   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rivdsi           !: river input fields
70   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   nitdep           !: atmospheric N deposition
71   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ironsed          !: Coastal supply of iron
72   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hydrofe          !: Hydrothermal vent supply of iron
73
74   REAL(wp), PUBLIC :: sedsilfrac, sedcalfrac
75   REAL(wp), PUBLIC :: rivalkinput, rivdicinput
76   REAL(wp), PUBLIC :: rivdininput, rivdipinput, rivdsiinput
77
78   !! * Substitutions
79#  include "vectopt_loop_substitute.h90"
80   !!----------------------------------------------------------------------
81   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
82   !! $Id$
83   !! Software governed by the CeCILL license (see ./LICENSE)
84   !!----------------------------------------------------------------------
85CONTAINS
86
87   SUBROUTINE p4z_sbc( kt )
88      !!----------------------------------------------------------------------
89      !!                  ***  routine p4z_sbc  ***
90      !!
91      !! ** purpose :   read and interpolate the external sources of nutrients
92      !!
93      !! ** method  :   read the files and interpolate the appropriate variables
94      !!
95      !! ** input   :   external netcdf files
96      !!
97      !!----------------------------------------------------------------------
98      INTEGER, INTENT(in) ::   kt   ! ocean time step
99      !
100      INTEGER  ::   ji, jj 
101      REAL(wp) ::   zcoef, zyyss
102      !!---------------------------------------------------------------------
103      !
104      IF( ln_timing )   CALL timing_start('p4z_sbc')
105      !
106      ! Compute dust at nit000 or only if there is more than 1 time record in dust file
107      IF( ln_dust ) THEN
108         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_dust > 1 ) ) THEN
109            CALL fld_read( kt, 1, sf_dust )
110            dust(:,:) = MAX( rtrn, sf_dust(1)%fnow(:,:,1) ) * ( 1.0 - fr_i(:,:) )
111         ENDIF
112      ENDIF
113      !
114      IF( ll_solub ) THEN
115         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_solub > 1 ) ) THEN
116            CALL fld_read( kt, 1, sf_solub )
117            solub(:,:) = sf_solub(1)%fnow(:,:,1)
118         ENDIF
119      ENDIF
120
121      ! N/P and Si releases due to coastal rivers
122      ! Compute river at nit000 or only if there is more than 1 time record in river file
123      ! -----------------------------------------
124      IF( ln_river ) THEN
125         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_riv > 1 ) ) THEN
126            CALL fld_read( kt, 1, sf_river )
127            IF( ln_p4z ) THEN
128               DO jj = 1, jpj
129                  DO ji = 1, jpi
130                     zcoef = ryyss * e1e2t(ji,jj) * h_rnf(ji,jj) 
131                     rivalk(ji,jj) =   sf_river(jr_dic)%fnow(ji,jj,1)  &
132                        &              * 1.E3        / ( 12. * zcoef + rtrn )
133                     rivdic(ji,jj) =   sf_river(jr_dic)%fnow(ji,jj,1)  &
134                        &              * 1.E3         / ( 12. * zcoef + rtrn )
135                     rivdin(ji,jj) =   sf_river(jr_din)%fnow(ji,jj,1)  &
136                        &              * 1.E3 / rno3 / ( 14. * zcoef + rtrn )
137                     rivdip(ji,jj) =   sf_river(jr_dip)%fnow(ji,jj,1)  &
138                        &              * 1.E3 / po4r / ( 31. * zcoef + rtrn )
139                     rivdsi(ji,jj) =   sf_river(jr_dsi)%fnow(ji,jj,1)  &
140                        &              * 1.E3        / ( 28.1 * zcoef + rtrn )
141                     rivdoc(ji,jj) =   sf_river(jr_doc)%fnow(ji,jj,1)  &
142                        &              * 1.E3        / ( 12. * zcoef + rtrn ) 
143                  END DO
144               END DO
145            ELSE    !  ln_p5z
146               DO jj = 1, jpj
147                  DO ji = 1, jpi
148                     zcoef = ryyss * e1e2t(ji,jj) * h_rnf(ji,jj) 
149                     rivalk(ji,jj) =   sf_river(jr_dic)%fnow(ji,jj,1)                                    &
150                        &              * 1.E3        / ( 12. * zcoef + rtrn )
151                     rivdic(ji,jj) = ( sf_river(jr_dic)%fnow(ji,jj,1) ) &
152                        &              * 1.E3 / ( 12. * zcoef + rtrn ) * tmask(ji,jj,1)
153                     rivdin(ji,jj) = ( sf_river(jr_din)%fnow(ji,jj,1) ) &
154                        &              * 1.E3 / rno3 / ( 14. * zcoef + rtrn ) * tmask(ji,jj,1)
155                     rivdip(ji,jj) = ( sf_river(jr_dip)%fnow(ji,jj,1) ) &
156                        &              * 1.E3 / po4r / ( 31. * zcoef + rtrn ) * tmask(ji,jj,1)
157                     rivdon(ji,jj) = ( sf_river(jr_don)%fnow(ji,jj,1) ) &
158                        &              * 1.E3 / rno3 / ( 14. * zcoef + rtrn ) * tmask(ji,jj,1)
159                     rivdop(ji,jj) = ( sf_river(jr_dop)%fnow(ji,jj,1) ) &
160                        &              * 1.E3 / po4r / ( 31. * zcoef + rtrn ) * tmask(ji,jj,1)
161                     rivdsi(ji,jj) =   sf_river(jr_dsi)%fnow(ji,jj,1)  &
162                        &              * 1.E3        / ( 28.1 * zcoef + rtrn )
163                     rivdoc(ji,jj) =   sf_river(jr_doc)%fnow(ji,jj,1)  &
164                        &              * 1.E3        / ( 12. * zcoef + rtrn )
165                  END DO
166               END DO
167            ENDIF
168         ENDIF
169      ENDIF
170
171      ! Compute N deposition at nit000 or only if there is more than 1 time record in N deposition file
172      IF( ln_ndepo ) THEN
173         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_ndep > 1 ) ) THEN
174             zcoef = 14. * rno3
175             CALL fld_read( kt, 1, sf_ndepo )
176             nitdep(:,:) = MAX( rtrn, sf_ndepo(1)%fnow(:,:,1) / zcoef / e3t_n(:,:,1) )
177         ENDIF
178         IF( .NOT.ln_linssh ) THEN
179            zcoef = 14. * rno3
180            nitdep(:,:) = MAX( rtrn, sf_ndepo(1)%fnow(:,:,1) / zcoef / e3t_n(:,:,1) )
181         ENDIF
182      ENDIF
183      !
184      IF( ln_timing )  CALL timing_stop('p4z_sbc')
185      !
186   END SUBROUTINE p4z_sbc
187
188
189   SUBROUTINE p4z_sbc_init
190      !!----------------------------------------------------------------------
191      !!                  ***  routine p4z_sbc_init  ***
192      !!
193      !! ** purpose :   initialization of the external sources of nutrients
194      !!
195      !! ** method  :   read the files and compute the budget
196      !!                called at the first timestep (nittrc000)
197      !!
198      !! ** input   :   external netcdf files
199      !!
200      !!----------------------------------------------------------------------
201      INTEGER  :: ji, jj, jk, jm, ifpr
202      INTEGER  :: ii0, ii1, ij0, ij1
203      INTEGER  :: numdust, numsolub, numriv, numiron, numdepo, numhydro
204      INTEGER  :: ierr, ierr1, ierr2, ierr3
205      INTEGER  :: ios                 ! Local integer output status for namelist read
206      INTEGER  :: ik50                !  last level where depth less than 50 m
207      INTEGER  :: isrow             ! index for ORCA1 starting row
208      REAL(wp) :: zexpide, zdenitide, zmaskt, zsurfc, zsurfp,ze3t, ze3t2, zcslp
209      REAL(wp) :: ztimes_dust, ztimes_riv, ztimes_ndep 
210      REAL(wp), DIMENSION(:), ALLOCATABLE :: rivinput
211      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zriver, zcmask
212      !
213      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
214      TYPE(FLD_N), DIMENSION(jpriv) ::  slf_river    ! array of namelist informations on the fields to read
215      TYPE(FLD_N) ::   sn_dust, sn_solub, sn_ndepo, sn_ironsed, sn_hydrofe   ! informations about the fields to be read
216      TYPE(FLD_N) ::   sn_riverdoc, sn_riverdic, sn_riverdsi   ! informations about the fields to be read
217      TYPE(FLD_N) ::   sn_riverdin, sn_riverdon, sn_riverdip, sn_riverdop
218      !!
219      NAMELIST/nampissbc/cn_dir, sn_dust, sn_solub, sn_riverdic, sn_riverdoc, sn_riverdin, sn_riverdon,     &
220        &                sn_riverdip, sn_riverdop, sn_riverdsi, sn_ndepo, sn_ironsed, sn_hydrofe, &
221        &                ln_dust, ln_solub, ln_river, ln_ndepo, ln_ironsed, ln_ironice, ln_hydrofe,    &
222        &                sedfeinput, distcoast, dustsolub, icefeinput, wdust, mfrac, nitrfix, diazolight, concfediaz, &
223        &                hratio, lgw_rath
224      !!----------------------------------------------------------------------
225      !
226      IF(lwp) THEN
227         WRITE(numout,*)
228         WRITE(numout,*) 'p4z_sbc_init : initialization of the external sources of nutrients '
229         WRITE(numout,*) '~~~~~~~~~~~~ '
230      ENDIF
231      !                            !* set file information
232      READ  ( numnatp_ref, nampissbc, IOSTAT = ios, ERR = 901)
233901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampissbc in reference namelist' )
234      READ  ( numnatp_cfg, nampissbc, IOSTAT = ios, ERR = 902 )
235902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampissbc in configuration namelist' )
236      IF(lwm) WRITE ( numonp, nampissbc )
237
238      IF(lwp) THEN
239         WRITE(numout,*) '   Namelist : nampissbc '
240         WRITE(numout,*) '      dust input from the atmosphere           ln_dust     = ', ln_dust
241         WRITE(numout,*) '      Variable solubility of iron input        ln_solub    = ', ln_solub
242         WRITE(numout,*) '      river input of nutrients                 ln_river    = ', ln_river
243         WRITE(numout,*) '      atmospheric deposition of n              ln_ndepo    = ', ln_ndepo
244         WRITE(numout,*) '      Fe input from sediments                  ln_ironsed  = ', ln_ironsed
245         WRITE(numout,*) '      Fe input from seaice                     ln_ironice  = ', ln_ironice
246         WRITE(numout,*) '      fe input from hydrothermal vents         ln_hydrofe  = ', ln_hydrofe
247         WRITE(numout,*) '      coastal release of iron                  sedfeinput  = ', sedfeinput
248         WRITE(numout,*) '      distance off the coast                   distcoast   = ', distcoast
249         WRITE(numout,*) '      solubility of the dust                   dustsolub   = ', dustsolub
250         WRITE(numout,*) '      Mineral Fe content of the dust           mfrac       = ', mfrac
251         WRITE(numout,*) '      Iron concentration in sea ice            icefeinput  = ', icefeinput
252         WRITE(numout,*) '      sinking speed of the dust                wdust       = ', wdust
253         WRITE(numout,*) '      nitrogen fixation rate                   nitrfix     = ', nitrfix
254         WRITE(numout,*) '      nitrogen fixation sensitivty to light    diazolight  = ', diazolight
255         WRITE(numout,*) '      Fe half-saturation cste for diazotrophs  concfediaz  = ', concfediaz
256         WRITE(numout,*) '      Fe to 3He ratio assumed for vent iron supply hratio  = ', hratio
257         IF( ln_ligand ) THEN
258            WRITE(numout,*) '      Weak ligand ratio from sed hydro sources  lgw_rath   = ', lgw_rath
259         ENDIF
260      END IF
261
262      IF( ln_dust .OR. ln_river .OR. ln_ndepo ) THEN   ;   ll_sbc = .TRUE.
263      ELSE                                             ;   ll_sbc = .FALSE.
264      ENDIF
265
266      IF( ln_dust .AND. ln_solub ) THEN                ;   ll_solub = .TRUE.
267      ELSE                                             ;   ll_solub = .FALSE.
268      ENDIF
269
270      ! set the number of level over which river runoffs are applied
271      ! online configuration : computed in sbcrnf
272      IF( l_offline ) THEN
273        nk_rnf(:,:) = 1
274        h_rnf (:,:) = gdept_n(:,:,1)
275      ENDIF
276
277      ! dust input from the atmosphere
278      ! ------------------------------
279      IF( ln_dust ) THEN 
280         !
281         IF(lwp) WRITE(numout,*) '    initialize dust input from atmosphere '
282         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
283         !
284         ALLOCATE( dust(jpi,jpj) )    ! allocation
285         !
286         ALLOCATE( sf_dust(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst
287         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_dust structure' )
288         !
289         CALL fld_fill( sf_dust, (/ sn_dust /), cn_dir, 'p4z_sed_init', 'Atmospheric dust deposition', 'nampissed' )
290                                   ALLOCATE( sf_dust(1)%fnow(jpi,jpj,1)   )
291         IF( sn_dust%ln_tint )     ALLOCATE( sf_dust(1)%fdta(jpi,jpj,1,2) )
292         !
293         IF( Agrif_Root() ) THEN   !  Only on the master grid
294            ! Get total input dust ; need to compute total atmospheric supply of Si in a year
295            CALL iom_open (  TRIM( sn_dust%clname ) , numdust )
296            ntimes_dust = iom_getszuld( numdust )   ! get number of record in file
297         END IF
298      END IF
299
300      ! Solubility of dust deposition of iron
301      ! Only if ln_dust and ln_solubility set to true (ll_solub = .true.)
302      ! -----------------------------------------------------------------
303      IF( ll_solub ) THEN
304         !
305         IF(lwp) WRITE(numout,*)
306         IF(lwp) WRITE(numout,*) '   ==>>>   ll_solub=T , initialize variable solubility of Fe '
307         !
308         ALLOCATE( solub(jpi,jpj) )    ! allocation
309         !
310         ALLOCATE( sf_solub(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst
311         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_solub structure' )
312         !
313         CALL fld_fill( sf_solub, (/ sn_solub /), cn_dir, 'p4z_sed_init', 'Solubility of atm. iron ', 'nampissed' )
314                                   ALLOCATE( sf_solub(1)%fnow(jpi,jpj,1)   )
315         IF( sn_solub%ln_tint )    ALLOCATE( sf_solub(1)%fdta(jpi,jpj,1,2) )
316         ! get number of record in file
317         CALL iom_open (  TRIM( sn_solub%clname ) , numsolub )
318         ntimes_solub = iom_getszuld( numsolub )   ! get number of record in file
319         CALL iom_close( numsolub )
320      ENDIF
321
322      ! nutrient input from rivers
323      ! --------------------------
324      IF( ln_river ) THEN
325         !
326         slf_river(jr_dic) = sn_riverdic   ;   slf_river(jr_doc) = sn_riverdoc   ;   slf_river(jr_din) = sn_riverdin 
327         slf_river(jr_don) = sn_riverdon   ;   slf_river(jr_dip) = sn_riverdip   ;   slf_river(jr_dop) = sn_riverdop
328         slf_river(jr_dsi) = sn_riverdsi 
329         !
330         ALLOCATE( rivdic(jpi,jpj), rivalk(jpi,jpj), rivdin(jpi,jpj), rivdip(jpi,jpj), rivdsi(jpi,jpj), rivdoc(jpi,jpj) ) 
331         IF( ln_p5z )  ALLOCATE( rivdon(jpi,jpj), rivdop(jpi,jpj) )
332         !
333         ALLOCATE( sf_river(jpriv), rivinput(jpriv), STAT=ierr1 )    !* allocate and fill sf_river (forcing structure) with sn_river_
334         rivinput(:) = 0._wp
335
336         IF( ierr1 > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_irver structure' )
337         !
338         CALL fld_fill( sf_river, slf_river, cn_dir, 'p4z_sed_init', 'Input from river ', 'nampissed' )
339         DO ifpr = 1, jpriv
340                                          ALLOCATE( sf_river(ifpr)%fnow(jpi,jpj,1  ) )
341            IF( slf_river(ifpr)%ln_tint ) ALLOCATE( sf_river(ifpr)%fdta(jpi,jpj,1,2) )
342         END DO
343         IF( Agrif_Root() ) THEN   !  Only on the master grid
344            ! Get total input rivers ; need to compute total river supply in a year
345            DO ifpr = 1, jpriv
346               CALL iom_open ( TRIM( slf_river(ifpr)%clname ), numriv )
347               ntimes_riv = iom_getszuld( numriv )
348               ALLOCATE( zriver(jpi,jpj,ntimes_riv) )
349               DO jm = 1, ntimes_riv
350                  CALL iom_get( numriv, jpdom_data, TRIM( slf_river(ifpr)%clvar ), zriver(:,:,jm), jm )
351               END DO
352               CALL iom_close( numriv )
353               ztimes_riv = 1._wp / REAL(ntimes_riv, wp) 
354               DO jm = 1, ntimes_riv
355                  rivinput(ifpr) = rivinput(ifpr) + glob_sum( 'p4zsbc', zriver(:,:,jm) * tmask(:,:,1) * ztimes_riv ) 
356               END DO
357               DEALLOCATE( zriver)
358            END DO
359            ! N/P and Si releases due to coastal rivers
360            ! -----------------------------------------
361            rivdicinput = (rivinput(jr_dic) + rivinput(jr_doc) ) * 1E3 / 12._wp
362            rivdininput = (rivinput(jr_din) + rivinput(jr_don) ) * 1E3 / rno3 / 14._wp
363            rivdipinput = (rivinput(jr_dip) + rivinput(jr_dop) ) * 1E3 / po4r / 31._wp
364            rivdsiinput = rivinput(jr_dsi) * 1E3 / 28.1_wp
365            rivalkinput = rivinput(jr_dic) * 1E3 / 12._wp
366            !
367         ENDIF
368      ELSE
369         rivdicinput = 0._wp
370         rivdininput = 0._wp
371         rivdipinput = 0._wp
372         rivdsiinput = 0._wp
373         rivalkinput = 0._wp
374      END IF 
375      ! nutrient input from dust
376      ! ------------------------
377      IF( ln_ndepo ) THEN
378         !
379         IF(lwp) WRITE(numout,*)
380         IF(lwp) WRITE(numout,*) '   ==>>>   ln_ndepo=T , initialize the nutrient input by dust from NetCDF file'
381         !
382         ALLOCATE( nitdep(jpi,jpj) )    ! allocation
383         !
384         ALLOCATE( sf_ndepo(1), STAT=ierr3 )           !* allocate and fill sf_sst (forcing structure) with sn_sst
385         IF( ierr3 > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_ndepo structure' )
386         !
387         CALL fld_fill( sf_ndepo, (/ sn_ndepo /), cn_dir, 'p4z_sed_init', 'Nutrient atmospheric depositon ', 'nampissed' )
388                                   ALLOCATE( sf_ndepo(1)%fnow(jpi,jpj,1)   )
389         IF( sn_ndepo%ln_tint )    ALLOCATE( sf_ndepo(1)%fdta(jpi,jpj,1,2) )
390         !
391         IF( Agrif_Root() ) THEN   !  Only on the master grid
392            ! Get total input dust ; need to compute total atmospheric supply of N in a year
393            CALL iom_open ( TRIM( sn_ndepo%clname ), numdepo )
394            ntimes_ndep = iom_getszuld( numdepo )
395         ENDIF
396      ENDIF
397
398      ! coastal and island masks
399      ! ------------------------
400      IF( ln_ironsed ) THEN     
401         !
402         IF(lwp) WRITE(numout,*)
403         IF(lwp) WRITE(numout,*) '   ==>>>   ln_ironsed=T , computation of an island mask to enhance coastal supply of iron'
404         !
405         ALLOCATE( ironsed(jpi,jpj,jpk) )    ! allocation
406         !
407         CALL iom_open ( TRIM( sn_ironsed%clname ), numiron )
408         ALLOCATE( zcmask(jpi,jpj,jpk) )
409         CALL iom_get  ( numiron, jpdom_data, TRIM( sn_ironsed%clvar ), zcmask(:,:,:), 1 )
410         CALL iom_close( numiron )
411         !
412         ik50 = 5        !  last level where depth less than 50 m
413         DO jk = jpkm1, 1, -1
414            IF( gdept_1d(jk) > 50. )   ik50 = jk - 1
415         END DO
416         IF(lwp) WRITE(numout,*)
417         IF(lwp) WRITE(numout,*) ' Level corresponding to 50m depth ',  ik50,' ', gdept_1d(ik50+1)
418         DO jk = 1, ik50
419            DO jj = 2, jpjm1
420               DO ji = fs_2, fs_jpim1
421                  ze3t   = e3t_0(ji,jj,jk)
422                  zsurfc =  e1u(ji,jj) * ( 1. - umask(ji  ,jj  ,jk) )   &
423                          + e1u(ji,jj) * ( 1. - umask(ji-1,jj  ,jk) )   &
424                          + e2v(ji,jj) * ( 1. - vmask(ji  ,jj  ,jk) )   &
425                          + e2v(ji,jj) * ( 1. - vmask(ji  ,jj-1,jk) )
426                  zsurfp = zsurfc * ze3t / e1e2t(ji,jj)
427                  ! estimation of the coastal slope : 5 km off the coast
428                  ze3t2 = ze3t * ze3t
429                  zcslp = SQRT( ( distcoast*distcoast + ze3t2 ) / ze3t2 )
430                  !
431                  zcmask(ji,jj,jk) = zcmask(ji,jj,jk) + zcslp * zsurfp
432               END DO
433            END DO
434         END DO
435         !
436         CALL lbc_lnk( 'p4zsbc', zcmask , 'T', 1. )      ! lateral boundary conditions on cmask   (sign unchanged)
437         !
438         DO jk = 1, jpk
439            DO jj = 1, jpj
440               DO ji = 1, jpi
441                  zexpide   = MIN( 8.,( gdept_n(ji,jj,jk) / 500. )**(-1.5) )
442                  zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2
443                  zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 )
444               END DO
445            END DO
446         END DO
447         ! Coastal supply of iron
448         ! -------------------------
449         ironsed(:,:,jpk) = 0._wp
450         DO jk = 1, jpkm1
451            ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( e3t_0(:,:,jk) * rday )
452         END DO
453         DEALLOCATE( zcmask)
454      ENDIF
455      !
456      ! Iron from Hydrothermal vents
457      ! ------------------------
458      IF( ln_hydrofe ) THEN
459         !
460         IF(lwp) WRITE(numout,*)
461         IF(lwp) WRITE(numout,*) '   ==>>>   ln_hydrofe=T , Input of iron from hydrothermal vents'
462         !
463         ALLOCATE( hydrofe(jpi,jpj,jpk) )    ! allocation
464         !
465         CALL iom_open ( TRIM( sn_hydrofe%clname ), numhydro )
466         CALL iom_get  ( numhydro, jpdom_data, TRIM( sn_hydrofe%clvar ), hydrofe(:,:,:), 1 )
467         CALL iom_close( numhydro )
468         !
469         DO jk = 1, jpk
470            hydrofe(:,:,jk) = ( hydrofe(:,:,jk) * hratio ) / ( e1e2t(:,:) * e3t_0(:,:,jk) * ryyss + rtrn ) / 1000._wp
471         ENDDO
472         !
473      ENDIF
474      !
475      IF( ll_sbc ) CALL p4z_sbc( nit000 ) 
476      !
477      IF(lwp) THEN
478         WRITE(numout,*)
479         WRITE(numout,*) '    Total input of elements from river supply'
480         WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
481         WRITE(numout,*) '    N Supply   : ', rivdininput*rno3*1E3/1E12*14.,' TgN/yr'
482         WRITE(numout,*) '    Si Supply  : ', rivdsiinput*1E3/1E12*28.1    ,' TgSi/yr'
483         WRITE(numout,*) '    P Supply   : ', rivdipinput*1E3*po4r/1E12*31.,' TgP/yr'
484         WRITE(numout,*) '    Alk Supply : ', rivalkinput*1E3/1E12         ,' Teq/yr'
485         WRITE(numout,*) '    DIC Supply : ', rivdicinput*1E3*12./1E12     ,' TgC/yr'
486         WRITE(numout,*) 
487      ENDIF
488      !
489      sedsilfrac = 0.03     ! percentage of silica loss in the sediments
490      sedcalfrac = 0.6      ! percentage of calcite loss in the sediments
491      !
492   END SUBROUTINE p4z_sbc_init
493
494   !!======================================================================
495END MODULE p4zsbc
Note: See TracBrowser for help on using the repository browser.