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.
trcrst.F90 in trunk/NEMO/TOP_SRC – NEMO

source: trunk/NEMO/TOP_SRC/trcrst.F90 @ 1284

Last change on this file since 1284 was 1254, checked in by cetlod, 15 years ago

update parameter files to take into account the new C14 bomb tracer model, see ticket:298
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 17.0 KB
Line 
1MODULE trcrst
2   !!======================================================================
3   !!                       ***  MODULE trcrst  ***
4   !! TOP :   create, write, read the restart files for passive tracers
5   !!======================================================================
6   !! History :   1.0  !  2007-02 (C. Ethe) adaptation from the ocean
7   !!----------------------------------------------------------------------
8#if defined key_top
9   !!----------------------------------------------------------------------
10   !!   'key_top'                                                TOP models
11   !!----------------------------------------------------------------------
12   !!   trc_rst_opn    : open  restart file
13   !!   trc_rst_read   : read  restart file
14   !!   trc_rst_wri    : write restart file
15   !!----------------------------------------------------------------------
16   USE oce_trc
17   USE trc
18   USE sms_lobster         ! LOBSTER variables
19   USE sms_pisces          ! PISCES variables
20   USE trcsms_cfc          ! CFC variables
21   USE trcsms_c14b         ! C14 variables
22   USE trcsms_my_trc       ! MY_TRC variables
23   USE trctrp_lec   
24   USE lib_mpp
25   USE iom
26   
27   IMPLICIT NONE
28   PRIVATE
29   
30   PUBLIC   trc_rst_opn       ! called by ???
31   PUBLIC   trc_rst_read      ! called by ???
32   PUBLIC   trc_rst_wri       ! called by ???
33   
34   INTEGER, PUBLIC ::   numrtr, numrtw   !: logical unit for trc restart (read and write)
35
36
37   !! * Substitutions
38#  include "top_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/TOP 1.0 , LOCEAN-IPSL (2005)
41   !! $Id$
42   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44   
45CONTAINS
46   
47   SUBROUTINE trc_rst_opn( kt )
48      !!----------------------------------------------------------------------
49      !!                    ***  trc_rst_opn  ***
50      !!
51      !! ** purpose  :   output of sea-trc variable in a netcdf file
52      !!----------------------------------------------------------------------
53      INTEGER, INTENT(in) ::   kt       ! number of iteration
54      !
55      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step define as a character
56      CHARACTER(LEN=50)   ::   clname   ! trc output restart file name
57      !!----------------------------------------------------------------------
58      !
59      IF( kt == nit000 ) lrst_trc = .FALSE. 
60# if defined key_off_tra
61      IF( MOD( kt - 1, nstock ) == 0 ) THEN
62         ! we use kt - 1 and not kt - nit000 to keep the same periodicity from the beginning of the experiment
63         nitrst = kt + nstock - 1                  ! define the next value of nitrst for restart writing
64         IF( nitrst > nitend )   nitrst = nitend   ! make sure we write a restart at the end of the run
65      ENDIF
66#endif
67     
68      IF( kt == nitrst - ndttrc .OR. nitend - nit000 + 1 < 2 * ndttrc ) THEN
69         ! beware if model runs less than 2*ndttrc time step
70         ! beware of the format used to write kt (default is i8.8, that should be large enough)
71         IF( nitrst > 1.0e9 ) THEN   ;   WRITE(clkt,*       ) nitrst
72         ELSE                        ;   WRITE(clkt,'(i8.8)') nitrst
73         ENDIF
74         ! create the file
75         IF(lwp) WRITE(numout,*)
76         clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_trcrst_out)
77         IF(lwp) WRITE(numout,*) '             open trc restart.output NetCDF file: '//clname
78         CALL iom_open( clname, numrtw, ldwrt = .TRUE., kiolib = jprstlib )
79         lrst_trc = .TRUE.
80      ENDIF
81      !
82   END SUBROUTINE trc_rst_opn
83
84
85   SUBROUTINE trc_rst_read 
86      !!----------------------------------------------------------------------
87      !!                    ***  trc_rst_opn  ***
88      !!
89      !! ** purpose  :   read passive tracer fields in restart files
90      !!----------------------------------------------------------------------
91      INTEGER  ::   jn 
92      INTEGER  ::   iarak0
93      REAL(wp) ::   zkt, zarak0
94      !!----------------------------------------------------------------------
95
96      IF(lwp) WRITE(numout,*)
97      IF(lwp) WRITE(numout,*) 'trc_rst_read : read the TOP restart file'
98      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
99
100      IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd ) THEN   ;   iarak0 = 1
101      ELSE                                           ;   iarak0 = 0
102      ENDIF
103
104      IF(lwp) WRITE(numout,*)
105      IF(lwp) WRITE(numout,*) '   the present run starts at the time step nit000 = ', nit000
106      IF(lwp .AND. iarak0 == 1 )   WRITE(numout,*) '   and needs previous fields for Arakawa sheme ??? '
107
108
109      ! Time domain : restart
110      ! -------------------------
111      IF(lwp) WRITE(numout,*)
112      IF(lwp) WRITE(numout,*) ' *** passive tracer restart option'
113      SELECT CASE ( nrsttr )
114      CASE ( 0 )
115         IF(lwp) WRITE(numout,*) '    nrsttr = 0 no control of nit000'
116      CASE ( 1 )
117         IF(lwp) WRITE(numout,*) '    nrsttr = 1 we control the date of nit000'
118      CASE ( 2 )
119         IF(lwp) WRITE(numout,*) '    nrsttr = 2 the date adatrj is read in restart file'
120      CASE DEFAULT
121         IF(lwp) WRITE(numout,*) '  ===>>>> nrsttr not equal 0, 1 or 2 : no control of the date'
122         IF(lwp) WRITE(numout,*) '  =======                  ========='
123      END SELECT
124
125      CALL iom_open( cn_trcrst_in, numrtr, kiolib = jprstlib )
126
127      CALL iom_get( numrtr, 'kt'   , zkt    )
128      CALL iom_get( numrtr, 'arak0', zarak0 )
129
130      IF(lwp) WRITE(numout,*)
131      IF(lwp) WRITE(numout,*) ' Info on the restart file read : '
132      IF(lwp) WRITE(numout,*) '    time-step           : ', NINT( zkt    )
133      IF(lwp) WRITE(numout,*) '    arakawa option      : ', NINT( zarak0 )
134
135
136      IF( nittrc000 - NINT( zkt ) /= 1 .AND. nrsttr /= 0 )  &      ! control of date
137         &   CALL ctl_stop( ' ===>>>> : problem with nit000 for the restart', &
138         &                  ' verify the restart file or rerun with nrstdt = 0 (namelist)' )
139
140      IF( iarak0 /= NINT( zarak0 ) )   &                           ! Control of the scheme
141         & CALL ctl_stop( ' ===>>>> : problem with advection scheme', &
142         & ' it must be the same type for both restart and previous run', &
143         & ' centered or euler '  )
144
145
146      ! READ prognostic variables and computes diagnostic variable
147      DO jn = 1, jptra
148         CALL iom_get( numrtr, jpdom_autoglo, 'TRN'//ctrcnm(jn), trn(:,:,:,jn) ) 
149      END DO
150
151     DO jn = 1, jptra
152         CALL iom_get( numrtr, jpdom_autoglo, 'TRB'//ctrcnm(jn), trb(:,:,:,jn) ) 
153      END DO
154
155#if defined key_lobster
156      CALL iom_get( numrtr, jpdom_autoglo, 'SEDB'//ctrcnm(jpdet), sedpocb(:,:) ) 
157      CALL iom_get( numrtr, jpdom_autoglo, 'SEDN'//ctrcnm(jpdet), sedpocn(:,:) ) 
158#endif
159
160#if defined key_pisces
161      CALL iom_get( numrtr, jpdom_autoglo, 'Silicalim', xksi(:,:) ) 
162      CALL iom_get( numrtr, jpdom_autoglo, 'Silicamax', xksimax(:,:) )
163      CALL trc_rst_ini  ! Initialisation of some variables
164#endif
165
166#if defined key_cfc
167      DO jn = jp_cfc0, jp_cfc1
168         CALL iom_get( numrtr, jpdom_autoglo, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jn) ) 
169      END DO
170#endif
171
172#if defined key_c14b
173      CALL iom_get( numrtr, jpdom_autoglo, 'qint_'//ctrcnm(jn) , qint_c14(:,:) ) 
174#endif
175
176#if defined key_my_trc
177#endif
178
179
180      CALL iom_close( numrtr )
181      !
182   END SUBROUTINE trc_rst_read
183
184
185   SUBROUTINE trc_rst_wri( kt )
186      !!----------------------------------------------------------------------
187      !!                    ***  trc_rst_wri  ***
188      !!
189      !! ** purpose  :   write passive tracer fields in restart files
190      !!----------------------------------------------------------------------
191      INTEGER, INTENT( in ) ::   kt    ! ocean time-step index
192      !!
193      INTEGER  :: ji, jj, jk, jn
194      !!----------------------------------------------------------------------
195
196      IF(  kt == nitrst ) THEN
197         IF(lwp) WRITE(numout,*)
198         IF(lwp) WRITE(numout,*) 'trc_wri : write the TOP restart file (NetCDF) at it= ', kt, ' date= ', ndastp
199         IF(lwp) WRITE(numout,*) '~~~~~~~'
200      ENDIF
201
202      CALL iom_rstput( kt, nitrst, numrtw, 'kt'   ,  REAL( kt, wp )  )
203      IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd ) THEN
204         CALL iom_rstput( kt, nitrst, numrtw, 'arak0', 1. )
205      ELSE
206         CALL iom_rstput( kt, nitrst, numrtw, 'arak0', 0. )
207      ENDIF
208
209
210         ! prognostic variables
211         ! --------------------
212      DO jn = 1, jptra
213         CALL iom_rstput( kt, nitrst, numrtw, 'TRN'//ctrcnm(jn), trn(:,:,:,jn) )
214      END DO
215
216      DO jn = 1, jptra
217         CALL iom_rstput( kt, nitrst, numrtw, 'TRB'//ctrcnm(jn), trb(:,:,:,jn) )
218      END DO
219
220#if defined key_lobster
221      CALL iom_rstput( kt, nitrst, numrtw, 'SEDB'//ctrcnm(jpdet), sedpocb(:,:) )
222      CALL iom_rstput( kt, nitrst, numrtw, 'SEDN'//ctrcnm(jpdet), sedpocn(:,:) )
223#endif
224
225#if defined key_pisces
226      CALL iom_rstput( kt, nitrst, numrtw, 'Silicalim', xksi(:,:) )
227      CALL iom_rstput( kt, nitrst, numrtw, 'Silicamax', xksimax(:,:) )
228#endif
229
230#if defined key_cfc
231      DO jn = jp_cfc0, jp_cfc1
232         CALL iom_rstput( kt, nitrst, numrtw, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jn) )
233      END DO
234#endif
235
236#if defined key_c14b
237      CALL iom_rstput( kt, nitrst, numrtw, 'qint_'//ctrcnm(jn), qint_c14(:,:) )
238#endif
239
240#if defined key_my_trc
241#endif
242
243
244       IF( kt == nitrst ) THEN
245          CALL trc_rst_stat            ! statistics
246          CALL iom_close( numrtw )     ! close the restart file (only at last time step)
247#if ! defined key_trdmld_trc
248          lrst_trc = .FALSE.
249#endif
250       ENDIF
251      !
252   END SUBROUTINE trc_rst_wri
253
254# if defined key_pisces 
255
256   SUBROUTINE trc_rst_ini 
257      !!----------------------------------------------------------------------
258      !!                    ***  trc_rst_ini  ***
259      !!
260      !! ** purpose  : Initialisation of some variables ( hi
261      !!----------------------------------------------------------------------
262      INTEGER  :: ji, jj, jk, jn 
263      REAL(wp) ::  &
264         alkmean = 2426. ,  & ! mean value of alkalinity ( Glodap ; for Goyet 2391. )
265         po4mean = 2.165 ,  & ! mean value of phosphates
266         no3mean = 30.90 ,  & ! mean value of nitrate
267         siomean = 91.51      ! mean value of silicate
268     
269      REAL(wp) ::   zvol, ztrasum
270      REAL(wp) ::   caralk, bicarb, co3
271
272      IF(lwp)  WRITE(numout,*)
273
274      IF( cp_cfg == "orca" .AND. .NOT. lk_trc_c1d ) THEN      ! ORCA condiguration (not 1D) !
275         !                                                    ! --------------------------- !
276         ! set total alkalinity, phosphate, NO3 & silicate
277
278         ! total alkalinity
279         ztrasum = 0.e0             
280         DO jk = 1, jpk
281            DO jj = 1, jpj
282               DO ji = 1, jpi
283                  zvol = cvol(ji,jj,jk)
284#  if defined key_off_degrad
285                  zvol = zvol * facvol(ji,jj,jk)
286#  endif
287                  ztrasum = ztrasum + trn(ji,jj,jk,jptal) * zvol
288               END DO
289            END DO
290         END DO
291         IF( lk_mpp )   CALL mpp_sum( ztrasum )     ! sum over the global domain 
292         
293         ztrasum = ztrasum / areatot * 1.e6
294         IF(lwp) WRITE(numout,*) '       TALK mean : ', ztrasum
295         trn(:,:,:,jptal) = trn(:,:,:,jptal) * alkmean / ztrasum
296           
297         ! phosphate
298         ztrasum = 0.e0
299         DO jk = 1, jpk
300            DO jj = 1, jpj
301               DO ji = 1, jpi
302                  zvol = cvol(ji,jj,jk)
303#  if defined key_off_degrad
304                  zvol = zvol * facvol(ji,jj,jk)
305#  endif
306                  ztrasum = ztrasum + trn(ji,jj,jk,jppo4) * zvol
307               END DO
308            END DO
309         END DO
310         IF( lk_mpp )   CALL mpp_sum( ztrasum )     ! sum over the global domain 
311         
312         ztrasum = ztrasum / areatot * 1.e6 / 122.
313         IF(lwp) WRITE(numout,*) '       PO4  mean : ', ztrasum
314         trn(:,:,:,jppo4) = trn(:,:,:,jppo4) * po4mean / ztrasum
315       
316         ! Nitrates         
317         ztrasum = 0.e0
318         DO jk = 1, jpk
319            DO jj = 1, jpj
320               DO ji = 1, jpi
321                  zvol = cvol(ji,jj,jk)
322#  if defined key_off_degrad
323                  zvol = zvol * facvol(ji,jj,jk)
324#  endif
325                  ztrasum = ztrasum + trn(ji,jj,jk,jpno3) * zvol
326               END DO
327            END DO
328         END DO
329         IF( lk_mpp )   CALL mpp_sum( ztrasum )     ! sum over the global domain 
330         
331         ztrasum = ztrasum / areatot * 1.e6 / 7.6
332         IF(lwp) WRITE(numout,*) '       NO3  mean : ', ztrasum
333         trn(:,:,:,jpno3) = trn(:,:,:,jpno3) * no3mean / ztrasum
334         
335         ! Silicate
336         ztrasum = 0.e0
337         DO jk = 1, jpk
338            DO jj = 1, jpj
339               DO ji = 1, jpi
340                  zvol = cvol(ji,jj,jk)
341#  if defined key_off_degrad
342                  zvol = zvol * facvol(ji,jj,jk)
343#  endif
344                  ztrasum = ztrasum + trn(ji,jj,jk,jpsil) * zvol
345               END DO
346            END DO
347         END DO
348         IF( lk_mpp )   CALL mpp_sum( ztrasum )     ! sum over the global domain 
349         ztrasum = ztrasum / areatot * 1.e6
350         IF(lwp) WRITE(numout,*) '       SiO3 mean : ', ztrasum
351         trn(:,:,:,jpsil) = MIN( 400.e-6,trn(:,:,:,jpsil) * siomean / ztrasum ) 
352         !
353      ENDIF
354
355!#if defined key_kriest
356!     !! Initialize number of particles from a standart restart file
357!     !! The name of big organic particles jpgoc has been only change
358!     !! and replace by jpnum but the values here are concentration
359!     trn(:,:,:,jppoc) = trn(:,:,:,jppoc) + trn(:,:,:,jpnum)
360!     trn(:,:,:,jpnum) = trn(:,:,:,jppoc) / ( 6. * xkr_massp )
361!#endif
362      !!  Set hi (???) from  total alcalinity, borat (???), akb3 (???) and ak23 (???)
363      !!  ---------------------------------------------------------------------
364      DO jk = 1, jpk
365         DO jj = 1, jpj
366            DO ji = 1,jpi
367               caralk = trn(ji,jj,jk,jptal) - borat(ji,jj,jk) / ( 1. + 1.e-8 / ( rtrn + akb3(ji,jj,jk) ) )
368               co3    = ( caralk - trn(ji,jj,jk,jpdic) ) *       tmask(ji,jj,jk)   &
369                  &   +             0.5e-3               * ( 1.- tmask(ji,jj,jk) )
370               bicarb = 2.* trn(ji,jj,jk,jpdic) - caralk
371               hi(ji,jj,jk) = ( ak23(ji,jj,jk) * bicarb / co3 )   *       tmask(ji,jj,jk)   &
372                  &         +             1.0e-9                  * ( 1.- tmask(ji,jj,jk) )
373            END DO
374         END DO
375      END DO
376
377   END SUBROUTINE trc_rst_ini
378
379#endif
380      !!----------------------------------------------------------------------
381
382   SUBROUTINE trc_rst_stat
383      !!----------------------------------------------------------------------
384      !!                    ***  trc_rst_stat  ***
385      !!
386      !! ** purpose  :   Compute tracers statistics
387      !!----------------------------------------------------------------------
388
389      INTEGER  :: ji, jj, jk, jn
390      REAL(wp) :: zdiag_var, zdiag_varmin, zdiag_varmax, zdiag_tot
391      REAL(wp) :: zder, zvol
392      !!----------------------------------------------------------------------
393
394
395      IF( lwp ) THEN
396         WRITE(numout,*) 
397         WRITE(numout,*) '           ----TRACER STAT----             '
398         WRITE(numout,*) 
399      ENDIF
400     
401      zdiag_tot = 0.e0
402      DO jn = 1, jptra
403         zdiag_var    = 0.e0
404         zdiag_varmin = 0.e0
405         zdiag_varmax = 0.e0
406         DO ji = 1, jpi
407            DO jj = 1, jpj
408               DO jk = 1,jpk
409                  zvol = cvol(ji,jj,jk)
410#  if defined key_off_degrad
411                  zvol = zvol * facvol(ji,jj,jk)
412#  endif
413                  zdiag_var = zdiag_var + trn(ji,jj,jk,jn) * zvol
414               END DO
415            END DO
416         END DO
417         
418         zdiag_varmin = MINVAL( trn(:,:,:,jn), mask= ((tmask*SPREAD(tmask_i,DIM=3,NCOPIES=jpk).NE.0.)) )
419         zdiag_varmax = MAXVAL( trn(:,:,:,jn), mask= ((tmask*SPREAD(tmask_i,DIM=3,NCOPIES=jpk).NE.0.)) )
420         IF( lk_mpp ) THEN
421            CALL mpp_min( zdiag_varmin )      ! min over the global domain
422            CALL mpp_max( zdiag_varmax )      ! max over the global domain
423            CALL mpp_sum( zdiag_var    )      ! sum over the global domain
424         END IF
425         zdiag_tot = zdiag_tot + zdiag_var
426         zdiag_var = zdiag_var / areatot
427         IF(lwp) WRITE(numout,*) '   MEAN NO ', jn, ctrcnm(jn), ' = ', zdiag_var,   &
428            &                    ' MIN = ', zdiag_varmin, ' MAX = ', zdiag_varmax
429      END DO
430     
431      zder = ( ( zdiag_tot - trai ) / ( trai + 1.e-12 )  ) * 100._wp
432      IF(lwp) WRITE(numout,*) '   Integral of all tracers over the full domain  = ', zdiag_tot
433      IF(lwp) WRITE(numout,*) '   Drift of the sum of all tracers =', zder, ' %'
434     
435   END SUBROUTINE trc_rst_stat
436
437#else
438   !!----------------------------------------------------------------------
439   !!  Dummy module :                                     No passive tracer
440   !!----------------------------------------------------------------------
441CONTAINS
442   SUBROUTINE trc_rst_read                      ! Empty routines
443   END SUBROUTINE trc_rst_read
444   SUBROUTINE trc_rst_wri( kt )
445      INTEGER, INTENT ( in ) :: kt
446      WRITE(*,*) 'trc_rst_wri: You should not have seen this print! error?', kt
447   END SUBROUTINE trc_rst_wri   
448#endif
449
450   !!======================================================================
451END MODULE trcrst
Note: See TracBrowser for help on using the repository browser.