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

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

repair dust forcing field reading freq

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