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

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

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

File size: 13.6 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      IF (bdustfer) THEN
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      ELSE
338         dust(:,:) = 0.0
339      ENDIF
340   END SUBROUTINE trc_sed_medusa_sbc
341
342#else
343   !!======================================================================
344   !!  Dummy module :                                   No MEDUSA bio-model
345   !!======================================================================
346CONTAINS
347   SUBROUTINE trc_sed_medusa( kt )                   ! Empty routine
348      INTEGER, INTENT( in ) ::   kt
349      WRITE(*,*) 'trc_sed_medusa: You should not have seen this print! error?', kt
350   END SUBROUTINE trc_sed_medusa
351#endif 
352
353   !!======================================================================
354END MODULE  trcsed_medusa
Note: See TracBrowser for help on using the repository browser.