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

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

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