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 @ 15449

Last change on this file since 15449 was 15449, checked in by cetlod, 2 years ago

dev_PISCO : minor bugfixes

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