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

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

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