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

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

Last change on this file since 7894 was 7894, checked in by jpalmier, 7 years ago

JPALM -- 11-04-2017 -- MEDUSA spring tidy-up refreshning session
File size: 14.9 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      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index     
94      !! AXY (10/02/09)
95      INTEGER  ::   jnt
96      !!
97      INTEGER  ::   ji, jj, jk
98      REAL(wp) ::   ztra
99      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwork
100
101      !! AXY (10/02/09)
102      REAL(wp) ::   rfact2
103
104      CHARACTER (len=25) :: charout
105     
106      !! JPALM - 26-11-2015 -add iom_use for diagnostic
107       REAL(wp), POINTER, DIMENSION(:,:  ) :: zw2d
108      !!---------------------------------------------------------------------
109      !!
110      IF( lk_iomput) THEN 
111           IF( med_diag%DSED%dgsave ) THEN
112               CALL wrk_alloc( jpi, jpj,      zw2d )
113                zw2d(:,:)      = 0.0      !!
114           ENDIF
115      ENDIF
116     
117      !! AXY (10/02/09)
118      jnt = 1
119      rfact2 = 1.0
120
121      ! Number of seconds per year and per month
122      ryyss = nyear_len(1) * rday
123      rmtss = ryyss / raamo
124
125      !! AXY (20/11/14): alter this to report on first MEDUSA call
126      IF( kt == nittrc000 ) THEN
127         IF(lwp) WRITE(numout,*)
128         IF(lwp) WRITE(numout,*) ' trc_sed_medusa: MEDUSA sedimentation'
129         IF(lwp) WRITE(numout,*) ' ~~~~~~~'
130    IF(lwp) WRITE(numout,*) ' kt =',kt
131      ENDIF
132
133      ! sedimentation of detrital nitrogen : upstream scheme
134      ! ----------------------------------------------------
135      !
136      zwork(:,:,:) = 0.e0        ! initialisation of sinking variable
137      ! for detrital nitrogen sedimentation only - jpdet
138      zwork(:,:,1  ) = 0.e0      ! surface value set to zero
139      zwork(:,:,jpk) = 0.e0      ! bottom value  set to zero
140      !
141      ! tracer flux at w-point: we use -vsed (downward flux)  with simplification : no e1*e2
142      DO jk = 2, jpk
143         ! AXY (17/07/14): change "0.d0" to "0."
144         ! zwork(:,:,jk) = -vsed * max(trn(:,:,jk-1,jpdet),0.d0) * tmask(:,:,jk-1)
145         zwork(:,:,jk) = -vsed * max(trn(:,:,jk-1,jpdet),0.) * tmask(:,:,jk-1)
146      END DO
147      !
148      ! tracer flux divergence at t-point added to the general trend
149      DO jk = 1, jpkm1
150         DO jj = 1, jpj
151            DO ji = 1,jpi
152               ztra  = - ( zwork(ji,jj,jk) - zwork(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
153               tra(ji,jj,jk,jpdet) = tra(ji,jj,jk,jpdet) + ztra
154               IF (lk_iomput .AND. .NOT. ln_diatrc) THEN
155                     IF( med_diag%DSED%dgsave ) THEN
156                         zw2d(ji,jj) = zw2d(ji,jj) + ztra * fse3t(ji,jj,jk) * 86400.
157                      ENDIF   
158               ENDIF   
159               
160            END DO
161         END DO
162      END DO
163      !
164      IF (lk_iomput .AND. .NOT. ln_diatrc) THEN
165           IF( med_diag%DSED%dgsave ) THEN
166                CALL iom_put( "DSED"  ,  zw2d)
167                CALL wrk_dealloc( jpi, jpj,    zw2d  )
168            ENDIF
169      ENDIF
170      !!
171# if defined key_roam
172
173      ! sedimentation of detrital carbon : upstream scheme
174      ! --------------------------------------------------
175      !
176      zwork(:,:,:) = 0.e0        ! initialisation of sinking variable
177      ! for detrital carbon sedimentation only - jpdtc
178      zwork(:,:,1  ) = 0.e0      ! surface value set to zero
179      zwork(:,:,jpk) = 0.e0      ! bottom value  set to zero
180      !
181      ! tracer flux at w-point: we use -vsed (downward flux)  with simplification : no e1*e2
182      DO jk = 2, jpk
183         ! AXY (17/07/14): change "0.d0" to "0."
184         ! zwork(:,:,jk) = -vsed * max(trn(:,:,jk-1,jpdtc),0.d0) * tmask(:,:,jk-1)
185         zwork(:,:,jk) = -vsed * max(trn(:,:,jk-1,jpdtc),0.) * tmask(:,:,jk-1)
186      END DO
187      !
188      ! tracer flux divergence at t-point added to the general trend
189      DO jk = 1, jpkm1
190         DO jj = 1, jpj
191            DO ji = 1,jpi
192               ztra  = - ( zwork(ji,jj,jk) - zwork(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
193               tra(ji,jj,jk,jpdtc) = tra(ji,jj,jk,jpdtc) + ztra
194            END DO
195         END DO
196      END DO
197      !
198# endif
199
200      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
201         WRITE(charout, FMT="('sed')")
202         CALL prt_ctl_trc_info(charout)
203         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
204      ENDIF
205
206   END SUBROUTINE trc_sed_medusa
207
208   !! ======================================================================
209   !! ======================================================================
210   !! ======================================================================
211
212   !! AXY (10/02/09)
213   !! JPALM -- 31-03-2016 -- Completely change trc_sed_medusa_sbc.
214   !!                     -- We now need to read dust file through a namelist.
215   !!                     To be able to use time varying dust depositions from
216   !!                     -- copy and adapt the PISCES p4z_sbc_ini subroutine
217   !!                     -- Only use the dust related part.     
218   SUBROUTINE trc_sed_medusa_sbc(kt)
219
220      !!----------------------------------------------------------------------
221      !!                  ***  ROUTINE trc_sed_medusa_sbc  ***
222      !!
223      !! ** Purpose :   Read and dust namelist and files.
224      !!                The interpolation is done in trc_sed through
225      !!                "CALL fld_read( kt, 1, sf_dust )"
226      !!
227      !! ** Method  :   Read the sbc namelist, and the adapted dust file, if required
228      !!                called at the first timestep (nittrc000)
229      !!
230      !! ** input   :   -- namelist sbc ref and cfg
231      !!                -- external netcdf files
232      !!
233      !!----------------------------------------------------------------------
234      !! * arguments
235      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
236      INTEGER  :: ji, jj, jk, jm, ifpr
237      INTEGER  :: ii0, ii1, ij0, ij1
238      INTEGER  :: numdust
239      INTEGER  :: ierr 
240      INTEGER  :: jfld                ! dummy loop arguments
241      INTEGER  :: ios                 ! Local integer output status for namelist read
242      INTEGER  :: isrow               ! index for ORCA1 starting row
243      REAL(wp) :: ztimes_dust
244      REAL(wp), DIMENSION(nbtimes) :: zsteps                 ! times records
245      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zdust
246      !
247      CHARACTER(len=100)         ::  cn_dir       ! Root directory for location of ssr files
248      TYPE(FLD_N), DIMENSION(1)  ::   slf_d       ! array of namelist informations on the fields to read
249      TYPE(FLD_N)                ::   sn_dust     ! informations about the fields to be read
250      !
251      NAMELIST/nammedsbc/cn_dir, sn_dust, bdustfer 
252
253      !!---------------------------------------------------------------------
254      !
255      IF( nn_timing == 1 )  CALL timing_start('trc_sed_medusa_sbc')
256      !
257      !                            !* set file information
258      REWIND( numnatp_ref )        ! Namelist nammedsbc in reference namelist : MEDUSA external sources of Dust
259      READ  ( numnatp_ref, nammedsbc, IOSTAT = ios, ERR = 901)
260901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammedsbc in reference namelist', lwp )
261
262      REWIND( numnatp_cfg )        ! Namelist nammedsbc in configuration namelist : MEDUSA external sources of Dust
263      READ  ( numnatp_cfg, nammedsbc, IOSTAT = ios, ERR = 902 )
264902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammedsbc in configuration namelist', lwp )
265      IF(lwm) WRITE ( numonp, nammedsbc )
266
267      IF(lwp) THEN
268          WRITE(numout,*) ' '
269          WRITE(numout,*) ' namelist : nammedsbc '
270          WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ '
271          WRITE(numout,*) '    dust input from the atmosphere bdustfer     = ', bdustfer
272      END IF
273
274      ! dust input from the atmosphere
275      ! ------------------------------
276      IF( bdustfer ) THEN
277         !
278         IF(lwp) WRITE(numout,*) '    initialize dust input from atmosphere '
279         IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
280         !
281         !! already allocated in sms_medusa
282         !!ALLOCATE( dust(jpi,jpj) )    ! allocation
283         !
284         slf_d(1) = sn_dust                          ! store namelist information in an array
285         !
286         ALLOCATE( sf_dust(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst
287         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'trc_sed_medusa_sbc: unable to allocate sf_dust structure' )
288         ALLOCATE( sf_dust(1)%fnow(jpi,jpj,1))
289         IF( slf_d(1)%ln_tint )     ALLOCATE( sf_dust(1)%fdta(jpi,jpj,1,2) )
290         !
291         CALL fld_fill( sf_dust, slf_d, cn_dir, 'trc_sed_medusa_sbc', 'Atmospheric dust deposition', 'nammedsed' )
292         !
293         CALL fld_read( kt, 1, sf_dust )
294         dust(:,:) = sf_dust(1)%fnow(:,:,1)
295         !
296      ELSEIF (lk_oasis) THEN
297         dust = Dust_in_cpl
298      ELSE
299         dust(:,:) = 0.0
300      END IF
301      !
302      zirondep(:,:) = 0.e0     !! Initialisation of deposition variables
303      zirondep(:,:) = dust(:,:) * Fe_dust_mratio / xfe_mass * 1.e6 * 86400.      !! mmol-Fe/m2/d
304      !
305      IF( nn_timing == 1 )  CALL timing_stop('trc_sed_medusa_sbc')
306      !
307   END SUBROUTINE trc_sed_medusa_sbc
308
309   !! ======================================================================
310   !! ======================================================================
311   !! ======================================================================
312
313   !! AXY & JPALM (28/02/17)
314
315   SUBROUTINE trc_sed_medusa_dust( kt )
316      !!---------------------------------------------------------------------
317      !!                     ***  ROUTINE trc_sed_medusa_dust  ***
318      !!
319      !! ** Purpose : compute current dust *before* trc_bio_medusa call
320      !!
321      !! ** Method  : does what it says on the tin
322      !!---------------------------------------------------------------------
323      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index     
324
325      !! AXY (20/11/14): alter this to report on first MEDUSA call
326      IF( kt == nittrc000 ) THEN
327         IF(lwp) WRITE(numout,*)
328         IF(lwp) WRITE(numout,*) ' trc_sed_medusa_dust: MEDUSA dust timestep'
329         IF(lwp) WRITE(numout,*) ' ~~~~~~~'
330    IF(lwp) WRITE(numout,*) ' kt =',kt
331      ENDIF
332
333      IF( bdustfer ) THEN
334            CALL fld_read( kt, 1, sf_dust )
335            dust(:,:) = sf_dust(1)%fnow(:,:,1)
336      ELSEIF (lk_oasis) THEN
337         dust = Dust_in_cpl
338      ELSE
339         dust(:,:) = 0.0
340      ENDIF
341      !!
342      zirondep(:,:) = 0.e0     !! Initialisation of deposition variables
343      zirondep(:,:) = dust(:,:) * Fe_dust_mratio / xfe_mass * 1.e6 * 86400.  !! mmol-Fe/m2/d
344     
345   END SUBROUTINE trc_sed_medusa_dust
346
347#else
348   !!======================================================================
349   !!  Dummy module :                                   No MEDUSA bio-model
350   !!======================================================================
351CONTAINS
352   SUBROUTINE trc_sed_medusa( kt )                   ! Empty routine
353      INTEGER, INTENT( in ) ::   kt
354      WRITE(*,*) 'trc_sed_medusa: You should not have seen this print! error?', kt
355   END SUBROUTINE trc_sed_medusa
356#endif 
357
358   !!======================================================================
359END MODULE  trcsed_medusa
360
Note: See TracBrowser for help on using the repository browser.