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

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

JPALM —28-12-2015— cleaning

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