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 @ 5931

Last change on this file since 5931 was 5931, checked in by jpalmier, 8 years ago

JPALM --26-11-2015 -- Update MEDUSA diagnostics with iom_use

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                  .
188            END DO
189         END DO
190      END DO
191      !
192# if defined key_trc_diabio
193      CALL lbc_lnk (trbio(:,:,1,8), 'T', 1. )                    ! Lateral boundary conditions on trcbio
194# endif
195      IF( ln_diatrc ) CALL lbc_lnk( trc2d(:,:,8), 'T', 1. )      ! Lateral boundary conditions on trc2d
196      !!
197      IF (lk_iomput .AND. .NOT. ln_diatrc) THEN
198           IF( med_diag%DSED%dgsave ) THEN
199                CALL iom_put( "DSED"  ,  zw2d)
200                CALL wrk_dealloc( jpi, jpj,    zw2d  )
201            ENDIF
202      ELSE IF (lk_iomput .AND. ln_diatrc)  THEN   
203          CALL iom_put( "DSED",trc2d(:,:,8) )
204      ENDIF
205      !!
206# if defined key_roam
207
208      ! sedimentation of detrital carbon : upstream scheme
209      ! --------------------------------------------------
210      !
211      zwork(:,:,:) = 0.e0        ! initialisation of sinking variable
212      ! for detrital carbon sedimentation only - jpdtc
213      zwork(:,:,1  ) = 0.e0      ! surface value set to zero
214      zwork(:,:,jpk) = 0.e0      ! bottom value  set to zero
215      !
216      ! tracer flux at w-point: we use -vsed (downward flux)  with simplification : no e1*e2
217      DO jk = 2, jpk
218         ! AXY (17/07/14): change "0.d0" to "0."
219         ! zwork(:,:,jk) = -vsed * max(trn(:,:,jk-1,jpdtc),0.d0) * tmask(:,:,jk-1)
220         zwork(:,:,jk) = -vsed * max(trn(:,:,jk-1,jpdtc),0.) * tmask(:,:,jk-1)
221         !
222         ! AXY (16/01/14): stop sinking in upper 10m to reduce model instability
223         !                 in shallower grid cells
224         ! if ( jk .lt. 9 ) zwork(:,:,jk) = 0.e0
225      END DO
226      !
227      ! tracer flux divergence at t-point added to the general trend
228      DO jk = 1, jpkm1
229         DO jj = 1, jpj
230            DO ji = 1,jpi
231               ztra  = - ( zwork(ji,jj,jk) - zwork(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
232               tra(ji,jj,jk,jpdtc) = tra(ji,jj,jk,jpdtc) + ztra
233!! #  if defined key_trc_diabio
234!!                trbio(ji,jj,jk,8) = ztra
235!! #  endif
236!!             IF( ln_diatrc ) &
237!!                &  trc2d(ji,jj,8) = trc2d(ji,jj,8) + ztra * fse3t(ji,jj,jk) * 86400.
238            END DO
239         END DO
240      END DO
241      !
242!! #  if defined key_trc_diabio
243!!       CALL lbc_lnk (trbio(:,:,1,8), 'T', 1. )                    ! Lateral boundary conditions on trcbio
244!! #  endif
245!!       IF( ln_diatrc ) CALL lbc_lnk( trc2d(:,:,8), 'T', 1. )      ! Lateral boundary conditions on trc2d
246!! #  if defined key_iomput
247!!       CALL iom_put( "DSED",trc2d(:,:,8) )
248!! #  endif
249
250# endif
251
252      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
253         WRITE(charout, FMT="('sed')")
254         CALL prt_ctl_trc_info(charout)
255         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
256      ENDIF
257
258   END SUBROUTINE trc_sed_medusa
259
260   !! ======================================================================
261   !! ======================================================================
262   !! ======================================================================
263
264   !! AXY (10/02/09)
265   SUBROUTINE trc_sed_medusa_sbc(kt)
266
267      !!----------------------------------------------------------------------
268      !!                  ***  ROUTINE trc_sed_medusa_sbc  ***
269      !!
270      !! ** Purpose :   Read and interpolate the external sources of
271      !!                nutrients
272      !!
273      !! ** Method  :   Read the files and interpolate the appropriate variables
274      !!
275      !! ** input   :   external netcdf files
276      !!
277      !!----------------------------------------------------------------------
278      !! * arguments
279      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
280
281      !! * Local declarations
282      INTEGER ::   &
283         imois, imois2,       &  ! temporary integers
284         i15  , iman             !    "          "
285      REAL(wp) ::   &
286         zxy                     !    "         "
287
288      !!---------------------------------------------------------------------
289      IF (bdustfer) THEN
290         !! Initialization
291         !! --------------
292         !!
293         i15 = nday / 16
294         iman  = INT( raamo )
295         imois = nmonth + i15 - 1
296         IF( imois == 0 ) imois = iman
297         imois2 = nmonth
298
299         !! 1. first call kt=nittrc000
300         !! -----------------------
301         !!
302         IF (kt == nittrc000) THEN
303            ! initializations
304            nflx1  = 0
305            nflx11 = 0
306            ! open the file
307            IF(lwp) THEN
308               WRITE(numout,*) ' '
309               WRITE(numout,*) ' **** Routine trc_sed_medusa_sbc'
310            ENDIF
311            CALL iom_open ( 'dust.orca.nc', numdust )
312            IF(lwp) WRITE(numout,*) 'trc_sed_medusa_sbc: dust.orca.nc opened'
313         ENDIF
314   
315         !! Read monthly file
316         !! ----------------
317         !!
318         IF( kt == nittrc000 .OR. imois /= nflx1 ) THEN
319
320            !! Calendar computation
321            !!
322            !! nflx1 number of the first file record used in the simulation
323            !! nflx2 number of the last  file record
324            !!
325            nflx1 = imois
326            nflx2 = nflx1+1
327            nflx1 = MOD( nflx1, iman )
328            nflx2 = MOD( nflx2, iman )
329            IF( nflx1 == 0 )   nflx1 = iman
330            IF( nflx2 == 0 )   nflx2 = iman
331            IF(lwp) WRITE(numout,*) 'trc_sed_medusa_sbc: first record file used nflx1 ',nflx1
332            IF(lwp) WRITE(numout,*) 'trc_sed_medusa_sbc: last  record file used nflx2 ',nflx2
333
334            !! Read monthly fluxes data
335            !!
336            !! humidity
337            !!
338            CALL iom_get ( numdust, jpdom_data, 'dust', dustmo(:,:,1), nflx1 )
339            CALL iom_get ( numdust, jpdom_data, 'dust', dustmo(:,:,2), nflx2 )
340
341            IF(lwp .AND. nitend-nit000 <= 100 ) THEN
342               WRITE(numout,*)
343               WRITE(numout,*) ' read clio flx ok'
344               WRITE(numout,*)
345               WRITE(numout,*)
346               WRITE(numout,*) 'Clio month: ',nflx1,'  field: dust'
347               CALL prihre( dustmo(:,:,1),jpi,jpj,1,jpi,20,1,jpj,10,1e9,numout )
348            ENDIF
349
350         ENDIF
351
352         !! 3. at every time step interpolation of fluxes
353         !! ---------------------------------------------
354         !!
355         zxy = FLOAT( nday + 15 - 30 * i15 ) / 30
356         dust(:,:) = ( (1.-zxy) * dustmo(:,:,1) + zxy * dustmo(:,:,2) )
357
358         IF( kt == nitend ) THEN
359            CALL iom_close (numdust)
360            IF(lwp) WRITE(numout,*) 'trc_sed_medusa_sbc: dust.orca.nc closed'
361         ENDIF
362      ELSE
363         dust(:,:) = 0.0
364      ENDIF
365   END SUBROUTINE trc_sed_medusa_sbc
366
367#else
368   !!======================================================================
369   !!  Dummy module :                                   No MEDUSA bio-model
370   !!======================================================================
371CONTAINS
372   SUBROUTINE trc_sed_medusa( kt )                   ! Empty routine
373      INTEGER, INTENT( in ) ::   kt
374      WRITE(*,*) 'trc_sed_medusa: You should not have seen this print! error?', kt
375   END SUBROUTINE trc_sed_medusa
376#endif 
377
378   !!======================================================================
379END MODULE  trcsed_medusa
Note: See TracBrowser for help on using the repository browser.