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 branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z – NEMO

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90 @ 4148

Last change on this file since 4148 was 4148, checked in by cetlod, 10 years ago

merge in trunk changes between r3853 and r3940 and commit the changes, see ticket #1169
File size: 19.8 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 ( kt == nittrc000 ) CALL FLUSH    ( numonp )     ! flush output namelist PISCES
147      IF( nn_timing == 1 )  CALL timing_stop('p4z_sms')
148      !
149      !
150   END SUBROUTINE p4z_sms
151
152   SUBROUTINE p4z_sms_init
153      !!----------------------------------------------------------------------
154      !!                     ***  p4z_sms_init  *** 
155      !!
156      !! ** Purpose :   read PISCES namelist
157      !!
158      !! ** input   :   file 'namelist.trc.s' containing the following
159      !!             namelist: natext, natbio, natsms
160      !!                       natkriest ("key_kriest")
161      !!----------------------------------------------------------------------
162      NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, ferat3, wsbio2, niter1max, niter2max
163#if defined key_kriest
164      NAMELIST/nampiskrp/ xkr_eta, xkr_zeta, xkr_ncontent, xkr_mass_min, xkr_mass_max
165#endif
166      NAMELIST/nampisdmp/ ln_pisdmp, nn_pisdmp
167      NAMELIST/nampismass/ ln_check_mass
168      INTEGER :: ios                 ! Local integer output status for namelist read
169      !!----------------------------------------------------------------------
170
171      REWIND( numnatp_ref )              ! Namelist nampisbio in reference namelist : Pisces variables
172      READ  ( numnatp_ref, nampisbio, IOSTAT = ios, ERR = 901)
173901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisbio in reference namelist', lwp )
174
175      REWIND( numnatp_cfg )              ! Namelist nampisbio in configuration namelist : Pisces variables
176      READ  ( numnatp_cfg, nampisbio, IOSTAT = ios, ERR = 902 )
177902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisbio in configuration namelist', lwp )
178      WRITE ( numonp, nampisbio )
179
180      IF(lwp) THEN                         ! control print
181         WRITE(numout,*) ' Namelist : nampisbio'
182         WRITE(numout,*) '    frequence pour la biologie                nrdttrc   =', nrdttrc
183         WRITE(numout,*) '    POC sinking speed                         wsbio     =', wsbio
184         WRITE(numout,*) '    half saturation constant for mortality    xkmort    =', xkmort
185         WRITE(numout,*) '    Fe/C in zooplankton                       ferat3    =', ferat3
186         WRITE(numout,*) '    Big particles sinking speed               wsbio2    =', wsbio2
187         WRITE(numout,*) '    Maximum number of iterations for POC      niter1max =', niter1max
188         WRITE(numout,*) '    Maximum number of iterations for GOC      niter2max =', niter2max
189      ENDIF
190
191#if defined key_kriest
192
193      !                               ! nampiskrp : kriest parameters
194      !                               ! -----------------------------
195      REWIND( numnatp_ref )              ! Namelist nampiskrp in reference namelist : Pisces Kriest
196      READ  ( numnatp_ref, nampiskrp, IOSTAT = ios, ERR = 903)
197903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampiskrp in reference namelist', lwp )
198
199      REWIND( numnatp_cfg )              ! Namelist nampiskrp in configuration namelist : Pisces Kriest
200      READ  ( numnatp_cfg, nampiskrp, IOSTAT = ios, ERR = 904 )
201904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampiskrp in configuration namelist', lwp )
202      WRITE ( numonp, nampiskrp )
203
204      IF(lwp) THEN
205         WRITE(numout,*)
206         WRITE(numout,*) ' Namelist : nampiskrp'
207         WRITE(numout,*) '    Sinking  exponent                        xkr_eta      = ', xkr_eta
208         WRITE(numout,*) '    N content exponent                       xkr_zeta     = ', xkr_zeta
209         WRITE(numout,*) '    N content factor                         xkr_ncontent = ', xkr_ncontent
210         WRITE(numout,*) '    Minimum mass for Aggregates              xkr_mass_min = ', xkr_mass_min
211         WRITE(numout,*) '    Maximum mass for Aggregates              xkr_mass_max = ', xkr_mass_max
212         WRITE(numout,*)
213     ENDIF
214
215
216     ! Computation of some variables
217     xkr_massp = xkr_ncontent * 7.625 * xkr_mass_min**xkr_zeta
218
219#endif
220
221      REWIND( numnatp_ref )              ! Namelist nampisdmp in reference namelist : Pisces damping
222      READ  ( numnatp_ref, nampisdmp, IOSTAT = ios, ERR = 905)
223905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisdmp in reference namelist', lwp )
224
225      REWIND( numnatp_cfg )              ! Namelist nampisdmp in configuration namelist : Pisces damping
226      READ  ( numnatp_cfg, nampisdmp, IOSTAT = ios, ERR = 906 )
227906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisdmp in configuration namelist', lwp )
228      WRITE ( numonp, nampisdmp )
229
230      IF(lwp) THEN                         ! control print
231         WRITE(numout,*)
232         WRITE(numout,*) ' Namelist : nampisdmp'
233         WRITE(numout,*) '    Relaxation of tracer to glodap mean value             ln_pisdmp      =', ln_pisdmp
234         WRITE(numout,*) '    Frequency of Relaxation                               nn_pisdmp      =', nn_pisdmp
235         WRITE(numout,*) ' '
236      ENDIF
237
238      REWIND( numnatp_ref )              ! Namelist nampismass in reference namelist : Pisces mass conservation check
239      READ  ( numnatp_ref, nampismass, IOSTAT = ios, ERR = 907)
240907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampismass in reference namelist', lwp )
241
242      REWIND( numnatp_cfg )              ! Namelist nampismass in configuration namelist : Pisces mass conservation check
243      READ  ( numnatp_cfg, nampismass, IOSTAT = ios, ERR = 908 )
244908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampismass in configuration namelist', lwp )
245      WRITE ( numonp, nampismass )
246
247      IF(lwp) THEN                         ! control print
248         WRITE(numout,*) ' '
249         WRITE(numout,*) ' Namelist parameter for mass conservation checking'
250         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
251         WRITE(numout,*) '    Flag to check mass conservation of NO3/Si/TALK ln_check_mass = ', ln_check_mass
252      ENDIF
253
254   END SUBROUTINE p4z_sms_init
255
256   SUBROUTINE p4z_rst( kt, cdrw )
257      !!---------------------------------------------------------------------
258      !!                   ***  ROUTINE p4z_rst  ***
259      !!
260      !!  ** Purpose : Read or write variables in restart file:
261      !!
262      !!  WRITE(READ) mode:
263      !!       kt        : number of time step since the begining of the experiment at the
264      !!                   end of the current(previous) run
265      !!---------------------------------------------------------------------
266      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
267      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
268      !
269      INTEGER  ::  ji, jj, jk
270      REAL(wp) ::  zcaralk, zbicarb, zco3
271      REAL(wp) ::  ztmas, ztmas1
272      !!---------------------------------------------------------------------
273
274      IF( TRIM(cdrw) == 'READ' ) THEN
275         !
276         IF(lwp) WRITE(numout,*)
277         IF(lwp) WRITE(numout,*) ' p4z_rst : Read specific variables from pisces model '
278         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
279         !
280         IF( iom_varid( numrtr, 'PH', ldstop = .FALSE. ) > 0 ) THEN
281            CALL iom_get( numrtr, jpdom_autoglo, 'PH' , hi(:,:,:)  )
282         ELSE
283!            hi(:,:,:) = 1.e-9
284            ! Set PH from  total alkalinity, borat (???), akb3 (???) and ak23 (???)
285            ! --------------------------------------------------------
286            DO jk = 1, jpk
287               DO jj = 1, jpj
288                  DO ji = 1, jpi
289                     ztmas   = tmask(ji,jj,jk)
290                     ztmas1  = 1. - tmask(ji,jj,jk)
291                     zcaralk = trn(ji,jj,jk,jptal) - borat(ji,jj,jk) / (  1. + 1.E-8 / ( rtrn + akb3(ji,jj,jk) )  )
292                     zco3    = ( zcaralk - trn(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1
293                     zbicarb = ( 2. * trn(ji,jj,jk,jpdic) - zcaralk )
294                     hi(ji,jj,jk) = ( ak23(ji,jj,jk) * zbicarb / zco3 ) * ztmas + 1.e-9 * ztmas1
295                  END DO
296               END DO
297            END DO
298         ENDIF
299         CALL iom_get( numrtr, jpdom_autoglo, 'Silicalim', xksi(:,:) )
300         IF( iom_varid( numrtr, 'Silicamax', ldstop = .FALSE. ) > 0 ) THEN
301            CALL iom_get( numrtr, jpdom_autoglo, 'Silicamax' , xksimax(:,:)  )
302         ELSE
303            xksimax(:,:) = xksi(:,:)
304         ENDIF
305         !
306      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
307         IF( kt == nitrst ) THEN
308            IF(lwp) WRITE(numout,*)
309            IF(lwp) WRITE(numout,*) 'p4z_rst : write pisces restart file  kt =', kt
310            IF(lwp) WRITE(numout,*) '~~~~~~~'
311         ENDIF
312         CALL iom_rstput( kt, nitrst, numrtw, 'PH', hi(:,:,:) )
313         CALL iom_rstput( kt, nitrst, numrtw, 'Silicalim', xksi(:,:) )
314         CALL iom_rstput( kt, nitrst, numrtw, 'Silicamax', xksimax(:,:) )
315      ENDIF
316      !
317   END SUBROUTINE p4z_rst
318
319   SUBROUTINE p4z_dmp( kt )
320      !!----------------------------------------------------------------------
321      !!                    ***  p4z_dmp  ***
322      !!
323      !! ** purpose  : Relaxation of some tracers
324      !!----------------------------------------------------------------------
325      !
326      INTEGER, INTENT( in )  ::     kt ! time step
327      !
328      REAL(wp) ::  alkmean = 2426.     ! mean value of alkalinity ( Glodap ; for Goyet 2391. )
329      REAL(wp) ::  po4mean = 2.165     ! mean value of phosphates
330      REAL(wp) ::  no3mean = 30.90     ! mean value of nitrate
331      REAL(wp) ::  silmean = 91.51     ! mean value of silicate
332      !
333      REAL(wp) :: zarea, zalksum, zpo4sum, zno3sum, zsilsum
334      !!---------------------------------------------------------------------
335
336
337      IF(lwp)  WRITE(numout,*)
338      IF(lwp)  WRITE(numout,*) ' p4z_dmp : Restoring of nutrients at time-step kt = ', kt
339      IF(lwp)  WRITE(numout,*)
340
341      IF( cp_cfg == "orca" .AND. .NOT. lk_c1d ) THEN      ! ORCA configuration (not 1D) !
342         !                                                    ! --------------------------- !
343         ! set total alkalinity, phosphate, nitrate & silicate
344         zarea          = 1._wp / glob_sum( cvol(:,:,:) ) * 1e6             
345
346         zalksum = glob_sum( trn(:,:,:,jptal) * cvol(:,:,:)  ) * zarea
347         zpo4sum = glob_sum( trn(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea * po4r
348         zno3sum = glob_sum( trn(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea * rno3
349         zsilsum = glob_sum( trn(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea
350 
351         IF(lwp) WRITE(numout,*) '       TALK mean : ', zalksum
352         trn(:,:,:,jptal) = trn(:,:,:,jptal) * alkmean / zalksum
353
354         IF(lwp) WRITE(numout,*) '       PO4  mean : ', zpo4sum
355         trn(:,:,:,jppo4) = trn(:,:,:,jppo4) * po4mean / zpo4sum
356
357         IF(lwp) WRITE(numout,*) '       NO3  mean : ', zno3sum
358         trn(:,:,:,jpno3) = trn(:,:,:,jpno3) * no3mean / zno3sum
359
360         IF(lwp) WRITE(numout,*) '       SiO3 mean : ', zsilsum
361         trn(:,:,:,jpsil) = MIN( 400.e-6,trn(:,:,:,jpsil) * silmean / zsilsum )
362         !
363      ENDIF
364
365   END SUBROUTINE p4z_dmp
366
367
368   SUBROUTINE p4z_chk_mass( kt )
369      !!----------------------------------------------------------------------
370      !!                  ***  ROUTINE p4z_chk_mass  ***
371      !!
372      !! ** Purpose :  Mass conservation check
373      !!
374      !!---------------------------------------------------------------------
375      !
376      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index     
377      !!
378      !!---------------------------------------------------------------------
379
380      IF( kt == nittrc000 ) THEN
381         IF( ln_check_mass .AND. lwp) THEN      !   Open budget file of NO3, ALK, Si, Fer
382            CALL ctl_opn( numco2, 'carbon.budget'  , 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
383            CALL ctl_opn( numnut, 'nutrient.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
384         ENDIF
385      ENDIF
386
387      IF( ln_check_mass .AND. kt == nitend ) THEN      !   Compute the budget of NO3, ALK, Si, Fer
388         no3budget = glob_sum( (   trn(:,:,:,jpno3) + trn(:,:,:,jpnh4)  &
389            &                    + trn(:,:,:,jpphy) + trn(:,:,:,jpdia)  &
390            &                    + trn(:,:,:,jpzoo) + trn(:,:,:,jpmes)  &
391            &                    + trn(:,:,:,jppoc)                     &
392#if ! defined key_kriest
393            &                    + trn(:,:,:,jpgoc)                     &
394#endif
395            &                    + trn(:,:,:,jpdoc)                     ) * cvol(:,:,:)  ) 
396         !
397         silbudget = glob_sum( (   trn(:,:,:,jpsil) + trn(:,:,:,jpgsi)  &
398            &                    + trn(:,:,:,jpdsi)                     ) * cvol(:,:,:)  )
399         !
400         alkbudget = glob_sum( (   trn(:,:,:,jpno3) * rno3              &
401            &                    + trn(:,:,:,jptal)                     &
402            &                    + trn(:,:,:,jpcal) * 2.                ) * cvol(:,:,:)  )
403         !
404         ferbudget = glob_sum( (   trn(:,:,:,jpfer) + trn(:,:,:,jpnfe)  &
405            &                    + trn(:,:,:,jpdfe)                     &
406#if ! defined key_kriest
407            &                    + trn(:,:,:,jpbfe)                     &
408#endif
409            &                    + trn(:,:,:,jpsfe)                     &
410            &                    + trn(:,:,:,jpzoo)                     &
411            &                    + trn(:,:,:,jpmes) * ferat3            ) * cvol(:,:,:)  )
412
413         !
414         t_atm_co2_flx  = t_atm_co2_flx / glob_sum( e1e2t(:,:) )
415         t_oce_co2_flx  = t_oce_co2_flx         * 12. / 1.e15 * (-1 )
416         tpp            = tpp           * 1000. * 12. / 1.E15
417         t_oce_co2_exp  = t_oce_co2_exp * 1000. * 12. / 1.E15
418         !
419         no3budget = no3budget / areatot
420         silbudget = silbudget / areatot
421         alkbudget = alkbudget / areatot
422         ferbudget = ferbudget / areatot
423         !
424         IF(lwp) THEN
425            WRITE(numco2,9000) ndastp, t_atm_co2_flx, t_oce_co2_flx, tpp, t_oce_co2_exp
426            WRITE(numnut,9500) ndastp, alkbudget, no3budget, silbudget, ferbudget
427         ENDIF
428         !
429      ENDIF
430       !
431 9000  FORMAT(i8,f10.5,e18.10,f10.5,f10.5)
432 9500  FORMAT(i8,4e18.10)     
433       !
434   END SUBROUTINE p4z_chk_mass
435
436#else
437   !!======================================================================
438   !!  Dummy module :                                   No PISCES bio-model
439   !!======================================================================
440CONTAINS
441   SUBROUTINE p4z_sms( kt )                   ! Empty routine
442      INTEGER, INTENT( in ) ::   kt
443      WRITE(*,*) 'p4z_sms: You should not have seen this print! error?', kt
444   END SUBROUTINE p4z_sms
445#endif 
446
447   !!======================================================================
448END MODULE p4zsms 
Note: See TracBrowser for help on using the repository browser.