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

Last change on this file since 5740 was 5740, checked in by jpalmier, 5 years ago

JPALM —14-09-2015 — correct mistake in trcsed_medusa

File size: 13.3 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
26   USE lbclnk
27   USE prtctl_trc      ! Print control for debbuging
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   trc_sed_medusa     ! called in ???
33   PUBLIC   trc_sed_medusa_sbc     
34
35   !! * Module variables
36   INTEGER ::                   &
37     ryyss,                     &  !: number of seconds per year
38     rmtss                         !: number of seconds per month
39
40   !! AXY (10/02/09)
41   LOGICAL, PUBLIC ::                  &
42      bdustfer = .FALSE.
43   REAL(wp), PUBLIC ::                 &
44      sedfeinput = 1.e-9_wp  ,         &
45      dustsolub  = 0.014_wp
46   INTEGER ::                          &
47      numdust,                         &
48      nflx1,  nflx2,                   &
49      nflx11, nflx12
50   !!* Substitution
51#  include "domzgr_substitute.h90"
52   !!----------------------------------------------------------------------
53   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
54   !! $Id$
55   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
56   !!----------------------------------------------------------------------
57
58CONTAINS
59
60   SUBROUTINE trc_sed_medusa( kt )
61      !!---------------------------------------------------------------------
62      !!                     ***  ROUTINE trc_sed_medusa  ***
63      !!
64      !! ** Purpose :   compute the now trend due to the vertical sedimentation of
65      !!              detritus and add it to the general trend of detritus equations
66      !!
67      !! ** Method  :   this ROUTINE compute not exactly the advection but the
68      !!              transport term, i.e.  dz(wt) and dz(ws)., dz(wtr)
69      !!              using an upstream scheme
70      !!              the now vertical advection of tracers is given by:
71      !!                      dz(trn wn) = 1/bt dk+1( e1t e2t vsed (trn) )
72      !!              add this trend now to the general trend of tracer (ta,sa,tra):
73      !!                             tra = tra + dz(trn wn)
74      !!       
75      !!              IF 'key_trc_diabio' is defined, the now vertical advection
76      !!              trend of passive tracers is saved for futher diagnostics.
77      !!---------------------------------------------------------------------
78      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index     
79      !! AXY (10/02/09)
80      INTEGER  ::   jnt
81      !!
82      INTEGER  ::   ji, jj, jk
83      REAL(wp) ::   ztra
84      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwork
85
86      !! AXY (10/02/09)
87      REAL(wp), DIMENSION(jpi,jpj)     ::   zsidep    !! Si deposition
88      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zirondep  !! Fe deposition
89      REAL(wp) ::   rfact2
90
91      CHARACTER (len=25) :: charout
92      !!---------------------------------------------------------------------
93
94      !! AXY (10/02/09)
95      jnt = 1
96      rfact2 = 1.0
97
98      ! Number of seconds per year and per month
99      ryyss = nyear_len(1) * rday
100      rmtss = ryyss / raamo
101
102      !! AXY (20/11/14): alter this to report on first MEDUSA call
103      !! IF( kt == nit000 ) THEN
104      IF( kt == nittrc000 ) THEN
105         IF(lwp) WRITE(numout,*)
106         IF(lwp) WRITE(numout,*) ' trc_sed_medusa: MEDUSA sedimentation'
107         IF(lwp) WRITE(numout,*) ' ~~~~~~~'
108    IF(lwp) WRITE(numout,*) ' kt =',kt
109      ENDIF
110
111      !! AXY (04/11/13): replace this with a call in trc_ini_medusa
112      !! AXY (25/02/10)
113      !! call routine for populating CCD array if this is the first time-step
114      !! IF( kt == nittrc000 ) CALL medusa_ccd( kt )
115
116      !! AXY (04/11/13): replace this with a call in trc_ini_medusa
117      !! AXY (26/01/12)
118      !! call routine for populating river arrays if this is the first time-step
119      !! IF( kt == nittrc000 ) CALL medusa_river( kt )
120
121      !! AXY (10/02/09)
122      IF( (jnt == 1) .and. (bdustfer) )  CALL trc_sed_medusa_sbc( kt )
123      !!
124      zirondep(:,:,:) = 0.e0     !! Initialisation of deposition variables
125      zsidep  (:,:)   = 0.e0
126      !!
127      !! Iron and Si deposition at the surface
128      !! -------------------------------------
129      !!
130      DO jj = 1, jpj
131         DO ji = 1, jpi
132            zirondep(ji,jj,1) = (dustsolub * dust(ji,jj) / (55.85 * rmtss) + 3.e-10 / ryyss) &
133            & * rfact2 / fse3t(ji,jj,1)
134            zsidep  (ji,jj)   = 8.8 * 0.075 * dust(ji,jj) * rfact2 / &
135            & (fse3t(ji,jj,1) * 28.1 * rmtss)
136         END DO
137      END DO
138
139      ! sedimentation of detrital nitrogen : upstream scheme
140      ! ----------------------------------------------------
141      !
142      zwork(:,:,:) = 0.e0        ! initialisation of sinking variable
143      ! for detrital nitrogen sedimentation only - jpdet
144      zwork(:,:,1  ) = 0.e0      ! surface value set to zero
145      zwork(:,:,jpk) = 0.e0      ! bottom value  set to zero
146      !
147      ! tracer flux at w-point: we use -vsed (downward flux)  with simplification : no e1*e2
148      DO jk = 2, jpk
149         ! AXY (17/07/14): change "0.d0" to "0."
150         ! zwork(:,:,jk) = -vsed * max(trn(:,:,jk-1,jpdet),0.d0) * tmask(:,:,jk-1)
151         zwork(:,:,jk) = -vsed * max(trn(:,:,jk-1,jpdet),0.) * tmask(:,:,jk-1)
152         !
153         ! AXY (16/01/14): stop sinking in upper 10m to reduce model instability
154         !                 in shallower grid cells
155         ! if ( jk .lt. 9 ) zwork(:,:,jk) = 0.e0
156      END DO
157      !
158      ! tracer flux divergence at t-point added to the general trend
159      DO jk = 1, jpkm1
160         DO jj = 1, jpj
161            DO ji = 1,jpi
162               ztra  = - ( zwork(ji,jj,jk) - zwork(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
163               tra(ji,jj,jk,jpdet) = tra(ji,jj,jk,jpdet) + ztra
164# if defined key_trc_diabio
165               trbio(ji,jj,jk,8) = ztra
166# endif
167               IF( ln_diatrc ) &
168                  &  trc2d(ji,jj,8) = trc2d(ji,jj,8) + ztra * fse3t(ji,jj,jk) * 86400.
169            END DO
170         END DO
171      END DO
172      !
173# if defined key_trc_diabio
174      CALL lbc_lnk (trbio(:,:,1,8), 'T', 1. )                    ! Lateral boundary conditions on trcbio
175# endif
176      IF( ln_diatrc ) CALL lbc_lnk( trc2d(:,:,8), 'T', 1. )      ! Lateral boundary conditions on trc2d
177# if defined key_iomput
178      CALL iom_put( "DSED",trc2d(:,:,8) )
179# endif
180
181# if defined key_roam
182
183      ! sedimentation of detrital carbon : upstream scheme
184      ! --------------------------------------------------
185      !
186      zwork(:,:,:) = 0.e0        ! initialisation of sinking variable
187      ! for detrital carbon sedimentation only - jpdtc
188      zwork(:,:,1  ) = 0.e0      ! surface value set to zero
189      zwork(:,:,jpk) = 0.e0      ! bottom value  set to zero
190      !
191      ! tracer flux at w-point: we use -vsed (downward flux)  with simplification : no e1*e2
192      DO jk = 2, jpk
193         ! AXY (17/07/14): change "0.d0" to "0."
194         ! zwork(:,:,jk) = -vsed * max(trn(:,:,jk-1,jpdtc),0.d0) * tmask(:,:,jk-1)
195         zwork(:,:,jk) = -vsed * max(trn(:,:,jk-1,jpdtc),0.) * tmask(:,:,jk-1)
196         !
197         ! AXY (16/01/14): stop sinking in upper 10m to reduce model instability
198         !                 in shallower grid cells
199         ! if ( jk .lt. 9 ) zwork(:,:,jk) = 0.e0
200      END DO
201      !
202      ! tracer flux divergence at t-point added to the general trend
203      DO jk = 1, jpkm1
204         DO jj = 1, jpj
205            DO ji = 1,jpi
206               ztra  = - ( zwork(ji,jj,jk) - zwork(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
207               tra(ji,jj,jk,jpdtc) = tra(ji,jj,jk,jpdtc) + ztra
208!! #  if defined key_trc_diabio
209!!                trbio(ji,jj,jk,8) = ztra
210!! #  endif
211!!             IF( ln_diatrc ) &
212!!                &  trc2d(ji,jj,8) = trc2d(ji,jj,8) + ztra * fse3t(ji,jj,jk) * 86400.
213            END DO
214         END DO
215      END DO
216      !
217!! #  if defined key_trc_diabio
218!!       CALL lbc_lnk (trbio(:,:,1,8), 'T', 1. )                    ! Lateral boundary conditions on trcbio
219!! #  endif
220!!       IF( ln_diatrc ) CALL lbc_lnk( trc2d(:,:,8), 'T', 1. )      ! Lateral boundary conditions on trc2d
221!! #  if defined key_iomput
222!!       CALL iom_put( "DSED",trc2d(:,:,8) )
223!! #  endif
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   SUBROUTINE trc_sed_medusa_sbc(kt)
241
242      !!----------------------------------------------------------------------
243      !!                  ***  ROUTINE trc_sed_medusa_sbc  ***
244      !!
245      !! ** Purpose :   Read and interpolate the external sources of
246      !!                nutrients
247      !!
248      !! ** Method  :   Read the files and interpolate the appropriate variables
249      !!
250      !! ** input   :   external netcdf files
251      !!
252      !!----------------------------------------------------------------------
253      !! * arguments
254      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
255
256      !! * Local declarations
257      INTEGER ::   &
258         imois, imois2,       &  ! temporary integers
259         i15  , iman             !    "          "
260      REAL(wp) ::   &
261         zxy                     !    "         "
262
263      !!---------------------------------------------------------------------
264
265      !! Initialization
266      !! --------------
267      !!
268      i15 = nday / 16
269      iman  = INT( raamo )
270      imois = nmonth + i15 - 1
271      IF( imois == 0 ) imois = iman
272      imois2 = nmonth
273
274      !! 1. first call kt=nittrc000
275      !! -----------------------
276      !!
277      IF( kt == nittrc000 ) THEN
278         ! initializations
279         nflx1  = 0
280         nflx11 = 0
281         ! open the file
282         IF(lwp) THEN
283            WRITE(numout,*) ' '
284            WRITE(numout,*) ' **** Routine trc_sed_medusa_sbc'
285         ENDIF
286         CALL iom_open ( 'dust.orca.nc', numdust )
287    IF(lwp) WRITE(numout,*) 'trc_sed_medusa_sbc: dust.orca.nc opened'
288      ENDIF
289
290      !! Read monthly file
291      !! ----------------
292      !!
293      IF( kt == nittrc000 .OR. imois /= nflx1 ) THEN
294
295         !! Calendar computation
296         !!
297         !! nflx1 number of the first file record used in the simulation
298         !! nflx2 number of the last  file record
299         !!
300         nflx1 = imois
301         nflx2 = nflx1+1
302         nflx1 = MOD( nflx1, iman )
303         nflx2 = MOD( nflx2, iman )
304         IF( nflx1 == 0 )   nflx1 = iman
305         IF( nflx2 == 0 )   nflx2 = iman
306         IF(lwp) WRITE(numout,*) 'trc_sed_medusa_sbc: first record file used nflx1 ',nflx1
307         IF(lwp) WRITE(numout,*) 'trc_sed_medusa_sbc: last  record file used nflx2 ',nflx2
308
309         !! Read monthly fluxes data
310         !!
311         !! humidity
312         !!
313         CALL iom_get ( numdust, jpdom_data, 'dust', dustmo(:,:,1), nflx1 )
314         CALL iom_get ( numdust, jpdom_data, 'dust', dustmo(:,:,2), nflx2 )
315
316         IF(lwp .AND. nitend-nit000 <= 100 ) THEN
317            WRITE(numout,*)
318            WRITE(numout,*) ' read clio flx ok'
319            WRITE(numout,*)
320            WRITE(numout,*)
321            WRITE(numout,*) 'Clio month: ',nflx1,'  field: dust'
322            CALL prihre( dustmo(:,:,1),jpi,jpj,1,jpi,20,1,jpj,10,1e9,numout )
323         ENDIF
324
325      ENDIF
326
327      !! 3. at every time step interpolation of fluxes
328      !! ---------------------------------------------
329      !!
330      zxy = FLOAT( nday + 15 - 30 * i15 ) / 30
331      dust(:,:) = ( (1.-zxy) * dustmo(:,:,1) + zxy * dustmo(:,:,2) )
332
333      IF( kt == nitend ) THEN
334         CALL iom_close (numdust)
335         IF(lwp) WRITE(numout,*) 'trc_sed_medusa_sbc: dust.orca.nc closed'
336      ENDIF
337
338   END SUBROUTINE trc_sed_medusa_sbc
339
340#else
341   !!======================================================================
342   !!  Dummy module :                                   No MEDUSA bio-model
343   !!======================================================================
344CONTAINS
345   SUBROUTINE trc_sed_medusa( kt )                   ! Empty routine
346      INTEGER, INTENT( in ) ::   kt
347      WRITE(*,*) 'trc_sed_medusa: You should not have seen this print! error?', kt
348   END SUBROUTINE trc_sed_medusa
349#endif 
350
351   !!======================================================================
352END MODULE  trcsed_medusa
Note: See TracBrowser for help on using the repository browser.