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.
p4zsms.F90 in trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z – NEMO

source: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90 @ 3972

Last change on this file since 3972 was 3972, checked in by cetlod, 11 years ago

trunk: bugfix on iron budget calculation in PISCES, see ticket #1131
File size: 18.0 KB
Line 
1MODULE p4zsms
2   !!======================================================================
3   !!                         ***  MODULE p4zsms  ***
4   !! TOP :   PISCES Source Minus Sink manager
5   !!======================================================================
6   !! History :   1.0  !  2004-03 (O. Aumont) Original code
7   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
8   !!----------------------------------------------------------------------
9#if defined key_pisces
10   !!----------------------------------------------------------------------
11   !!   'key_pisces'                                       PISCES bio-model
12   !!----------------------------------------------------------------------
13   !!   p4zsms        :  Time loop of passive tracers sms
14   !!----------------------------------------------------------------------
15   USE oce_trc         !  shared variables between ocean and passive tracers
16   USE trc             !  passive tracers common variables
17   USE trcdta
18   USE sms_pisces      !  PISCES Source Minus Sink variables
19   USE p4zbio          !  Biological model
20   USE p4zche          !  Chemical model
21   USE p4zlys          !  Calcite saturation
22   USE p4zflx          !  Gas exchange
23   USE p4zsbc          !  External source of nutrients
24   USE p4zsed          !  Sedimentation
25   USE p4zint          !  time interpolation
26   USE iom             !  I/O manager
27   USE trdmod_oce      !  Ocean trends variables
28   USE trdmod_trc      !  TOP trends variables
29   USE sedmodel        !  Sediment model
30   USE prtctl_trc      !  print control for debugging
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   p4z_sms_init    ! called in p4zsms.F90
36   PUBLIC   p4z_sms    ! called in p4zsms.F90
37
38   REAL(wp) :: alkbudget, no3budget, silbudget, ferbudget
39   INTEGER ::  numco2, numnut  !: logical unit for co2 budget
40
41   !!----------------------------------------------------------------------
42   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
43   !! $Id: p4zsms.F90 3320 2012-03-05 16:37:52Z cetlod $
44   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46
47CONTAINS
48
49   SUBROUTINE p4z_sms( kt )
50      !!---------------------------------------------------------------------
51      !!                     ***  ROUTINE p4z_sms  ***
52      !!
53      !! ** Purpose :   Managment of the call to Biological sources and sinks
54      !!              routines of PISCES bio-model
55      !!
56      !! ** Method  : - at each new day ...
57      !!              - several calls of bio and sed ???
58      !!              - ...
59      !!---------------------------------------------------------------------
60      !
61      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index     
62      !!
63      INTEGER ::   jnt, jn, jl
64      CHARACTER (len=25) :: charout
65      REAL(wp), POINTER, DIMENSION(:,:,:,:)  :: ztrdpis
66      !!---------------------------------------------------------------------
67      !
68      IF( nn_timing == 1 )  CALL timing_start('p4z_sms')
69      !
70      IF( l_trdtrc )  THEN
71         CALL wrk_alloc( jpi, jpj, jpk, jp_pisces, ztrdpis ) 
72         DO jn = 1, jp_pisces
73            jl = jn + jp_pcs0 - 1
74            ztrdpis(:,:,:,jn) = trn(:,:,:,jl)
75         ENDDO
76      ENDIF
77      !
78      IF( ln_rsttr .AND. kt == nittrc000 )                         CALL p4z_rst( nittrc000, 'READ' )  !* read or initialize all required fields
79      IF( ln_pisdmp .AND. MOD( kt - nn_dttrc, nn_pisdmp ) == 0 )   CALL p4z_dmp( kt )      ! Relaxation of some tracers
80      !
81      IF( ndayflxtr /= nday_year ) THEN      ! New days
82         !
83         ndayflxtr = nday_year
84
85         IF(lwp) write(numout,*)
86         IF(lwp) write(numout,*) ' New chemical constants and various rates for biogeochemistry at new day : ', nday_year
87         IF(lwp) write(numout,*) '~~~~~~'
88
89         CALL p4z_che              ! computation of chemical constants
90         CALL p4z_int( kt )        ! computation of various rates for biogeochemistry
91         !
92      ENDIF
93
94      IF( ll_sbc ) CALL p4z_sbc( kt )   ! external sources of nutrients
95
96      DO jnt = 1, nrdttrc          ! Potential time splitting if requested
97         !
98         CALL p4z_bio (kt, jnt)    ! Biology
99         CALL p4z_sed (kt, jnt)    ! Sedimentation
100         !
101         DO jn = jp_pcs0, jp_pcs1
102            trb(:,:,:,jn) = trn(:,:,:,jn)
103         ENDDO
104         !
105      END DO
106
107      IF( l_trdtrc )  THEN
108         DO jn = 1, jp_pisces
109            jl = jn + jp_pcs0 - 1
110            ztrdpis(:,:,:,jn) = ( ztrdpis(:,:,:,jn) - trn(:,:,:,jl) ) * rfact2r
111         ENDDO
112      ENDIF
113
114      CALL p4z_lys( kt )             ! Compute CaCO3 saturation
115      CALL p4z_flx( kt )             ! Compute surface fluxes
116
117      DO jn = jp_pcs0, jp_pcs1
118        CALL lbc_lnk( trn(:,:,:,jn), 'T', 1. )
119        CALL lbc_lnk( trb(:,:,:,jn), 'T', 1. )
120        CALL lbc_lnk( tra(:,:,:,jn), 'T', 1. )
121      END DO
122      !
123      IF( lk_sed ) THEN 
124         !
125         CALL sed_model( kt )     !  Main program of Sediment model
126         !
127         DO jn = jp_pcs0, jp_pcs1
128           CALL lbc_lnk( trn(:,:,:,jn), 'T', 1. )
129         END DO
130         !
131      ENDIF
132      !
133      IF( lrst_trc )  CALL p4z_rst( kt, 'WRITE' )  !* Write PISCES informations in restart file
134      !
135      IF( l_trdtrc ) THEN
136         DO jn = 1, jp_pisces
137            jl = jn + jp_pcs0 - 1
138             ztrdpis(:,:,:,jn) = ztrdpis(:,:,:,jn) + tra(:,:,:,jl)
139             CALL trd_mod_trc( ztrdpis(:,:,:,jn), jn, jptra_trd_sms, kt )   ! save trends
140          END DO
141          CALL wrk_dealloc( jpi, jpj, jpk, jp_pisces, ztrdpis ) 
142      END IF
143      !
144      CALL p4z_chk_mass( kt ) ! Mass conservation checking
145
146      IF( nn_timing == 1 )  CALL timing_stop('p4z_sms')
147      !
148      !
149   END SUBROUTINE p4z_sms
150
151   SUBROUTINE p4z_sms_init
152      !!----------------------------------------------------------------------
153      !!                     ***  p4z_sms_init  *** 
154      !!
155      !! ** Purpose :   read PISCES namelist
156      !!
157      !! ** input   :   file 'namelist.trc.s' containing the following
158      !!             namelist: natext, natbio, natsms
159      !!                       natkriest ("key_kriest")
160      !!----------------------------------------------------------------------
161      NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, ferat3, wsbio2, niter1max, niter2max
162#if defined key_kriest
163      NAMELIST/nampiskrp/ xkr_eta, xkr_zeta, xkr_ncontent, xkr_mass_min, xkr_mass_max
164#endif
165      NAMELIST/nampisdmp/ ln_pisdmp, nn_pisdmp
166      NAMELIST/nampismass/ ln_check_mass
167      !!----------------------------------------------------------------------
168
169
170      REWIND( numnatp )
171      READ  ( numnatp, nampisbio )
172
173      IF(lwp) THEN                         ! control print
174         WRITE(numout,*) ' Namelist : nampisbio'
175         WRITE(numout,*) '    frequence pour la biologie                nrdttrc   =', nrdttrc
176         WRITE(numout,*) '    POC sinking speed                         wsbio     =', wsbio
177         WRITE(numout,*) '    half saturation constant for mortality    xkmort    =', xkmort
178         WRITE(numout,*) '    Fe/C in zooplankton                       ferat3    =', ferat3
179         WRITE(numout,*) '    Big particles sinking speed               wsbio2    =', wsbio2
180         WRITE(numout,*) '    Maximum number of iterations for POC      niter1max =', niter1max
181         WRITE(numout,*) '    Maximum number of iterations for GOC      niter2max =', niter2max
182      ENDIF
183
184#if defined key_kriest
185
186      !                               ! nampiskrp : kriest parameters
187      !                               ! -----------------------------
188      xkr_eta      = 0.62
189      xkr_zeta     = 1.62
190      xkr_ncontent = 5.7E-6
191      xkr_mass_min = 0.0002
192      xkr_mass_max = 1.
193
194      REWIND( numnatp )                     ! read natkriest
195      READ  ( numnatp, nampiskrp )
196
197      IF(lwp) THEN
198         WRITE(numout,*)
199         WRITE(numout,*) ' Namelist : nampiskrp'
200         WRITE(numout,*) '    Sinking  exponent                        xkr_eta      = ', xkr_eta
201         WRITE(numout,*) '    N content exponent                       xkr_zeta     = ', xkr_zeta
202         WRITE(numout,*) '    N content factor                         xkr_ncontent = ', xkr_ncontent
203         WRITE(numout,*) '    Minimum mass for Aggregates              xkr_mass_min = ', xkr_mass_min
204         WRITE(numout,*) '    Maximum mass for Aggregates              xkr_mass_max = ', xkr_mass_max
205         WRITE(numout,*)
206     ENDIF
207
208
209     ! Computation of some variables
210     xkr_massp = xkr_ncontent * 7.625 * xkr_mass_min**xkr_zeta
211
212#endif
213
214      ln_pisdmp = .true.
215      nn_pisdmp = 1
216
217      REWIND( numnatp )
218      READ  ( numnatp, nampisdmp )
219
220      IF(lwp) THEN                         ! control print
221         WRITE(numout,*)
222         WRITE(numout,*) ' Namelist : nampisdmp'
223         WRITE(numout,*) '    Relaxation of tracer to glodap mean value             ln_pisdmp      =', ln_pisdmp
224         WRITE(numout,*) '    Frequency of Relaxation                               nn_pisdmp      =', nn_pisdmp
225         WRITE(numout,*) ' '
226      ENDIF
227
228      ln_check_mass = .false.
229      REWIND( numnatp )       
230      READ  ( numnatp, nampismass )
231      IF(lwp) THEN                         ! control print
232         WRITE(numout,*) ' '
233         WRITE(numout,*) ' Namelist parameter for mass conservation checking'
234         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
235         WRITE(numout,*) '    Flag to check mass conservation of NO3/Si/TALK ln_check_mass = ', ln_check_mass
236      ENDIF
237
238   END SUBROUTINE p4z_sms_init
239
240   SUBROUTINE p4z_rst( kt, cdrw )
241      !!---------------------------------------------------------------------
242      !!                   ***  ROUTINE p4z_rst  ***
243      !!
244      !!  ** Purpose : Read or write variables in restart file:
245      !!
246      !!  WRITE(READ) mode:
247      !!       kt        : number of time step since the begining of the experiment at the
248      !!                   end of the current(previous) run
249      !!---------------------------------------------------------------------
250      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
251      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
252      !
253      INTEGER  ::  ji, jj, jk
254      REAL(wp) ::  zcaralk, zbicarb, zco3
255      REAL(wp) ::  ztmas, ztmas1
256      !!---------------------------------------------------------------------
257
258      IF( TRIM(cdrw) == 'READ' ) THEN
259         !
260         IF(lwp) WRITE(numout,*)
261         IF(lwp) WRITE(numout,*) ' p4z_rst : Read specific variables from pisces model '
262         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
263         !
264         IF( iom_varid( numrtr, 'PH', ldstop = .FALSE. ) > 0 ) THEN
265            CALL iom_get( numrtr, jpdom_autoglo, 'PH' , hi(:,:,:)  )
266         ELSE
267!            hi(:,:,:) = 1.e-9
268            ! Set PH from  total alkalinity, borat (???), akb3 (???) and ak23 (???)
269            ! --------------------------------------------------------
270            DO jk = 1, jpk
271               DO jj = 1, jpj
272                  DO ji = 1, jpi
273                     ztmas   = tmask(ji,jj,jk)
274                     ztmas1  = 1. - tmask(ji,jj,jk)
275                     zcaralk = trn(ji,jj,jk,jptal) - borat(ji,jj,jk) / (  1. + 1.E-8 / ( rtrn + akb3(ji,jj,jk) )  )
276                     zco3    = ( zcaralk - trn(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1
277                     zbicarb = ( 2. * trn(ji,jj,jk,jpdic) - zcaralk )
278                     hi(ji,jj,jk) = ( ak23(ji,jj,jk) * zbicarb / zco3 ) * ztmas + 1.e-9 * ztmas1
279                  END DO
280               END DO
281            END DO
282         ENDIF
283         CALL iom_get( numrtr, jpdom_autoglo, 'Silicalim', xksi(:,:) )
284         IF( iom_varid( numrtr, 'Silicamax', ldstop = .FALSE. ) > 0 ) THEN
285            CALL iom_get( numrtr, jpdom_autoglo, 'Silicamax' , xksimax(:,:)  )
286         ELSE
287            xksimax(:,:) = xksi(:,:)
288         ENDIF
289         !
290      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
291         IF( kt == nitrst ) THEN
292            IF(lwp) WRITE(numout,*)
293            IF(lwp) WRITE(numout,*) 'p4z_rst : write pisces restart file  kt =', kt
294            IF(lwp) WRITE(numout,*) '~~~~~~~'
295         ENDIF
296         CALL iom_rstput( kt, nitrst, numrtw, 'PH', hi(:,:,:) )
297         CALL iom_rstput( kt, nitrst, numrtw, 'Silicalim', xksi(:,:) )
298         CALL iom_rstput( kt, nitrst, numrtw, 'Silicamax', xksimax(:,:) )
299      ENDIF
300      !
301   END SUBROUTINE p4z_rst
302
303   SUBROUTINE p4z_dmp( kt )
304      !!----------------------------------------------------------------------
305      !!                    ***  p4z_dmp  ***
306      !!
307      !! ** purpose  : Relaxation of some tracers
308      !!----------------------------------------------------------------------
309      !
310      INTEGER, INTENT( in )  ::     kt ! time step
311      !
312      REAL(wp) ::  alkmean = 2426.     ! mean value of alkalinity ( Glodap ; for Goyet 2391. )
313      REAL(wp) ::  po4mean = 2.165     ! mean value of phosphates
314      REAL(wp) ::  no3mean = 30.90     ! mean value of nitrate
315      REAL(wp) ::  silmean = 91.51     ! mean value of silicate
316      !
317      REAL(wp) :: zarea, zalksum, zpo4sum, zno3sum, zsilsum
318      !!---------------------------------------------------------------------
319
320
321      IF(lwp)  WRITE(numout,*)
322      IF(lwp)  WRITE(numout,*) ' p4z_dmp : Restoring of nutrients at time-step kt = ', kt
323      IF(lwp)  WRITE(numout,*)
324
325      IF( cp_cfg == "orca" .AND. .NOT. lk_c1d ) THEN      ! ORCA configuration (not 1D) !
326         !                                                    ! --------------------------- !
327         ! set total alkalinity, phosphate, nitrate & silicate
328         zarea          = 1._wp / glob_sum( cvol(:,:,:) ) * 1e6             
329
330         zalksum = glob_sum( trn(:,:,:,jptal) * cvol(:,:,:)  ) * zarea
331         zpo4sum = glob_sum( trn(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea * po4r
332         zno3sum = glob_sum( trn(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea * rno3
333         zsilsum = glob_sum( trn(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea
334 
335         IF(lwp) WRITE(numout,*) '       TALK mean : ', zalksum
336         trn(:,:,:,jptal) = trn(:,:,:,jptal) * alkmean / zalksum
337
338         IF(lwp) WRITE(numout,*) '       PO4  mean : ', zpo4sum
339         trn(:,:,:,jppo4) = trn(:,:,:,jppo4) * po4mean / zpo4sum
340
341         IF(lwp) WRITE(numout,*) '       NO3  mean : ', zno3sum
342         trn(:,:,:,jpno3) = trn(:,:,:,jpno3) * no3mean / zno3sum
343
344         IF(lwp) WRITE(numout,*) '       SiO3 mean : ', zsilsum
345         trn(:,:,:,jpsil) = MIN( 400.e-6,trn(:,:,:,jpsil) * silmean / zsilsum )
346         !
347      ENDIF
348
349   END SUBROUTINE p4z_dmp
350
351
352   SUBROUTINE p4z_chk_mass( kt )
353      !!----------------------------------------------------------------------
354      !!                  ***  ROUTINE p4z_chk_mass  ***
355      !!
356      !! ** Purpose :  Mass conservation check
357      !!
358      !!---------------------------------------------------------------------
359      !
360      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index     
361      !!
362      !!---------------------------------------------------------------------
363
364      IF( kt == nittrc000 ) THEN
365         IF( ln_check_mass .AND. lwp) THEN      !   Open budget file of NO3, ALK, Si, Fer
366            CALL ctl_opn( numco2, 'carbon.budget'  , 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
367            CALL ctl_opn( numnut, 'nutrient.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
368         ENDIF
369      ENDIF
370
371      IF( ln_check_mass .AND. kt == nitend ) THEN      !   Compute the budget of NO3, ALK, Si, Fer
372         no3budget = glob_sum( (   trn(:,:,:,jpno3) + trn(:,:,:,jpnh4)  &
373            &                    + trn(:,:,:,jpphy) + trn(:,:,:,jpdia)  &
374            &                    + trn(:,:,:,jpzoo) + trn(:,:,:,jpmes)  &
375            &                    + trn(:,:,:,jppoc)                     &
376#if ! defined key_kriest
377            &                    + trn(:,:,:,jpgoc)                     &
378#endif
379            &                    + trn(:,:,:,jpdoc)                     ) * cvol(:,:,:)  ) 
380         !
381         silbudget = glob_sum( (   trn(:,:,:,jpsil) + trn(:,:,:,jpgsi)  &
382            &                    + trn(:,:,:,jpdsi)                     ) * cvol(:,:,:)  )
383         !
384         alkbudget = glob_sum( (   trn(:,:,:,jpno3) * rno3              &
385            &                    + trn(:,:,:,jptal)                     &
386            &                    + trn(:,:,:,jpcal) * 2.                ) * cvol(:,:,:)  )
387         !
388         ferbudget = glob_sum( (   trn(:,:,:,jpfer) + trn(:,:,:,jpnfe)  &
389            &                    + trn(:,:,:,jpdfe)                     &
390#if ! defined key_kriest
391            &                    + trn(:,:,:,jpbfe)                     &
392#endif
393            &                    + trn(:,:,:,jpsfe)                     &
394            &                    + trn(:,:,:,jpzoo) * ferat3            &
395            &                    + trn(:,:,:,jpmes) * ferat3            ) * cvol(:,:,:)  )
396
397         !
398         t_atm_co2_flx  = t_atm_co2_flx / glob_sum( e1e2t(:,:) )
399         t_oce_co2_flx  = t_oce_co2_flx         * 12. / 1.e15 * (-1 )
400         tpp            = tpp           * 1000. * 12. / 1.E15
401         t_oce_co2_exp  = t_oce_co2_exp * 1000. * 12. / 1.E15
402         !
403         no3budget = no3budget / areatot
404         silbudget = silbudget / areatot
405         alkbudget = alkbudget / areatot
406         ferbudget = ferbudget / areatot
407         !
408         IF(lwp) THEN
409            WRITE(numco2,9000) ndastp, t_atm_co2_flx, t_oce_co2_flx, tpp, t_oce_co2_exp
410            WRITE(numnut,9500) ndastp, alkbudget, no3budget, silbudget, ferbudget
411         ENDIF
412         !
413      ENDIF
414       !
415 9000  FORMAT(i8,f10.5,e18.10,f10.5,f10.5)
416 9500  FORMAT(i8,4e18.10)     
417       !
418   END SUBROUTINE p4z_chk_mass
419
420#else
421   !!======================================================================
422   !!  Dummy module :                                   No PISCES bio-model
423   !!======================================================================
424CONTAINS
425   SUBROUTINE p4z_sms( kt )                   ! Empty routine
426      INTEGER, INTENT( in ) ::   kt
427      WRITE(*,*) 'p4z_sms: You should not have seen this print! error?', kt
428   END SUBROUTINE p4z_sms
429#endif 
430
431   !!======================================================================
432END MODULE p4zsms 
Note: See TracBrowser for help on using the repository browser.