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

Last change on this file since 1177 was 1177, checked in by cetlod, 16 years ago

update modules to take into account new trends organization, see ticket:248

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