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.
trcsed_medusa.F90 in branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsed_medusa.F90 @ 10045

Last change on this file since 10045 was 10045, checked in by jpalmier, 6 years ago

Andrew's changes to add the OMIP double_DIC (activated with key_omip_dic)

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