source: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsed_medusa.F90 @ 8074

Last change on this file since 8074 was 8074, checked in by jpalmier, 4 years ago

JPALM — reverse MEDUSA cleaning and update MOCSY

File size: 19.0 KB
Line 
1MODULE trcsed_medusa
2   !!======================================================================
3   !!                         ***  MODULE trcsed_medusa  ***
4   !! TOP :   MEDUSA Compute loss of organic matter in the sediments
5   !!======================================================================
6   !! History :    -   !  1995-06 (M. Levy)  original code
7   !!              -   !  2000-12 (E. Kestenare)  clean up
8   !!             2.0  !  2007-12  (C. Deltel, G. Madec)  F90 + simplifications
9   !!              -   !  2008-08  (K. Popova) adaptation for MEDUSA
10   !!              -   !  2008-11  (A. Yool) continuing adaptation for MEDUSA
11   !!              -   !  2010-03  (A. Yool) updated for branch inclusion
12   !!              -   !  2011-04  (A. Yool) updated for ROAM project
13   !!----------------------------------------------------------------------
14#if defined key_medusa
15   !!----------------------------------------------------------------------
16   !!   'key_medusa'                                      MEDUSA bio-model
17   !!----------------------------------------------------------------------
18   !!   trc_sed_medusa        :  Compute loss of organic matter in the sediments
19   !!----------------------------------------------------------------------
20   USE oce_trc         !
21   USE trc
22   USE sms_medusa
23   !! AXY (10/02/09)
24   USE iom
25   !! USE trc_nam_dia                   ! JPALM 13-11-2015 -- if iom_use for diag
26   !! USE trc_nam_iom_medusa            ! JPALM 13-11-2015 -- if iom_use for diag
27   USE fldread                          !  time interpolation
28   USE lbclnk
29   USE prtctl_trc                       ! Print control for debbuging
30   !! JPALM (27-06-2016): add lk_oasis for CO2 and DMS coupling with atm
31   USE sbc_oce, ONLY: lk_oasis
32   USE oce,     ONLY: Dust_in_cpl
33   !! Check Dust dep
34# if defined key_debug_medusa   
35   !! USE trcrst, ONLY: trc_rst_dia_stat   !! variable stat
36# endif
37
38
39   IMPLICIT NONE
40   PRIVATE
41
42   PUBLIC   trc_sed_medusa     ! called in ???
43   PUBLIC   trc_sed_medusa_sbc     
44   PUBLIC   trc_sed_medusa_dust 
45
46   !! * Module variables
47   INTEGER ::                   &
48     ryyss,                     &  !: number of seconds per year
49     rmtss                         !: number of seconds per month
50
51   !! AXY (10/02/09)
52   LOGICAL, PUBLIC  ::   bdustfer  !: boolean for dust input from the atmosphere
53   REAL(wp), PUBLIC ::                 &
54      sedfeinput = 1.e-9_wp  ,         &
55      dustsolub  = 0.014_wp
56   REAL(wp), PARAMETER :: Fe_dust_mratio = 0.035   !! Fe:dust mass ratio = 0.035
57   INTEGER , PARAMETER ::        nbtimes = 365     !: maximum number of times record in a file
58   INTEGER  :: ntimes_dust                         ! number of time steps in a file
59
60   INTEGER ::                          &
61      numdust,                         &
62      nflx1,  nflx2,                   &
63      nflx11, nflx12
64   
65   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_dust      ! structure of input dust
66
67   
68   !!* Substitution
69#  include "domzgr_substitute.h90"
70   !!----------------------------------------------------------------------
71   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
72   !! $Id$
73   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
74   !!----------------------------------------------------------------------
75
76CONTAINS
77
78   SUBROUTINE trc_sed_medusa( kt )
79      !!---------------------------------------------------------------------
80      !!                     ***  ROUTINE trc_sed_medusa  ***
81      !!
82      !! ** Purpose :   compute the now trend due to the vertical sedimentation of
83      !!              detritus and add it to the general trend of detritus equations
84      !!
85      !! ** Method  :   this ROUTINE compute not exactly the advection but the
86      !!              transport term, i.e.  dz(wt) and dz(ws)., dz(wtr)
87      !!              using an upstream scheme
88      !!              the now vertical advection of tracers is given by:
89      !!                      dz(trn wn) = 1/bt dk+1( e1t e2t vsed (trn) )
90      !!              add this trend now to the general trend of tracer (ta,sa,tra):
91      !!                             tra = tra + dz(trn wn)
92      !!       
93      !!              IF 'key_trc_diabio' is defined, the now vertical advection
94      !!              trend of passive tracers is saved for futher diagnostics.
95      !!---------------------------------------------------------------------
96      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index     
97      !! AXY (10/02/09)
98      INTEGER  ::   jnt
99      !!
100      INTEGER  ::   ji, jj, jk
101      REAL(wp) ::   ztra
102      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwork
103
104      !! AXY (10/02/09)
105      REAL(wp) ::   rfact2
106
107      CHARACTER (len=25) :: charout
108     
109      !! JPALM - 26-11-2015 -add iom_use for diagnostic
110       REAL(wp), POINTER, DIMENSION(:,:  ) :: zw2d
111      !!---------------------------------------------------------------------
112      !!
113      IF( lk_iomput) THEN 
114           IF( med_diag%DSED%dgsave ) THEN
115               CALL wrk_alloc( jpi, jpj,      zw2d )
116                zw2d(:,:)      = 0.0      !!
117           ENDIF
118      ENDIF
119     
120      !! AXY (10/02/09)
121      jnt = 1
122      rfact2 = 1.0
123
124      ! Number of seconds per year and per month
125      ryyss = nyear_len(1) * rday
126      rmtss = ryyss / raamo
127
128      !! AXY (20/11/14): alter this to report on first MEDUSA call
129      IF( kt == nittrc000 ) THEN
130         IF(lwp) WRITE(numout,*)
131         IF(lwp) WRITE(numout,*) ' trc_sed_medusa: MEDUSA sedimentation'
132         IF(lwp) WRITE(numout,*) ' ~~~~~~~'
133    IF(lwp) WRITE(numout,*) ' kt =',kt
134      ENDIF
135
136      ! sedimentation of detrital nitrogen : upstream scheme
137      ! ----------------------------------------------------
138      !
139      zwork(:,:,:) = 0.e0        ! initialisation of sinking variable
140      ! for detrital nitrogen sedimentation only - jpdet
141      zwork(:,:,1  ) = 0.e0      ! surface value set to zero
142      !!    DO ji = 1, jpi
143      !!       zirondep(ji,jj,1) = (dustsolub * dust(ji,jj) / (55.85 * rmtss) + 3.e-10 / ryyss) &
144      !!       & * rfact2 / fse3t(ji,jj,1)
145      !!       zsidep  (ji,jj)   = 8.8 * 0.075 * dust(ji,jj) * rfact2 / &
146      !!       & (fse3t(ji,jj,1) * 28.1 * rmtss)
147      !!    END DO
148      !! END DO
149
150      ! sedimentation of detrital nitrogen : upstream scheme
151      ! ----------------------------------------------------
152      !
153      zwork(:,:,:) = 0.e0        ! initialisation of sinking variable
154      ! for detrital nitrogen sedimentation only - jpdet
155      zwork(:,:,1  ) = 0.e0      ! surface value set to zero
156      zwork(:,:,jpk) = 0.e0      ! bottom value  set to zero
157      !
158      ! tracer flux at w-point: we use -vsed (downward flux)  with simplification : no e1*e2
159      DO jk = 2, jpk
160         ! AXY (17/07/14): change "0.d0" to "0."
161         ! zwork(:,:,jk) = -vsed * max(trn(:,:,jk-1,jpdet),0.d0) * tmask(:,:,jk-1)
162         zwork(:,:,jk) = -vsed * max(trn(:,:,jk-1,jpdet),0.) * tmask(:,:,jk-1)
163         !
164         ! AXY (16/01/14): stop sinking in upper 10m to reduce model instability
165         !                 in shallower grid cells
166         ! if ( jk .lt. 9 ) zwork(:,:,jk) = 0.e0
167      END DO
168      !
169      ! tracer flux divergence at t-point added to the general trend
170      DO jk = 1, jpkm1
171         DO jj = 1, jpj
172            DO ji = 1,jpi
173               ztra  = - ( zwork(ji,jj,jk) - zwork(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
174               tra(ji,jj,jk,jpdet) = tra(ji,jj,jk,jpdet) + ztra
175# if defined key_trc_diabio
176               trbio(ji,jj,jk,8) = ztra
177# endif
178               IF (lk_iomput .AND. .NOT. ln_diatrc) THEN
179                     IF( med_diag%DSED%dgsave ) THEN
180                         zw2d(ji,jj) = zw2d(ji,jj) + ztra * fse3t(ji,jj,jk) * 86400.
181                      ENDIF   
182               ELSE IF( ln_diatrc )  THEN
183                    trc2d(ji,jj,8) = trc2d(ji,jj,8) + ztra * fse3t(ji,jj,jk) * 86400.
184               ENDIF   
185               
186            END DO
187         END DO
188      END DO
189      !
190# if defined key_trc_diabio
191      CALL lbc_lnk (trbio(:,:,1,8), 'T', 1. )                    ! Lateral boundary conditions on trcbio
192# endif
193      IF( ln_diatrc ) CALL lbc_lnk( trc2d(:,:,8), 'T', 1. )      ! Lateral boundary conditions on trc2d
194      !!
195      IF (lk_iomput .AND. .NOT. ln_diatrc) THEN
196           IF( med_diag%DSED%dgsave ) THEN
197                CALL iom_put( "DSED"  ,  zw2d)
198                CALL wrk_dealloc( jpi, jpj,    zw2d  )
199            ENDIF
200      ELSE IF (lk_iomput .AND. ln_diatrc)  THEN   
201          CALL iom_put( "DSED",trc2d(:,:,8) )
202      ENDIF
203      !!
204# if defined key_roam
205
206      ! sedimentation of detrital carbon : upstream scheme
207      ! --------------------------------------------------
208      !
209      zwork(:,:,:) = 0.e0        ! initialisation of sinking variable
210      ! for detrital carbon sedimentation only - jpdtc
211      zwork(:,:,1  ) = 0.e0      ! surface value set to zero
212      zwork(:,:,jpk) = 0.e0      ! bottom value  set to zero
213      !
214      ! tracer flux at w-point: we use -vsed (downward flux)  with simplification : no e1*e2
215      DO jk = 2, jpk
216         ! AXY (17/07/14): change "0.d0" to "0."
217         ! zwork(:,:,jk) = -vsed * max(trn(:,:,jk-1,jpdtc),0.d0) * tmask(:,:,jk-1)
218         zwork(:,:,jk) = -vsed * max(trn(:,:,jk-1,jpdtc),0.) * tmask(:,:,jk-1)
219         !
220         ! AXY (16/01/14): stop sinking in upper 10m to reduce model instability
221         !                 in shallower grid cells
222         ! if ( jk .lt. 9 ) zwork(:,:,jk) = 0.e0
223      END DO
224      !
225      ! tracer flux divergence at t-point added to the general trend
226      DO jk = 1, jpkm1
227         DO jj = 1, jpj
228            DO ji = 1,jpi
229               ztra  = - ( zwork(ji,jj,jk) - zwork(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
230               tra(ji,jj,jk,jpdtc) = tra(ji,jj,jk,jpdtc) + ztra
231!! #  if defined key_trc_diabio
232!!                trbio(ji,jj,jk,8) = ztra
233!! #  endif
234!!             IF( ln_diatrc ) &
235!!                &  trc2d(ji,jj,8) = trc2d(ji,jj,8) + ztra * fse3t(ji,jj,jk) * 86400.
236            END DO
237         END DO
238      END DO
239      !
240!! #  if defined key_trc_diabio
241!!       CALL lbc_lnk (trbio(:,:,1,8), 'T', 1. )                    ! Lateral boundary conditions on trcbio
242!! #  endif
243!!       IF( ln_diatrc ) CALL lbc_lnk( trc2d(:,:,8), 'T', 1. )      ! Lateral boundary conditions on trc2d
244!! #  if defined key_iomput
245!!       CALL iom_put( "DSED",trc2d(:,:,8) )
246!! #  endif
247
248# endif
249
250      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
251         WRITE(charout, FMT="('sed')")
252         CALL prt_ctl_trc_info(charout)
253         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
254      ENDIF
255
256   END SUBROUTINE trc_sed_medusa
257
258   !! ======================================================================
259   !! ======================================================================
260   !! ======================================================================
261
262   !! AXY (10/02/09)
263   !! JPALM -- 31-03-2016 -- Completely change trc_sed_medusa_sbc.
264   !!                     -- We now need to read dust file through a namelist.
265   !!                     To be able to use time varying dust depositions from
266   !!                     -- copy and adapt the PISCES p4z_sbc_ini subroutine
267   !!                     -- Only use the dust related part.     
268   SUBROUTINE trc_sed_medusa_sbc(kt)
269
270      !!----------------------------------------------------------------------
271      !!                  ***  ROUTINE trc_sed_medusa_sbc  ***
272      !!
273      !! ** Purpose :   Read and dust namelist and files.
274      !!                The interpolation is done in trc_sed through
275      !!                "CALL fld_read( kt, 1, sf_dust )"
276      !!
277      !! ** Method  :   Read the sbc namelist, and the adapted dust file, if required
278      !!                called at the first timestep (nittrc000)
279      !!
280      !! ** input   :   -- namelist sbc ref and cfg
281      !!                -- external netcdf files
282      !!
283      !!----------------------------------------------------------------------
284      !! * arguments
285      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
286      INTEGER  :: ji, jj, jk, jm, ifpr
287      INTEGER  :: ii0, ii1, ij0, ij1
288      INTEGER  :: numdust
289      INTEGER  :: ierr 
290      INTEGER  :: jfld                ! dummy loop arguments
291      INTEGER  :: ios                 ! Local integer output status for namelist read
292      INTEGER  :: isrow               ! index for ORCA1 starting row
293      REAL(wp) :: ztimes_dust
294      REAL(wp), DIMENSION(nbtimes) :: zsteps                 ! times records
295      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zdust
296      !
297      CHARACTER(len=100)         ::  cn_dir       ! Root directory for location of ssr files
298      TYPE(FLD_N), DIMENSION(1)  ::   slf_d       ! array of namelist informations on the fields to read
299      TYPE(FLD_N)                ::   sn_dust     ! informations about the fields to be read
300      !
301      NAMELIST/nammedsbc/cn_dir, sn_dust, bdustfer 
302
303      !!---------------------------------------------------------------------
304      !
305      IF( nn_timing == 1 )  CALL timing_start('trc_sed_medusa_sbc')
306      !
307      !                            !* set file information
308      REWIND( numnatp_ref )        ! Namelist nammedsbc in reference namelist : MEDUSA external sources of Dust
309      READ  ( numnatp_ref, nammedsbc, IOSTAT = ios, ERR = 901)
310901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammedsbc in reference namelist', lwp )
311
312      REWIND( numnatp_cfg )        ! Namelist nammedsbc in configuration namelist : MEDUSA external sources of Dust
313      READ  ( numnatp_cfg, nammedsbc, IOSTAT = ios, ERR = 902 )
314902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammedsbc in configuration namelist', lwp )
315      IF(lwm) WRITE ( numonp, nammedsbc )
316
317      IF(lwp) THEN
318          WRITE(numout,*) ' '
319          WRITE(numout,*) ' namelist : nammedsbc '
320          WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ '
321          WRITE(numout,*) '    dust input from the atmosphere bdustfer     = ', bdustfer
322      END IF
323
324      ! dust input from the atmosphere
325      ! ------------------------------
326      IF( bdustfer ) THEN
327         !
328         IF(lwp) WRITE(numout,*) '    initialize dust input from atmosphere '
329         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
330         !
331         !! already allocated in sms_medusa
332         !!ALLOCATE( dust(jpi,jpj) )    ! allocation
333         !
334         slf_d(1) = sn_dust                          ! store namelist information in an array
335         !
336         ALLOCATE( sf_dust(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst
337         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'trc_sed_medusa_sbc: unable to allocate sf_dust structure' )
338         ALLOCATE( sf_dust(1)%fnow(jpi,jpj,1))
339         IF( slf_d(1)%ln_tint )     ALLOCATE( sf_dust(1)%fdta(jpi,jpj,1,2) )
340         !
341         CALL fld_fill( sf_dust, slf_d, cn_dir, 'trc_sed_medusa_sbc', 'Atmospheric dust deposition', 'nammedsed' )
342         !
343         CALL fld_read( kt, 1, sf_dust )
344         dust(:,:) = sf_dust(1)%fnow(:,:,1)
345         !
346      ELSEIF (lk_oasis) THEN
347         dust = Dust_in_cpl
348      ELSE
349         dust(:,:) = 0.0
350      END IF
351      !
352      zirondep(:,:) = 0.e0     !! Initialisation of deposition variables
353      zirondep(:,:) = dust(:,:) * Fe_dust_mratio / xfe_mass * 1.e6 * 86400.      !! mmol-Fe/m2/d
354      !
355      IF( nn_timing == 1 )  CALL timing_stop('trc_sed_medusa_sbc')
356      !
357   END SUBROUTINE trc_sed_medusa_sbc
358
359   !! ======================================================================
360   !! ======================================================================
361   !! ======================================================================
362
363   !! AXY & JPALM (28/02/17)
364
365   SUBROUTINE trc_sed_medusa_dust( kt )
366      !!---------------------------------------------------------------------
367      !!                     ***  ROUTINE trc_sed_medusa_dust  ***
368      !!
369      !! ** Purpose : compute current dust *before* trc_bio_medusa call
370      !!
371      !! ** Method  : does what it says on the tin
372      !!---------------------------------------------------------------------
373      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index     
374
375      !! AXY (20/11/14): alter this to report on first MEDUSA call
376      IF( kt == nittrc000 ) THEN
377         IF(lwp) WRITE(numout,*)
378         IF(lwp) WRITE(numout,*) ' trc_sed_medusa_dust: MEDUSA dust timestep'
379         IF(lwp) WRITE(numout,*) ' ~~~~~~~'
380    IF(lwp) WRITE(numout,*) ' kt =',kt
381      ENDIF
382
383      !! AXY (04/11/13): replace this with a call in trc_ini_medusa
384      !! AXY (25/02/10)
385      !! call routine for populating CCD array if this is the first time-step
386      !! IF( kt == nittrc000 ) CALL medusa_ccd( kt )
387
388      !! AXY (04/11/13): replace this with a call in trc_ini_medusa
389      !! AXY (26/01/12)
390      !! call routine for populating river arrays if this is the first time-step
391      !! IF( kt == nittrc000 ) CALL medusa_river( kt )
392
393      !! AXY (10/02/09)
394      !! IF( (jnt == 1) .and. (bdustfer) )  CALL trc_sed_medusa_sbc( kt )
395
396      !! JPALM -- 31-03-2016 -- rewrite trc_sed_medusa_sbc.
397      !! IF (kt == nittrc000 ) CALL trc_sed_medusa_sbc
398
399      !! JPALM -- 20-07-2016 -- adapt dust forcing fields reading and conversion
400      !!                     To read dust dep in kg-dust/m2/s instead of g-Fe/m2/month
401      !!                     So all forcings and coupling dust dep are in the same SI units
402      !!                     and then convert in mmol-Fe/m2/day
403
404      IF( bdustfer ) THEN
405            CALL fld_read( kt, 1, sf_dust )
406            dust(:,:) = sf_dust(1)%fnow(:,:,1)
407      ELSEIF (lk_oasis) THEN
408         dust = Dust_in_cpl
409      ELSE
410         dust(:,:) = 0.0
411      ENDIF
412      !!
413      zirondep(:,:) = 0.e0     !! Initialisation of deposition variables
414      zirondep(:,:) = dust(:,:) * Fe_dust_mratio / xfe_mass * 1.e6 * 86400.  !! mmol-Fe/m2/d
415     
416      !! JPALM -- 20-07-2016 -- Zirondep and zsidep are not used.
417      !!                     So comment out the following lines. but keep them
418      !!                     as we may want to used them later on
419      !!================================================     
420      !!
421      !! zirondep(:,:,:) = 0.e0     !! Initialisation of deposition variables
422      !! zsidep  (:,:)   = 0.e0
423      !!
424      !! Iron and Si deposition at the surface
425      !! -------------------------------------
426      !!
427      !! DO jj = 1, jpj
428      !!    DO ji = 1, jpi
429      !!       zirondep(ji,jj,1) = (dustsolub * dust(ji,jj) / (55.85 * rmtss) + 3.e-10 / ryyss) &
430      !!       & * rfact2 / fse3t(ji,jj,1)
431      !!       zsidep  (ji,jj)   = 8.8 * 0.075 * dust(ji,jj) * rfact2 / &
432      !!       & (fse3t(ji,jj,1) * 28.1 * rmtss)
433      !!    END DO
434      !! END DO
435
436   END SUBROUTINE trc_sed_medusa_dust
437
438#else
439   !!======================================================================
440   !!  Dummy module :                                   No MEDUSA bio-model
441   !!======================================================================
442CONTAINS
443   SUBROUTINE trc_sed_medusa( kt )                   ! Empty routine
444      INTEGER, INTENT( in ) ::   kt
445      WRITE(*,*) 'trc_sed_medusa: You should not have seen this print! error?', kt
446   END SUBROUTINE trc_sed_medusa
447#endif 
448
449   !!======================================================================
450END MODULE  trcsed_medusa
451
Note: See TracBrowser for help on using the repository browser.