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 branches/dev_001_GM/NEMO/TOP_SRC – NEMO

source: branches/dev_001_GM/NEMO/TOP_SRC/trcrst.F90 @ 771

Last change on this file since 771 was 771, checked in by gm, 16 years ago

dev_001_GM - small error corrections
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.8 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_passivetrc
9   !!----------------------------------------------------------------------
10   !!   'key_passivetrc'                                    Passive tracers
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
19   USE trccfc          ! CFC variables
20   USE trctrp_lec   
21   USE lib_mpp
22   USE iom
23   
24   IMPLICIT NONE
25   PRIVATE
26   
27   PUBLIC   trc_rst_opn       ! called by ???
28   PUBLIC   trc_rst_read      ! called by ???
29   PUBLIC   trc_rst_wri       ! called by ???
30   
31   LOGICAL, PUBLIC ::   lrst_trc         !: logical to control the trc restart write
32   INTEGER, PUBLIC ::   numrtr, numrtw   !: logical unit for trc restart (read and write)
33
34   !! * Substitutions
35#  include "passivetrc_substitute.h90"
36   !!----------------------------------------------------------------------
37   !! NEMO/TOP 1.0 , LOCEAN-IPSL (2005)
38   !! $Id$
39   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
40   !!----------------------------------------------------------------------
41   
42CONTAINS
43   
44   SUBROUTINE trc_rst_opn( kt )
45      !!----------------------------------------------------------------------
46      !!                    ***  trc_rst_opn  ***
47      !!
48      !! ** purpose  :   output of sea-trc variable in a netcdf file
49      !!----------------------------------------------------------------------
50      INTEGER, INTENT(in) ::   kt       ! number of iteration
51      !
52      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step define as a character
53      CHARACTER(LEN=50)   ::   clname   ! trc output restart file name
54      !!----------------------------------------------------------------------
55      !
56      IF( kt == nit000 )  THEN
57         lrst_trc = .FALSE.
58# if defined key_off_tra
59         nitrst = nitend  ! in online version, already done in rst_opn routine defined in restart.F90 module
60# endif
61      ENDIF
62     
63      IF( kt == nitrst - ndttrc .OR. nitend - nit000 + 1 < 2 * ndttrc ) THEN
64         ! beware if model runs less than 2*ndttrc time step
65         ! beware of the format used to write kt (default is i8.8, that should be large enough)
66         IF( nitrst > 1.0e9 ) THEN   ;   WRITE(clkt,*       ) nitrst
67         ELSE                        ;   WRITE(clkt,'(i8.8)') nitrst
68         ENDIF
69         ! create the file
70         IF(lwp) WRITE(numout,*)
71         clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_trc"
72         IF(lwp) WRITE(numout,*) '             open trc restart.output NetCDF file: '//clname
73         CALL iom_open( clname, numrtw, ldwrt = .TRUE., kiolib = jprstlib )
74         lrst_trc = .TRUE.
75      ENDIF
76      !
77   END SUBROUTINE trc_rst_opn
78
79
80   SUBROUTINE trc_rst_read 
81      !!----------------------------------------------------------------------
82      !!                    ***  trc_rst_opn  ***
83      !!
84      !! ** purpose  :   read passive tracer fields in restart files
85      !!----------------------------------------------------------------------
86      INTEGER  ::   jn 
87      INTEGER  ::   iarak0
88      REAL(wp) ::   zkt, zarak0
89# if defined key_trc_pisces 
90      REAL(wp) ::   ztrasum
91      INTEGER  ::   ji, jj, jk
92      REAL(wp) ::   caralk, bicarb, co3
93# endif
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( 'restart.trc', 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_local, 'TRN'//ctrcnm(jn), trn(:,:,:,jn) ) 
149         CALL iom_get( numrtr, jpdom_local, 'TRB'//ctrcnm(jn), trb(:,:,:,jn) ) 
150      END DO
151# if defined key_trc_lobster1
152      CALL iom_get( numrtr, jpdom_local, 'SEDB'//ctrcnm(jpdet), sedpocb(:,:) ) 
153      CALL iom_get( numrtr, jpdom_local, 'SEDN'//ctrcnm(jpdet), sedpocn(:,:) ) 
154# elif defined key_trc_pisces
155      CALL iom_get( numrtr, jpdom_local, 'Silicalim', xksi(:,:) ) 
156      xksimax = xksi
157# elif defined key_cfc
158      DO jn = 1, jptra
159         CALL iom_get( numrtr, jpdom_local, 'qint'//ctrcnm(jn), qint(:,:,jn) ) 
160         CALL iom_get( numrtr, jpdom_local, 'qtr'//ctrcnm(jn) , qtr( :,:,jn) ) 
161      END DO
162# endif
163
164# if defined key_trc_pisces 
165      !                                                         ! --------------------------- !
166      IF( cp_cfg == "orca" .AND. .NOT. lk_trccfg_1d ) THEN      ! ORCA condiguration (not 1D) !
167         !                                                      ! --------------------------- !
168         !                                                      ! set total alkalinity, phosphate, NO3 & silicate
169         !                          ! total alkalinity
170         ztrasum = 0.e0             ! ----------------
171         DO jk = 1, jpk
172            DO jj = 1, jpj
173               DO ji = 1, jpi
174                  ztrasum = ztrasum + trn(ji,jj,jk,jptal) * tmask(ji,jj,jk) * tmask_i(ji,jj)  &
175#  if defined key_off_degrad
176                     &              * facvol(ji,jj,jk)   &
177#  endif
178                     &              * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
179               END DO
180            END DO
181         END DO
182         IF( lk_mpp )   CALL mpp_sum( ztrasum )     ! sum over the global domain 
183
184         IF(lwp) WRITE(numout,*) 'TALK moyen ', ztrasum / areatot * 1.e6
185         ztrasum = ztrasum / areatot * 1.e6
186         trn(:,:,:,jptal) = trn(:,:,:,jptal) * 2391. / ztrasum
187
188         !                          ! phosphate
189         ztrasum = 0.e0             ! ---------
190         DO jk = 1, jpk
191            DO jj = 1, jpj
192               DO ji = 1, jpi
193                  ztrasum = ztrasum + trn(ji,jj,jk,jppo4) * tmask(ji,jj,jk) * tmask_i(ji,jj)  &
194#  if defined key_off_degrad
195                     &              * facvol(ji,jj,jk)   &
196#  endif
197                     &              * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
198               END DO
199            END DO
200         END DO
201         IF( lk_mpp )   CALL mpp_sum( ztrasum )     ! sum over the global domain 
202
203         IF(lwp) WRITE(numout,*) 'PO4 moyen ', ztrasum/areatot*1E6/122.
204         ztrasum = ztrasum / areatot * 1.e6 / 122.
205         trn(:,:,:,jppo4) = trn(:,:,:,jppo4) * 2.165 / ztrasum
206
207         !                          ! NO3
208         ztrasum = 0.e0             ! ---
209         DO jk = 1, jpk
210            DO jj = 1, jpj
211               DO ji = 1, jpi
212                  ztrasum = ztrasum + trn(ji,jj,jk,jpno3) * tmask(ji,jj,jk) * tmask_i(ji,jj)  &
213#  if defined key_off_degrad
214                     &              * facvol(ji,jj,jk)   &
215#  endif
216                     &              * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
217               END DO
218            END DO
219         END DO
220         IF( lk_mpp )   CALL mpp_sum( ztrasum )     ! sum over the global domain 
221
222         IF(lwp) WRITE(numout,*) 'NO3 moyen ', ztrasum / areatot*1.e6 / 7.6
223         ztrasum = ztrasum / areatot * 1.e6 / 7.6
224         trn(:,:,:,jpno3) = trn(:,:,:,jpno3) * 30.9 / ztrasum
225
226         !                          ! Silicate
227         ztrasum = 0.e0             ! --------
228         DO jk = 1, jpk
229            DO jj = 1, jpj
230               DO ji = 1, jpi
231                  ztrasum = ztrasum + trn(ji,jj,jk,jpsil) * tmask(ji,jj,jk) * tmask_i(ji,jj)   &
232#  if defined key_off_degrad
233                     &              * facvol(ji,jj,jk)   &
234#  endif
235                     &              * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
236               END DO
237            END DO
238         END DO
239         IF( lk_mpp )   CALL mpp_sum( ztrasum )     ! sum over the global domain 
240
241         IF(lwp) WRITE(numout,*) 'SiO3 moyen ', ztrasum / areatot * 1.e6
242         ztrasum = ztrasum / areatot * 1.e6
243         trn(:,:,:,jpsil) = MIN( 400.e-6,trn(:,:,:,jpsil) * 91.51 / ztrasum ) 
244         !
245      ENDIF
246
247!#if defined key_trc_kriest
248!      !! Initialize number of particles from a standart restart file
249!      !! The name of big organic particles jpgoc has been only change
250!      !! and replace by jpnum but the values here are concentration
251!      trn(:,:,:,jppoc) = trn(:,:,:,jppoc) + trn(:,:,:,jpnum)
252!      trn(:,:,:,jpnum) = trn(:,:,:,jppoc) / ( 6. * xkr_massp )
253!#endif
254      !!  Set hi (???) from  total alcalinity, borat (???), akb3 (???) and ak23 (???)
255      !!  ---------------------------------------------------------------------
256      DO jk = 1, jpk
257         DO jj = 1, jpj
258            DO ji = 1,jpi
259               caralk = trn(ji,jj,jk,jptal) - borat(ji,jj,jk) / ( 1. + 1.e-8 / ( rtrn + akb3(ji,jj,jk) ) )
260               co3    = ( caralk - trn(ji,jj,jk,jpdic) ) *       tmask(ji,jj,jk)   &
261                  &   +             0.5e-3               * ( 1.- tmask(ji,jj,jk) )
262               bicarb = 2.* trn(ji,jj,jk,jpdic) - caralk
263               hi(ji,jj,jk) = ( ak23(ji,jj,jk) * bicarb / co3 )   *       tmask(ji,jj,jk)   &
264                  &         +             1.0e-9                  * ( 1.- tmask(ji,jj,jk) )
265            END DO
266         END DO
267      END DO
268# endif
269      trb(:,:,:,:) = trn(:,:,:,:)
270
271      CALL iom_close( numrtr )
272      !
273   END SUBROUTINE trc_rst_read
274
275
276   SUBROUTINE trc_rst_wri( kt )
277      !!----------------------------------------------------------------------
278      !!                    ***  trc_rst_wri  ***
279      !!
280      !! ** purpose  :   write passive tracer fields in restart files
281      !!----------------------------------------------------------------------
282      INTEGER, INTENT( in ) ::   kt    ! ocean time-step index
283      !!
284      INTEGER  :: ji, jj, jk, jn
285      REAL(wp) :: zdiag_var, zdiag_varmin, zdiag_varmax, zdiag_tot
286      REAL(wp) :: zder
287      !!----------------------------------------------------------------------
288
289      IF( MOD( kt, nstock ) == 0 .OR. kt == nitend ) THEN
290
291         ! 0. initialisations
292         ! ------------------
293         IF(lwp) WRITE(numout,*)
294         IF(lwp) WRITE(numout,*) 'trc_wri : write the TOP restart file (NetCDF) at it= ', kt, ' date= ', ndastp
295         IF(lwp) WRITE(numout,*) '~~~~~~~'
296
297
298         CALL iom_rstput( kt, nitrst, numrtw, 'kt'   ,  REAL( kt, wp )  )
299
300         IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd ) THEN
301            CALL iom_rstput( kt, nitrst, numrtw, 'arak0', 1. )
302         ELSE
303            CALL iom_rstput( kt, nitrst, numrtw, 'arak0', 0. )
304         ENDIF
305
306
307         ! prognostic variables
308         ! --------------------
309         DO jn = 1, jptra
310            CALL iom_rstput( kt, nitrst, numrtw, 'TRN'//ctrcnm(jn), trn(:,:,:,jn) )
311            CALL iom_rstput( kt, nitrst, numrtw, 'TRB'//ctrcnm(jn), trb(:,:,:,jn) )
312         END DO
313
314#if defined key_trc_lobster1
315         CALL iom_rstput( kt, nitrst, numrtw, 'SEDB'//ctrcnm(jpdet), sedpocb(:,:) )
316         CALL iom_rstput( kt, nitrst, numrtw, 'SEDN'//ctrcnm(jpdet), sedpocn(:,:) )
317#elif defined key_trc_pisces
318         CALL iom_rstput( kt, nitrst, numrtw, 'Silicalim', xksi(:,:) )
319
320#elif defined key_cfc
321         DO jn = 1, jptra
322            CALL iom_rstput( kt, nitrst, numrtw, 'qint'//ctrcnm(jn), qint(:,:,jn) )
323            CALL iom_rstput( kt, nitrst, numrtw, 'qtr'//ctrcnm(jn) , qtr (:,:,jn) )
324         END DO
325#endif
326
327         IF(lwp) WRITE(numout,*) '----TRACER STAT----'
328
329         zdiag_tot = 0.e0
330         DO jn = 1, jptra
331            zdiag_var    = 0.e0
332            zdiag_varmin = 0.e0
333            zdiag_varmax = 0.e0
334            DO ji = 1, jpi
335               DO jj = 1, jpj
336                  DO jk = 1,jpk
337                     zdiag_var = zdiag_var + trn(ji,jj,jk,jn) * tmask(ji,jj,jk) * tmask_i(ji,jj)   &
338#if defined key_off_degrad
339                        &   * facvol(ji,jj,jk)   &
340#endif
341                        &   * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
342                  END DO
343               END DO
344            END DO
345
346            zdiag_varmin = MINVAL( trn(:,:,:,jn), mask= ((tmask*SPREAD(tmask_i,DIM=3,NCOPIES=jpk).NE.0.)) )
347            zdiag_varmax = MAXVAL( trn(:,:,:,jn), mask= ((tmask*SPREAD(tmask_i,DIM=3,NCOPIES=jpk).NE.0.)) )
348
349            IF( lk_mpp ) THEN
350               CALL mpp_min( zdiag_varmin )      ! min over the global domain
351               CALL mpp_max( zdiag_varmax )      ! max over the global domain
352               CALL mpp_sum( zdiag_var    )      ! sum over the global domain
353            END IF
354
355            zdiag_tot = zdiag_tot + zdiag_var
356            zdiag_var = zdiag_var / areatot
357
358            IF(lwp) WRITE(numout,*) '   MEAN NO ', jn, ctrcnm(jn), ' = ', zdiag_var,   &
359               &                    ' MIN = ', zdiag_varmin, ' MAX = ', zdiag_varmax
360         END DO
361
362         zder = ( ( zdiag_tot - trai ) / trai ) * 100._wp
363         IF(lwp) WRITE(numout,*) '   Integral of all tracers over the full domain  = ', zdiag_tot
364         IF(lwp) WRITE(numout,*) '   Drift of the sum of all tracers =', zder, ' %'
365
366         CALL iom_close(numrtw)
367         !
368      ENDIF
369      !
370   END SUBROUTINE trc_rst_wri
371
372#else
373   !!----------------------------------------------------------------------
374   !!  Dummy module :                                     No passive tracer
375   !!----------------------------------------------------------------------
376CONTAINS
377   SUBROUTINE trc_rst_read                      ! Empty routines
378   END SUBROUTINE trc_rst_read
379   SUBROUTINE trc_rst_wri( kt )
380      INTEGER, INTENT ( in ) :: kt
381      WRITE(*,*) 'trc_rst_wri: You should not have seen this print! error?', kt
382   END SUBROUTINE trc_rst_wri   
383#endif
384
385   !!======================================================================
386END MODULE trcrst
Note: See TracBrowser for help on using the repository browser.