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

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

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