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/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z/p4zbc.F90 @ 12258

Last change on this file since 12258 was 12258, checked in by cetlod, 4 years ago

dev_r11943_MERGE_2019 : : Cleaning PISCES diagnostics output ; fully sette tested. run.stat and tracer.stat are identical to the ones of dev_r11943_MERGE_2019
File size: 16.8 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 "vectopt_loop_substitute.h90"
50   !!----------------------------------------------------------------------
51   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
52   !! $Id: p4zbc.F90 10869 2019-04-15 10:34:03Z cetlod $
53   !! Software governed by the CeCILL license (see ./LICENSE)
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   SUBROUTINE p4z_bc( kt, Kbb, Kmm, Krhs )
58      !!----------------------------------------------------------------------
59      !!                  ***  routine p4z_bc  ***
60      !!
61      !! ** purpose :   read and interpolate the external sources of nutrients
62      !!
63      !! ** method  :   read the files and interpolate the appropriate variables
64      !!
65      !! ** input   :   external netcdf files
66      !!
67      !!----------------------------------------------------------------------
68      INTEGER, INTENT(in) ::   kt              ! ocean time step
69      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs  ! time level index
70      !
71      INTEGER  ::  ji, jj, jk, jl 
72      REAL(wp) ::  zcoef, zyyss
73      REAL(wp) ::  zdep, ztrfer, zwdust, zwflux, zrivdin
74      !
75      CHARACTER (len=25) :: charout
76      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zirondep
77      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zironice, zndep
78      !!---------------------------------------------------------------------
79      !
80      IF( ln_timing )   CALL timing_start('p4z_bc')
81      !
82      IF( ll_dust )  THEN
83         ALLOCATE(  zirondep(jpi,jpj,jpk) )
84         !
85         CALL fld_read( kt, 1, sf_dust )
86         dust(:,:) = MAX( rtrn, sf_dust(1)%fnow(:,:,1) )
87         !
88         jl = n_trc_indsbc(jpfer)
89         zirondep(:,:,1) = rf_trsfac(jl) * sf_trcsbc(jl)%fnow(:,:,1) / e3t(:,:,1,Kmm) / rn_sbc_time
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/j
92         zwdust = 0.03 / ( wdust / rday ) / ( 270. * rday )
93         DO jk = 2, jpkm1
94            zirondep(:,:,jk) = ( mfrac * dust(:,:) * zwdust / mMass_Fe ) * rfact * EXP( -gdept(:,:,jk,Kmm) / 540. )
95            tr(:,:,jk,jpfer,Krhs) = tr(:,:,jk,jpfer,Krhs) + zirondep(:,:,jk)
96            tr(:,:,jk,jppo4,Krhs) = tr(:,:,jk,jppo4,Krhs) + zirondep(:,:,jk) * 0.023
97         ENDDO
98         !
99         IF( lk_iomput ) THEN
100             CALL iom_put( "Irondep", zirondep(:,:,1) * 1.e+3 * rfactr * e3t(:,:,1,Kmm) * tmask(:,:,1) ) ! surface downward dust depo of iron
101             CALL iom_put( "pdust"  , dust(:,:) / ( wdust * rday ) * tmask(:,:,1) ) ! dust concentration at surface
102         ENDIF
103         DEALLOCATE( zirondep )
104      ENDIF
105
106      ! N/P and Si releases due to coastal rivers
107      ! Compute river at nit000 or only if there is more than 1 time record in river file
108      ! -----------------------------------------
109            ! Add the external input of nutrients from river
110      ! ----------------------------------------------------------
111      IF( ll_river ) THEN
112          jl = n_trc_indcbc(jpno3)
113          DO jj = 1, jpj
114             DO ji = 1, jpi
115                DO jk = 1, nk_rnf(ji,jj)
116                   zcoef = rn_rfact / ( e1e2t(ji,jj) * h_rnf(ji,jj) * rn_cbc_time ) * tmask(ji,jj,1)
117                   zrivdin = rf_trcfac(jl) * sf_trccbc(jl)%fnow(ji,jj,1) * zcoef
118                   tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) - rno3 * zrivdin * rfact
119               ENDDO
120             END DO
121          END DO
122      ENDIF
123     
124      ! Add the external input of nutrients from nitrogen deposition
125      ! ----------------------------------------------------------
126      IF( ll_ndepo ) THEN
127         ALLOCATE( zndep(jpi,jpj) )
128         IF( ln_trc_sbc(jpno3) ) THEN
129            jl = n_trc_indsbc(jpno3)
130            zndep(:,:) = rf_trsfac(jl) * sf_trcsbc(jl)%fnow(:,:,1) / e3t(:,:,1,Kmm) / rn_sbc_time
131            tr(:,:,1,jptal,Krhs) = tr(:,:,1,jptal,Krhs) - rno3 * zndep(:,:) * rfact
132         ENDIF
133         IF( ln_trc_sbc(jpnh4) ) THEN
134            jl = n_trc_indsbc(jpnh4)
135            zndep(:,:) = rf_trsfac(jl) * sf_trcsbc(jl)%fnow(:,:,1) / e3t(:,:,1,Kmm) / rn_sbc_time
136            tr(:,:,1,jptal,Krhs) = tr(:,:,1,jptal,Krhs) - rno3 * zndep(:,:) * rfact
137         ENDIF
138         DEALLOCATE( zndep )
139      ENDIF
140      !
141      ! Iron input/uptake due to sea ice : Crude parameterization based on
142      ! Lancelot et al.
143      ! ----------------------------------------------------
144      IF( ln_ironice ) THEN
145         !
146         ALLOCATE( zironice(jpi,jpj) )
147         !
148         DO jj = 1, jpj
149            DO ji = 1, jpi
150               zdep    = rfact / e3t(ji,jj,1,Kmm)
151               zwflux  = fmmflx(ji,jj) / 1000._wp
152               zironice(ji,jj) =  MAX( -0.99 * tr(ji,jj,1,jpfer,Kbb), -zwflux * icefeinput * zdep )
153            END DO
154         END DO
155         !
156         tr(:,:,1,jpfer,Krhs) = tr(:,:,1,jpfer,Krhs) + zironice(:,:)
157         !
158         IF( lk_iomput )  CALL iom_put( "Ironice", zironice(:,:) * 1.e+3 * rfactr * e3t(:,:,1,Kmm) * tmask(:,:,1) ) ! iron flux from ice
159         !
160         DEALLOCATE( zironice )
161         !
162      ENDIF
163
164      ! Add the external input of iron from sediment mobilization
165      ! ------------------------------------------------------
166      IF( ln_ironsed .AND. .NOT.lk_sed ) THEN
167          tr(:,:,:,jpfer,Krhs) = tr(:,:,:,jpfer,Krhs) + ironsed(:,:,:) * rfact
168          !
169          IF( lk_iomput )  CALL iom_put( "Ironsed", ironsed(:,:,:) * 1.e+3 * tmask(:,:,:) ) 
170      ENDIF
171
172      ! Add the external input of iron from hydrothermal vents
173      ! ------------------------------------------------------
174      IF( ln_hydrofe ) THEN
175         CALL fld_read( kt, 1, sf_hydrofe )
176         DO jk = 1, jpk
177            hydrofe(:,:,jk) = ( MAX( rtrn, sf_hydrofe(1)%fnow(:,:,jk) ) * hratio ) &
178              &              / ( e1e2t(:,:) * e3t(:,:,jk,Kmm) * ryyss + rtrn ) / 1000._wp &
179              &              * tmask(:,:,jk)
180         ENDDO
181                         tr(:,:,:,jpfer,Krhs) = tr(:,:,:,jpfer,Krhs) + hydrofe(:,:,:) * rfact
182         IF( ln_ligand ) tr(:,:,:,jplgw,Krhs) = tr(:,:,:,jplgw,Krhs) + ( hydrofe(:,:,:) * lgw_rath ) * rfact
183         !
184         IF( lk_iomput ) CALL iom_put( "HYDR", hydrofe(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! hydrothermal iron input
185      ENDIF
186      IF( ln_timing )  CALL timing_stop('p4z_bc')
187      !
188   END SUBROUTINE p4z_bc
189
190
191   SUBROUTINE p4z_bc_init( Kmm ) 
192      !!----------------------------------------------------------------------
193      !!                  ***  routine p4z_bc_init  ***
194      !!
195      !! ** purpose :   initialization of the external sources of nutrients
196      !!
197      !! ** method  :   read the files and compute the budget
198      !!                called at the first timestep (nittrc000)
199      !!
200      !! ** input   :   external netcdf files
201      !!
202      !!----------------------------------------------------------------------
203      INTEGER, INTENT( in ) ::   Kmm  ! time level index
204      INTEGER  :: ji, jj, jk, jm
205      INTEGER  :: ii0, ii1, ij0, ij1
206      INTEGER  :: numiron
207      INTEGER  :: ierr, ierr1, ierr2, ierr3
208      INTEGER  :: ios                 ! Local integer output status for namelist read
209      INTEGER  :: ik50                !  last level where depth less than 50 m
210      REAL(wp) :: zexpide, zdenitide, zmaskt, zsurfc, zsurfp,ze3t, ze3t2, zcslp
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) ::   sn_dust, sn_ironsed, sn_hydrofe   ! informations about the fields to be read
215      !!
216      NAMELIST/nampisbc/cn_dir, sn_dust, sn_ironsed, sn_hydrofe, &
217        &                ln_ironsed, ln_ironice, ln_hydrofe,    &
218        &                sedfeinput, distcoast, icefeinput, wdust, mfrac,  &
219        &                hratio, lgw_rath
220      !!----------------------------------------------------------------------
221      !
222      IF(lwp) THEN
223         WRITE(numout,*)
224         WRITE(numout,*) 'p4z_bc_init : initialization of the external sources of nutrients '
225         WRITE(numout,*) '~~~~~~~~~~~~ '
226      ENDIF
227      !                            !* set file information
228      READ  ( numnatp_ref, nampisbc, IOSTAT = ios, ERR = 901)
229901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampisbc in reference namelist' )
230      READ  ( numnatp_cfg, nampisbc, IOSTAT = ios, ERR = 902 )
231902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisbc in configuration namelist' )
232      IF(lwm) WRITE ( numonp, nampisbc )
233
234
235      IF(lwp) THEN
236         WRITE(numout,*) '   Namelist : nampissbc '
237         WRITE(numout,*) '      Fe input from sediments                  ln_ironsed  = ', ln_ironsed
238         WRITE(numout,*) '      Fe input from seaice                     ln_ironice  = ', ln_ironice
239         WRITE(numout,*) '      fe input from hydrothermal vents         ln_hydrofe  = ', ln_hydrofe
240         IF( ln_ironsed ) THEN
241            WRITE(numout,*) '      coastal release of iron                  sedfeinput  = ', sedfeinput
242            WRITE(numout,*) '      distance off the coast                   distcoast   = ', distcoast
243         ENDIF
244         IF( ln_ligand ) THEN
245            WRITE(numout,*) '      Weak ligand ratio from sed hydro sources  lgw_rath   = ', lgw_rath
246         ENDIF
247         IF( ln_ironice ) THEN
248            WRITE(numout,*) '      Iron concentration in sea ice            icefeinput  = ', icefeinput
249         ENDIF
250         IF( ln_trc_sbc(jpfer) ) THEN
251            WRITE(numout,*) '      Mineral Fe content of the dust           mfrac       = ', mfrac
252            WRITE(numout,*) '      sinking speed of the dust                wdust       = ', wdust
253         ENDIF
254         IF( ln_hydrofe ) THEN
255            WRITE(numout,*) '      Fe to 3He ratio assumed for vent iron supply hratio  = ', hratio
256         ENDIF
257      END IF
258
259      ll_bc    = ( ln_trcbc .AND. lltrcbc )  .OR. ln_hydrofe .OR. ln_ironsed .OR. ln_ironice
260      ll_dust  =  ln_trc_sbc(jpfer)   
261      ll_ndepo =  ln_trc_sbc(jpno3) .OR. ln_trc_sbc(jpnh4)   
262      ll_river =  ln_trc_cbc(jpno3) 
263
264      ! dust input from the atmosphere
265      ! ------------------------------
266      IF( ll_dust ) THEN
267         !
268         IF(lwp) WRITE(numout,*) '    initialize dust input from atmosphere '
269         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
270         !
271         ALLOCATE( dust(jpi,jpj) ) 
272         !
273         ALLOCATE( sf_dust(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst
274         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_dust structure' )
275         !
276         CALL fld_fill( sf_dust, (/ sn_dust /), cn_dir, 'p4z_sed_init', 'Atmospheric dust deposition', 'nampissed' )
277                                   ALLOCATE( sf_dust(1)%fnow(jpi,jpj,1)   )
278         IF( sn_dust%ln_tint )     ALLOCATE( sf_dust(1)%fdta(jpi,jpj,1,2) )
279         !
280      END IF
281
282      ! coastal and island masks
283      ! ------------------------
284      IF( ln_ironsed ) THEN     
285         !
286         IF(lwp) WRITE(numout,*)
287         IF(lwp) WRITE(numout,*) '   ==>>>   ln_ironsed=T , computation of an island mask to enhance coastal supply of iron'
288         !
289         ALLOCATE( ironsed(jpi,jpj,jpk) )    ! allocation
290         !
291         CALL iom_open ( TRIM( sn_ironsed%clname ), numiron )
292         ALLOCATE( zcmask(jpi,jpj,jpk) )
293         CALL iom_get  ( numiron, jpdom_data, TRIM( sn_ironsed%clvar ), zcmask(:,:,:), 1 )
294         CALL iom_close( numiron )
295         !
296         ik50 = 5        !  last level where depth less than 50 m
297         DO jk = jpkm1, 1, -1
298            IF( gdept_1d(jk) > 50. )   ik50 = jk - 1
299         END DO
300         IF(lwp) WRITE(numout,*)
301         IF(lwp) WRITE(numout,*) ' Level corresponding to 50m depth ',  ik50,' ', gdept_1d(ik50+1)
302         DO jk = 1, ik50
303            DO jj = 2, jpjm1
304               DO ji = fs_2, fs_jpim1
305                  ze3t   = e3t_0(ji,jj,jk)
306                  zsurfc =  e1u(ji,jj) * ( 1. - umask(ji  ,jj  ,jk) )   &
307                          + e1u(ji,jj) * ( 1. - umask(ji-1,jj  ,jk) )   &
308                          + e2v(ji,jj) * ( 1. - vmask(ji  ,jj  ,jk) )   &
309                          + e2v(ji,jj) * ( 1. - vmask(ji  ,jj-1,jk) )
310                  zsurfp = zsurfc * ze3t / e1e2t(ji,jj)
311                  ! estimation of the coastal slope : 5 km off the coast
312                  ze3t2 = ze3t * ze3t
313                  zcslp = SQRT( ( distcoast*distcoast + ze3t2 ) / ze3t2 )
314                  !
315                  zcmask(ji,jj,jk) = zcmask(ji,jj,jk) + zcslp * zsurfp
316               END DO
317            END DO
318         END DO
319         !
320         CALL lbc_lnk( 'p4zbc', zcmask , 'T', 1. )      ! lateral boundary conditions on cmask   (sign unchanged)
321         !
322         DO jk = 1, jpk
323            DO jj = 1, jpj
324               DO ji = 1, jpi
325                  zexpide   = MIN( 8.,( gdept(ji,jj,jk,Kmm) / 500. )**(-1.5) )
326                  zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2
327                  zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 )
328               END DO
329            END DO
330         END DO
331         ! Coastal supply of iron
332         ! -------------------------
333         ironsed(:,:,jpk) = 0._wp
334         DO jk = 1, jpkm1
335            ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( e3t_0(:,:,jk) * rday )
336         END DO
337         DEALLOCATE( zcmask)
338      ENDIF
339      !
340      ! Iron from Hydrothermal vents
341      ! ------------------------
342      IF( ln_hydrofe ) THEN
343         !
344         IF(lwp) WRITE(numout,*)
345         IF(lwp) WRITE(numout,*) '   ==>>>   ln_hydrofe=T , Input of iron from hydrothermal vents'
346         !
347         ALLOCATE( hydrofe(jpi,jpj,jpk) )    ! allocation
348         !
349         ALLOCATE( sf_hydrofe(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst
350         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_init: unable to allocate sf_hydro structure' )
351         !
352         CALL fld_fill( sf_hydrofe, (/ sn_hydrofe /), cn_dir, 'p4z_sed_init', 'Input of iron from hydrothermal vents', 'nampisbc' )
353                                   ALLOCATE( sf_hydrofe(1)%fnow(jpi,jpj,jpk)   )
354         IF( sn_hydrofe%ln_tint )    ALLOCATE( sf_hydrofe(1)%fdta(jpi,jpj,jpk,2) )
355         !
356      ENDIF
357      !
358   END SUBROUTINE p4z_bc_init
359
360   !!======================================================================
361END MODULE p4zbc
Note: See TracBrowser for help on using the repository browser.