source: branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsed_medusa.F90 @ 9309

Last change on this file since 9309 was 8442, checked in by frrh, 3 years ago

Commit changes relating to Met Office GMED ticket 340 for the
tidying of MEDUSA related code and debugging statements in the TOP code.

Only code introduced at revision 8434 of branch
http://fcm3/projects/NEMO.xm/log/branches/NERC/dev_r5518_GO6_split_trcbiomedusa
is included here, all previous revisions of that branch having been dealt with
under GMED ticket 339.

File size: 17.7 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( med_diag%DSED%dgsave ) THEN
176                   zw2d(ji,jj) = zw2d(ji,jj) + ztra * fse3t(ji,jj,jk) * 86400.
177               ENDIF   
178               
179            END DO
180         END DO
181      END DO
182      !
183      IF( med_diag%DSED%dgsave ) THEN
184           CALL iom_put( "DSED"  ,  zw2d)
185           CALL wrk_dealloc( jpi, jpj,    zw2d  )
186      ENDIF
187      !!
188# if defined key_roam
189
190      ! sedimentation of detrital carbon : upstream scheme
191      ! --------------------------------------------------
192      !
193      zwork(:,:,:) = 0.e0        ! initialisation of sinking variable
194      ! for detrital carbon sedimentation only - jpdtc
195      zwork(:,:,1  ) = 0.e0      ! surface value set to zero
196      zwork(:,:,jpk) = 0.e0      ! bottom value  set to zero
197      !
198      ! tracer flux at w-point: we use -vsed (downward flux)  with simplification : no e1*e2
199      DO jk = 2, jpk
200         ! AXY (17/07/14): change "0.d0" to "0."
201         ! zwork(:,:,jk) = -vsed * max(trn(:,:,jk-1,jpdtc),0.d0) * tmask(:,:,jk-1)
202         zwork(:,:,jk) = -vsed * max(trn(:,:,jk-1,jpdtc),0.) * tmask(:,:,jk-1)
203         !
204         ! AXY (16/01/14): stop sinking in upper 10m to reduce model instability
205         !                 in shallower grid cells
206         ! if ( jk .lt. 9 ) zwork(:,:,jk) = 0.e0
207      END DO
208      !
209      ! tracer flux divergence at t-point added to the general trend
210      DO jk = 1, jpkm1
211         DO jj = 1, jpj
212            DO ji = 1,jpi
213               ztra  = - ( zwork(ji,jj,jk) - zwork(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
214               tra(ji,jj,jk,jpdtc) = tra(ji,jj,jk,jpdtc) + ztra
215            END DO
216         END DO
217      END DO
218      !
219
220# endif
221
222      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
223         WRITE(charout, FMT="('sed')")
224         CALL prt_ctl_trc_info(charout)
225         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
226      ENDIF
227
228   END SUBROUTINE trc_sed_medusa
229
230   !! ======================================================================
231   !! ======================================================================
232   !! ======================================================================
233
234   !! AXY (10/02/09)
235   !! JPALM -- 31-03-2016 -- Completely change trc_sed_medusa_sbc.
236   !!                     -- We now need to read dust file through a namelist.
237   !!                     To be able to use time varying dust depositions from
238   !!                     -- copy and adapt the PISCES p4z_sbc_ini subroutine
239   !!                     -- Only use the dust related part.     
240   SUBROUTINE trc_sed_medusa_sbc(kt)
241
242      !!----------------------------------------------------------------------
243      !!                  ***  ROUTINE trc_sed_medusa_sbc  ***
244      !!
245      !! ** Purpose :   Read and dust namelist and files.
246      !!                The interpolation is done in trc_sed through
247      !!                "CALL fld_read( kt, 1, sf_dust )"
248      !!
249      !! ** Method  :   Read the sbc namelist, and the adapted dust file, if required
250      !!                called at the first timestep (nittrc000)
251      !!
252      !! ** input   :   -- namelist sbc ref and cfg
253      !!                -- external netcdf files
254      !!
255      !!----------------------------------------------------------------------
256      !! * arguments
257      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
258      INTEGER  :: ji, jj, jk, jm, ifpr
259      INTEGER  :: ii0, ii1, ij0, ij1
260      INTEGER  :: numdust
261      INTEGER  :: ierr 
262      INTEGER  :: jfld                ! dummy loop arguments
263      INTEGER  :: ios                 ! Local integer output status for namelist read
264      INTEGER  :: isrow               ! index for ORCA1 starting row
265      REAL(wp) :: ztimes_dust
266      REAL(wp), DIMENSION(nbtimes) :: zsteps                 ! times records
267      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zdust
268      !
269      CHARACTER(len=100)         ::  cn_dir       ! Root directory for location of ssr files
270      TYPE(FLD_N), DIMENSION(1)  ::   slf_d       ! array of namelist informations on the fields to read
271      TYPE(FLD_N)                ::   sn_dust     ! informations about the fields to be read
272      !
273      NAMELIST/nammedsbc/cn_dir, sn_dust, bdustfer 
274
275      !!---------------------------------------------------------------------
276      !
277      IF( nn_timing == 1 )  CALL timing_start('trc_sed_medusa_sbc')
278      !
279      !                            !* set file information
280      REWIND( numnatp_ref )        ! Namelist nammedsbc in reference namelist : MEDUSA external sources of Dust
281      READ  ( numnatp_ref, nammedsbc, IOSTAT = ios, ERR = 901)
282901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammedsbc in reference namelist', lwp )
283
284      REWIND( numnatp_cfg )        ! Namelist nammedsbc in configuration namelist : MEDUSA external sources of Dust
285      READ  ( numnatp_cfg, nammedsbc, IOSTAT = ios, ERR = 902 )
286902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammedsbc in configuration namelist', lwp )
287      IF(lwm) WRITE ( numonp, nammedsbc )
288
289      IF(lwp) THEN
290          WRITE(numout,*) ' '
291          WRITE(numout,*) ' namelist : nammedsbc '
292          WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ '
293          WRITE(numout,*) '    dust input from the atmosphere bdustfer     = ', bdustfer
294      END IF
295
296      ! dust input from the atmosphere
297      ! ------------------------------
298      IF( bdustfer ) THEN
299         !
300         IF(lwp) WRITE(numout,*) '    initialize dust input from atmosphere '
301         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
302         !
303         !! already allocated in sms_medusa
304         !!ALLOCATE( dust(jpi,jpj) )    ! allocation
305         !
306         slf_d(1) = sn_dust                          ! store namelist information in an array
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', 'trc_sed_medusa_sbc: unable to allocate sf_dust structure' )
310         ALLOCATE( sf_dust(1)%fnow(jpi,jpj,1))
311         IF( slf_d(1)%ln_tint )     ALLOCATE( sf_dust(1)%fdta(jpi,jpj,1,2) )
312         !
313         CALL fld_fill( sf_dust, slf_d, cn_dir, 'trc_sed_medusa_sbc', 'Atmospheric dust deposition', 'nammedsed' )
314         !
315         CALL fld_read( kt, 1, sf_dust )
316         dust(:,:) = sf_dust(1)%fnow(:,:,1)
317         !
318      ELSEIF (lk_oasis) THEN
319         dust = Dust_in_cpl
320      ELSE
321         dust(:,:) = 0.0
322      END IF
323      !
324      zirondep(:,:) = 0.e0     !! Initialisation of deposition variables
325      zirondep(:,:) = dust(:,:) * Fe_dust_mratio / xfe_mass * 1.e6 * 86400.      !! mmol-Fe/m2/d
326      !
327      IF( nn_timing == 1 )  CALL timing_stop('trc_sed_medusa_sbc')
328      !
329   END SUBROUTINE trc_sed_medusa_sbc
330
331   !! ======================================================================
332   !! ======================================================================
333   !! ======================================================================
334
335   !! AXY & JPALM (28/02/17)
336
337   SUBROUTINE trc_sed_medusa_dust( kt )
338      !!---------------------------------------------------------------------
339      !!                     ***  ROUTINE trc_sed_medusa_dust  ***
340      !!
341      !! ** Purpose : compute current dust *before* trc_bio_medusa call
342      !!
343      !! ** Method  : does what it says on the tin
344      !!---------------------------------------------------------------------
345      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index     
346
347      !! AXY (20/11/14): alter this to report on first MEDUSA call
348      IF( kt == nittrc000 ) THEN
349         IF(lwp) WRITE(numout,*)
350         IF(lwp) WRITE(numout,*) ' trc_sed_medusa_dust: MEDUSA dust timestep'
351         IF(lwp) WRITE(numout,*) ' ~~~~~~~'
352    IF(lwp) WRITE(numout,*) ' kt =',kt
353      ENDIF
354
355      !! AXY (04/11/13): replace this with a call in trc_ini_medusa
356      !! AXY (25/02/10)
357      !! call routine for populating CCD array if this is the first time-step
358      !! IF( kt == nittrc000 ) CALL medusa_ccd( kt )
359
360      !! AXY (04/11/13): replace this with a call in trc_ini_medusa
361      !! AXY (26/01/12)
362      !! call routine for populating river arrays if this is the first time-step
363      !! IF( kt == nittrc000 ) CALL medusa_river( kt )
364
365      !! AXY (10/02/09)
366      !! IF( (jnt == 1) .and. (bdustfer) )  CALL trc_sed_medusa_sbc( kt )
367
368      !! JPALM -- 31-03-2016 -- rewrite trc_sed_medusa_sbc.
369      !! IF (kt == nittrc000 ) CALL trc_sed_medusa_sbc
370
371      !! JPALM -- 20-07-2016 -- adapt dust forcing fields reading and conversion
372      !!                     To read dust dep in kg-dust/m2/s instead of g-Fe/m2/month
373      !!                     So all forcings and coupling dust dep are in the same SI units
374      !!                     and then convert in mmol-Fe/m2/day
375
376      IF( bdustfer ) THEN
377            CALL fld_read( kt, 1, sf_dust )
378            dust(:,:) = sf_dust(1)%fnow(:,:,1)
379      ELSEIF (lk_oasis) THEN
380         dust = Dust_in_cpl
381      ELSE
382         dust(:,:) = 0.0
383      ENDIF
384      !!
385      zirondep(:,:) = 0.e0     !! Initialisation of deposition variables
386      zirondep(:,:) = dust(:,:) * Fe_dust_mratio / xfe_mass * 1.e6 * 86400.  !! mmol-Fe/m2/d
387     
388      !! JPALM -- 20-07-2016 -- Zirondep and zsidep are not used.
389      !!                     So comment out the following lines. but keep them
390      !!                     as we may want to used them later on
391      !!================================================     
392      !!
393      !! zirondep(:,:,:) = 0.e0     !! Initialisation of deposition variables
394      !! zsidep  (:,:)   = 0.e0
395      !!
396      !! Iron and Si deposition at the surface
397      !! -------------------------------------
398      !!
399      !! DO jj = 1, jpj
400      !!    DO ji = 1, jpi
401      !!       zirondep(ji,jj,1) = (dustsolub * dust(ji,jj) / (55.85 * rmtss) + 3.e-10 / ryyss) &
402      !!       & * rfact2 / fse3t(ji,jj,1)
403      !!       zsidep  (ji,jj)   = 8.8 * 0.075 * dust(ji,jj) * rfact2 / &
404      !!       & (fse3t(ji,jj,1) * 28.1 * rmtss)
405      !!    END DO
406      !! END DO
407
408   END SUBROUTINE trc_sed_medusa_dust
409
410#else
411   !!======================================================================
412   !!  Dummy module :                                   No MEDUSA bio-model
413   !!======================================================================
414CONTAINS
415   SUBROUTINE trc_sed_medusa( kt )                   ! Empty routine
416      INTEGER, INTENT( in ) ::   kt
417      WRITE(*,*) 'trc_sed_medusa: You should not have seen this print! error?', kt
418   END SUBROUTINE trc_sed_medusa
419#endif 
420
421   !!======================================================================
422END MODULE  trcsed_medusa
423
Note: See TracBrowser for help on using the repository browser.