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.
p4zbc.F90 in NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/P4Z – NEMO

source: NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/P4Z/p4zbc.F90 @ 15075

Last change on this file since 15075 was 15075, checked in by aumont, 3 years ago

major update of the sediment module
File size: 18.1 KB
Line 
1MODULE p4zbc
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_bc        :  Read and interpolate time-varying nutrients fluxes
9   !!   p4z_bc_init   :  Initialization of p4z_bc
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   USE trcbc
17
18   IMPLICIT NONE
19   PRIVATE
20
21   PUBLIC   p4z_bc
22   PUBLIC   p4z_bc_init   
23
24   LOGICAL , PUBLIC ::   ln_ironsed   !: boolean for Fe input from sediments
25   LOGICAL , PUBLIC ::   ln_hydrofe   !: boolean for Fe input from hydrothermal vents
26   LOGICAL , PUBLIC ::   ln_dust_inp  !: boolean for Fe input from hydrothermal vents
27   REAL(wp), PUBLIC ::   sedfeinput   !: Coastal release of Iron
28   REAL(wp), PUBLIC ::   icefeinput   !: Iron concentration in sea ice
29   REAL(wp), PUBLIC ::   wdust        !: Sinking speed of the dust
30   REAL(wp), PUBLIC ::   mfrac        !: Mineral Content of the dust
31   REAL(wp)         ::   hratio       !: Fe:3He ratio assumed for vent iron supply
32   REAL(wp)         ::   distcoast    !: Distance off the coast for Iron from sediments
33   REAL(wp), PUBLIC ::   lgw_rath     !: Weak ligand ratio from hydro sources
34
35   LOGICAL , PUBLIC ::   ll_bc
36   LOGICAL , PUBLIC ::   ll_dust      !: boolean for dust input from the atmosphere
37   LOGICAL , PUBLIC ::   ll_river     !: boolean for river input of nutrients
38   LOGICAL , PUBLIC ::   ll_ndepo     !: boolean for atmospheric deposition of N
39   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_dust      ! structure of input dust
40   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_ironsed   ! structure of input iron from sediment
41   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_hydrofe   ! structure of input iron from sediment
42
43   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dust    !: dust fields
44   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ironsed          !: Coastal supply of iron
45   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hydrofe          !: Hydrothermal vent supply of iron
46
47   REAL(wp), PUBLIC :: sedsilfrac, sedcalfrac
48
49   !! * Substitutions
50#  include "do_loop_substitute.h90"
51#  include "domzgr_substitute.h90"
52   !!----------------------------------------------------------------------
53   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
54   !! $Id: p4zbc.F90 10869 2019-04-15 10:34:03Z cetlod $
55   !! Software governed by the CeCILL license (see ./LICENSE)
56   !!----------------------------------------------------------------------
57CONTAINS
58
59   SUBROUTINE p4z_bc( kt, Kbb, Kmm, Krhs )
60      !!----------------------------------------------------------------------
61      !!                  ***  routine p4z_bc  ***
62      !!
63      !! ** purpose :   read and interpolate the external sources of nutrients
64      !!
65      !! ** method  :   read the files and interpolate the appropriate variables
66      !!
67      !! ** input   :   external netcdf files
68      !!
69      !!----------------------------------------------------------------------
70      INTEGER, INTENT(in) ::   kt              ! ocean time step
71      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs  ! time level index
72      !
73      INTEGER  ::  ji, jj, jk, jl 
74      REAL(wp) ::  zdep, zwflux, zironice
75      REAL(wp) ::  zcoef, zwdust, zrivdin, zdustdep, zndep
76      !
77      CHARACTER (len=25) :: charout
78      !!---------------------------------------------------------------------
79      !
80      IF( ln_timing )   CALL timing_start('p4z_bc')
81     
82      ! Add the external input of nutrients from dust deposition in the water column
83      ! The inputs at surface have already been added
84      ! ----------------------------------------------------------
85      IF( ll_dust )  THEN
86         !
87         CALL fld_read( kt, 1, sf_dust )
88         dust(:,:) = MAX( rtrn, sf_dust(1)%fnow(:,:,1) )
89         !
90         ! Iron solubilization of particles in the water column
91         ! dust in kg/m2/s ---> 1/55.85 to put in mol/Fe ;  wdust in m/d
92         ! Dust are supposed to sink at wdust sinking speed. 3% of the iron
93         ! in dust is hypothesized to be soluble at a dissolution rate set to
94         ! 1/(250 days). The vertical distribution of iron in dust is computed
95         ! from a steady state assumption. Parameters are very uncertain and
96         ! are estimated from the literature quoted in Raiswell et al. (2011)
97         ! -------------------------------------------------------------------
98
99         zwdust = 0.03 / ( wdust / rday ) / ( 250. * rday )
100
101         ! Atmospheric input of Iron dissolves in the water column
102         IF ( ln_trc_sbc(jpfer) ) THEN
103            DO_3D( 1, 1, 1, 1, 2, jpkm1 )
104               zdustdep = dust(ji,jj) * zwdust * rfact * EXP( -gdept(ji,jj,jk,Kmm) /( 250. * wdust ) )
105               !
106               tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + zdustdep * mfrac / mMass_Fe 
107            END_3D
108
109            IF( lk_iomput ) THEN
110                ! surface downward dust depo of iron
111                jl = n_trc_indsbc(jpfer)
112                CALL iom_put( "Irondep", ( rf_trsfac(jl) * sf_trcsbc(jl)%fnow(:,:,1) / rn_sbc_time ) * 1.e+3 * tmask(:,:,1) )
113
114            ENDIF
115
116         ENDIF
117
118         ! Atmospheric input of PO4 dissolves in the water column
119         IF ( ln_trc_sbc(jppo4) ) THEN
120            DO_3D( 1, 1, 1, 1, 2, jpkm1 )
121               zdustdep = dust(ji,jj) * zwdust * rfact * EXP( -gdept(ji,jj,jk,Kmm) /( 250. * wdust ) )
122               !
123               tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) + zdustdep * 1.e-3 / mMass_P
124            END_3D
125         ENDIF
126
127         ! Atmospheric input of Si dissolves in the water column
128         IF ( ln_trc_sbc(jpsil) ) THEN
129            DO_3D( 1, 1, 1, 1, 2, jpkm1 )
130               zdustdep = dust(ji,jj) * zwdust * rfact * EXP( -gdept(ji,jj,jk,Kmm) /( 250. * wdust ) )
131               !
132               tr(ji,jj,jk,jpsil,Krhs) = tr(ji,jj,jk,jpsil,Krhs) + zdustdep * 0.269 / mMass_Si
133            END_3D
134         ENDIF
135
136         !
137         IF( lk_iomput ) THEN
138             ! dust concentration at surface
139             CALL iom_put( "pdust"  , dust(:,:) / ( wdust / rday ) * tmask(:,:,1) ) ! dust concentration at surface
140         ENDIF
141      ENDIF
142
143      ! -----------------------------------------
144      ! Add the external input of nutrients from river
145      ! ----------------------------------------------------------
146      IF( ll_river ) THEN
147          jl = n_trc_indcbc(jpno3)
148          DO_2D( 1, 1, 1, 1 )
149             DO jk = 1, nk_rnf(ji,jj)
150                zcoef = rn_rfact / ( e1e2t(ji,jj) * h_rnf(ji,jj) * rn_cbc_time ) * tmask(ji,jj,1)
151                zrivdin = rf_trcfac(jl) * sf_trccbc(jl)%fnow(ji,jj,1) * zcoef
152                tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) - rno3 * zrivdin * rfact
153            ENDDO
154          END_2D
155      ENDIF
156     
157      ! Add the external input of nutrients from nitrogen deposition
158      ! ----------------------------------------------------------
159      IF( ll_ndepo ) THEN
160         IF( ln_trc_sbc(jpno3) ) THEN
161            jl = n_trc_indsbc(jpno3)
162            DO_2D( 1, 1, 1, 1 )
163               zndep = rf_trsfac(jl) * sf_trcsbc(jl)%fnow(ji,jj,1) / e3t(ji,jj,1,Kmm) / rn_sbc_time
164               tr(ji,jj,1,jptal,Krhs) = tr(ji,jj,1,jptal,Krhs) - rno3 * zndep * rfact
165            END_2D
166         ENDIF
167         IF( ln_trc_sbc(jpnh4) ) THEN
168            jl = n_trc_indsbc(jpnh4)
169            DO_2D( 1, 1, 1, 1 )
170               zndep = rf_trsfac(jl) * sf_trcsbc(jl)%fnow(ji,jj,1) / e3t(ji,jj,1,Kmm) / rn_sbc_time
171               tr(ji,jj,1,jptal,Krhs) = tr(ji,jj,1,jptal,Krhs) + rno3 * zndep * rfact
172            END_2D
173         ENDIF
174      ENDIF
175      !
176      ! Iron input/uptake due to sea ice : Crude parameterization based on
177      ! Lancelot et al. Iron concentration in sea-ice is constant and set
178      ! in the namelist_pisces (icefeinput). ln_ironice is forced to false
179      ! when nn_ice_tr = 1
180      ! ----------------------------------------------------
181      IF( ln_ironice ) THEN
182         !
183         ! Compute the iron flux between sea ice and sea water
184         ! Simple parameterization assuming a fixed constant concentration in
185         ! sea-ice (icefeinput)
186         ! ------------------------------------------------------------------         
187         DO_2D( 1, 1, 1, 1 )
188            zdep     = rfact / e3t(ji,jj,1,Kmm)
189            zwflux   = fmmflx(ji,jj) / 1000._wp
190            zironice =  MAX( -0.99 * tr(ji,jj,1,jpfer,Kbb), -zwflux * icefeinput * zdep )
191            tr(ji,jj,1,jpfer,Krhs) = tr(ji,jj,1,jpfer,Krhs) + zironice
192         END_2D
193         !
194         ! iron flux from ice
195         IF( lk_iomput ) &
196         & CALL iom_put( "Ironice", MAX( -0.99 * tr(:,:,1,jpfer,Kbb), (-1.*fmmflx(:,:)/1000._wp )*icefeinput*1.e+3*tmask(:,:,1)) )
197         !
198      ENDIF
199
200      ! Add the external input of iron from sediment mobilization
201      ! ------------------------------------------------------
202      IF( ln_ironsed .AND. .NOT.lk_sed ) THEN
203          tr(:,:,:,jpfer,Krhs) = tr(:,:,:,jpfer,Krhs) + ironsed(:,:,:) * rfact
204          !
205          IF( lk_iomput )  CALL iom_put( "Ironsed", ironsed(:,:,:) * 1.e+3 * tmask(:,:,:) ) 
206      ENDIF
207
208      ! Add the external input of iron from hydrothermal vents
209      ! ------------------------------------------------------
210      IF( ln_hydrofe ) THEN
211         CALL fld_read( kt, 1, sf_hydrofe )
212         DO jk = 1, jpk
213            hydrofe(:,:,jk) = ( MAX( rtrn, sf_hydrofe(1)%fnow(:,:,jk) ) * hratio ) &
214              &              / ( e1e2t(:,:) * e3t(:,:,jk,Kmm) * ryyss + rtrn ) / 1000._wp &
215              &              * tmask(:,:,jk)
216         ENDDO
217                         tr(:,:,:,jpfer,Krhs) = tr(:,:,:,jpfer,Krhs) + hydrofe(:,:,:) * rfact
218         IF( ln_ligand ) tr(:,:,:,jplgw,Krhs) = tr(:,:,:,jplgw,Krhs) + ( hydrofe(:,:,:) * lgw_rath ) * rfact
219         !
220         IF( lk_iomput ) CALL iom_put( "HYDR", hydrofe(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! hydrothermal iron input
221      ENDIF
222      IF( ln_timing )  CALL timing_stop('p4z_bc')
223      !
224   END SUBROUTINE p4z_bc
225
226
227   SUBROUTINE p4z_bc_init( Kmm ) 
228      !!----------------------------------------------------------------------
229      !!                  ***  routine p4z_bc_init  ***
230      !!
231      !! ** purpose :   initialization of the external sources of nutrients
232      !!
233      !! ** method  :   read the files and compute the budget
234      !!                called at the first timestep (nittrc000)
235      !!
236      !! ** input   :   external netcdf files
237      !!
238      !!----------------------------------------------------------------------
239      INTEGER, INTENT( in ) ::   Kmm  ! time level index
240      INTEGER  :: ji, jj, jk, jm
241      INTEGER  :: ii0, ii1, ij0, ij1
242      INTEGER  :: numiron
243      INTEGER  :: ierr, ierr1, ierr2, ierr3
244      INTEGER  :: ios                 ! Local integer output status for namelist read
245      INTEGER  :: ik50                !  last level where depth less than 50 m
246      REAL(wp) :: zexpide, zdenitide, zmaskt, zsurfc, zsurfp,ze3t, ze3t2, zcslp
247      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zriver, zcmask
248      !
249      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
250      TYPE(FLD_N) ::   sn_dust, sn_ironsed, sn_hydrofe   ! informations about the fields to be read
251      !!
252      NAMELIST/nampisbc/cn_dir, sn_dust, sn_ironsed, sn_hydrofe, &
253        &                ln_ironsed, ln_ironice, ln_hydrofe, ln_dust_inp,   &
254        &                sedfeinput, distcoast, icefeinput, wdust, mfrac,   &
255        &                hratio, lgw_rath
256      !!----------------------------------------------------------------------
257      !
258      IF(lwp) THEN
259         WRITE(numout,*)
260         WRITE(numout,*) 'p4z_bc_init : initialization of the external sources of nutrients '
261         WRITE(numout,*) '~~~~~~~~~~~~ '
262      ENDIF
263      !                            !* set file information
264      READ  ( numnatp_ref, nampisbc, IOSTAT = ios, ERR = 901)
265901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampisbc in reference namelist' )
266      READ  ( numnatp_cfg, nampisbc, IOSTAT = ios, ERR = 902 )
267902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisbc in configuration namelist' )
268      IF(lwm) WRITE ( numonp, nampisbc )
269
270
271      IF(lwp) THEN
272         WRITE(numout,*) '   Namelist : nampissbc '
273         WRITE(numout,*) '      Fe input from sediments                  ln_ironsed  = ', ln_ironsed
274         WRITE(numout,*) '      Fe input from seaice                     ln_ironice  = ', ln_ironice
275         WRITE(numout,*) '      fe input from hydrothermal vents         ln_hydrofe  = ', ln_hydrofe
276         IF( ln_ironsed ) THEN
277            WRITE(numout,*) '      coastal release of iron                  sedfeinput  = ', sedfeinput
278            WRITE(numout,*) '      distance off the coast                   distcoast   = ', distcoast
279         ENDIF
280         IF( ln_ligand ) THEN
281            WRITE(numout,*) '      Weak ligand ratio from sed hydro sources  lgw_rath   = ', lgw_rath
282         ENDIF
283         IF( ln_ironice ) THEN
284            WRITE(numout,*) '      Iron concentration in sea ice            icefeinput  = ', icefeinput
285         ENDIF
286         IF( ln_trc_sbc(jpfer) ) THEN
287            WRITE(numout,*) '      Mineral Fe content of the dust           mfrac       = ', mfrac
288            WRITE(numout,*) '      sinking speed of the dust                wdust       = ', wdust
289         ENDIF
290         IF( ln_hydrofe ) THEN
291            WRITE(numout,*) '      Fe to 3He ratio assumed for vent iron supply hratio  = ', hratio
292         ENDIF
293      END IF
294
295      ll_bc    = ( ln_trcbc .AND. lltrcbc )  .OR. ln_hydrofe .OR. ln_ironsed .OR. ln_ironice .OR. ln_dust_inp
296      ll_dust  =  ln_dust_inp
297      ll_ndepo =  ln_trc_sbc(jpno3) .OR. ln_trc_sbc(jpnh4)   
298      ll_river =  ln_trc_cbc(jpno3) 
299
300      ! dust input from the atmosphere
301      ! ------------------------------
302      IF( ll_dust ) THEN
303         !
304         IF(lwp) WRITE(numout,*) '    initialize dust input from atmosphere '
305         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
306         !
307         ALLOCATE( dust(jpi,jpj) ) 
308         !
309         ALLOCATE( sf_dust(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst
310         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_bc_init: unable to allocate sf_dust structure' )
311         !
312         CALL fld_fill( sf_dust, (/ sn_dust /), cn_dir, 'p4z_bc_init', 'Atmospheric dust deposition', 'nampisbc' )
313                                   ALLOCATE( sf_dust(1)%fnow(jpi,jpj,1)   )
314         IF( sn_dust%ln_tint )     ALLOCATE( sf_dust(1)%fdta(jpi,jpj,1,2) )
315         !
316      END IF
317
318      ! coastal and island masks
319      ! ------------------------
320      IF( ln_ironsed ) THEN     
321         !
322         IF(lwp) WRITE(numout,*)
323         IF(lwp) WRITE(numout,*) '   ==>>>   ln_ironsed=T , computation of an island mask to enhance coastal supply of iron'
324         !
325         ALLOCATE( ironsed(jpi,jpj,jpk) )    ! allocation
326         !
327         CALL iom_open ( TRIM( sn_ironsed%clname ), numiron )
328         ALLOCATE( zcmask(jpi,jpj,jpk) )
329         CALL iom_get  ( numiron, jpdom_global, TRIM( sn_ironsed%clvar ), zcmask(:,:,:), 1 )
330         CALL iom_close( numiron )
331         !
332         ik50 = 5        !  last level where depth less than 50 m
333         DO jk = jpkm1, 1, -1
334            IF( gdept_1d(jk) > 50. )   ik50 = jk - 1
335         END DO
336         IF(lwp) WRITE(numout,*)
337         IF(lwp) WRITE(numout,*) ' Level corresponding to 50m depth ',  ik50,' ', gdept_1d(ik50+1)
338         DO_3D( 0, 0, 0, 0, 1, ik50 )
339            ze3t   = e3t_0(ji,jj,jk)
340            zsurfc =  e1u(ji,jj) * ( 1. - umask(ji  ,jj  ,jk) )   &
341                    + e1u(ji,jj) * ( 1. - umask(ji-1,jj  ,jk) )   &
342                    + e2v(ji,jj) * ( 1. - vmask(ji  ,jj  ,jk) )   &
343                    + e2v(ji,jj) * ( 1. - vmask(ji  ,jj-1,jk) )
344            zsurfp = zsurfc * ze3t / e1e2t(ji,jj)
345            ! estimation of the coastal slope : 5 km off the coast
346            ze3t2 = ze3t * ze3t
347            zcslp = SQRT( ( distcoast*distcoast + ze3t2 ) / ze3t2 )
348            !
349            zcmask(ji,jj,jk) = zcmask(ji,jj,jk) + zcslp * zsurfp
350         END_3D
351         !
352         CALL lbc_lnk( 'p4zbc', zcmask , 'T', 1.0_wp )      ! lateral boundary conditions on cmask   (sign unchanged)
353         !
354         DO_3D( 1, 1, 1, 1, 1, jpk )
355            zexpide   = MIN( 8.,( gdept(ji,jj,jk,Kmm) / 500. )**(-1.5) )
356            zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2
357            zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 )
358         END_3D
359         ! Coastal supply of iron
360         ! -------------------------
361         ironsed(:,:,jpk) = 0._wp
362         DO jk = 1, jpkm1
363            ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( e3t_0(:,:,jk) * rday )
364         END DO
365         DEALLOCATE( zcmask)
366      ENDIF
367      !
368      ! Iron from Hydrothermal vents
369      ! ------------------------
370      IF( ln_hydrofe ) THEN
371         !
372         IF(lwp) WRITE(numout,*)
373         IF(lwp) WRITE(numout,*) '   ==>>>   ln_hydrofe=T , Input of iron from hydrothermal vents'
374         !
375         ALLOCATE( hydrofe(jpi,jpj,jpk) )    ! allocation
376         !
377         ALLOCATE( sf_hydrofe(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst
378         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_bc_init: unable to allocate sf_hydro structure' )
379         !
380         CALL fld_fill( sf_hydrofe, (/ sn_hydrofe /), cn_dir, 'p4z_bc_init', 'Input of iron from hydrothermal vents', 'nampisbc' )
381                                   ALLOCATE( sf_hydrofe(1)%fnow(jpi,jpj,jpk)   )
382         IF( sn_hydrofe%ln_tint )    ALLOCATE( sf_hydrofe(1)%fdta(jpi,jpj,jpk,2) )
383         !
384      ENDIF
385      !
386   END SUBROUTINE p4z_bc_init
387
388   !!======================================================================
389END MODULE p4zbc
Note: See TracBrowser for help on using the repository browser.