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/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsed_medusa.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 5 years ago

The Dr Hook changes from my perl code.

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