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

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

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