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

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

JPALM — 11-04-2016 — add dust deposition input through namelist

— relax time to IDTRA surface flux

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