source: NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/TOP/trcstp.F90 @ 12397

Last change on this file since 12397 was 12397, checked in by davestorkey, 10 months ago

2020/KERNEL-03_Storkey_Coward_RK3_stage2 : Consolidation of code to
handle initial Euler timestep in the context of leapfrog
timestepping. This version passes all SETTE tests but fails to bit
compare with the control for several tests (ORCA2_ICE_PISCES, AMM12,
ISOMIP, AGRIF_DEMO, SPITZ12).

  • Property svn:keywords set to Id
File size: 12.6 KB
Line 
1MODULE trcstp
2   !!======================================================================
3   !!                       ***  MODULE trcstp  ***
4   !! Time-stepping    : time loop of opa for passive tracer
5   !!======================================================================
6   !! History :  1.0  !  2004-03  (C. Ethe)  Original
7   !!            4.1  !  2019-08  (A. Coward, D. Storkey) rewrite in preparation for new timestepping scheme
8   !!----------------------------------------------------------------------
9#if defined key_top
10   !!----------------------------------------------------------------------
11   !!   trc_stp       : passive tracer system time-stepping
12   !!----------------------------------------------------------------------
13   USE oce_trc        ! ocean dynamics and active tracers variables
14   USE sbc_oce
15   USE trc
16   USE trctrp         ! passive tracers transport
17   USE trcsms         ! passive tracers sources and sinks
18   USE trcwri
19   USE trcrst
20   USE trdtrc_oce
21   USE trdmxl_trc
22   USE sms_pisces,  ONLY : ln_check_mass
23   !
24   USE prtctl_trc     ! Print control for debbuging
25   USE iom            !
26   USE in_out_manager !
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   trc_stp    ! called by step
32
33   LOGICAL  ::   llnew                   ! ???
34   REAL(wp) ::   rdt_sampl               ! ???
35   INTEGER  ::   nb_rec_per_day, ktdcy   ! ???
36   REAL(wp) ::   rsecfst, rseclast       ! ???
37   REAL(wp), DIMENSION(:,:,:), SAVE, ALLOCATABLE ::   qsr_arr   ! save qsr during TOP time-step
38
39   !!----------------------------------------------------------------------
40   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
41   !! $Id$
42   !! Software governed by the CeCILL license (see ./LICENSE)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE trc_stp( kt, Kbb, Kmm, Krhs, Kaa )
47      !!-------------------------------------------------------------------
48      !!                     ***  ROUTINE trc_stp  ***
49      !!                     
50      !! ** Purpose :   Time loop of opa for passive tracer
51      !!
52      !! ** Method  :   Compute the passive tracers trends
53      !!                Update the passive tracers
54      !!-------------------------------------------------------------------
55      INTEGER, INTENT( in ) :: kt                  ! ocean time-step index
56      INTEGER, INTENT( in ) :: Kbb, Kmm, Krhs, Kaa ! time level indices
57      !
58      INTEGER ::   jk, jn   ! dummy loop indices
59      REAL(wp)::   ztrai    ! local scalar
60      LOGICAL ::   ll_trcstat ! local logical
61      CHARACTER (len=25) ::   charout   !
62      !!-------------------------------------------------------------------
63      !
64      IF( ln_timing )   CALL timing_start('trc_stp')
65      !
66      IF( l_1st_euler .OR. ln_top_euler ) THEN     ! at nittrc000
67         r2dttrc =  rdttrc           ! = rdttrc (use or restarting with Euler time stepping)
68      ELSEIF( kt <= nittrc000 + 1 ) THEN                                     ! at nittrc000 or nittrc000+1
69         r2dttrc = 2. * rdttrc       ! = 2 rdttrc (leapfrog)
70      ENDIF
71      !
72      ll_trcstat  = ( sn_cfctl%l_trcstat ) .AND. &
73     &              ( ( MOD( kt, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend ) )
74
75      IF( kt == nittrc000 )                      CALL trc_stp_ctl   ! control
76      IF( kt == nittrc000 .AND. lk_trdmxl_trc )  CALL trd_mxl_trc_init    ! trends: Mixed-layer
77      !
78      IF( .NOT.ln_linssh ) THEN                                           ! update ocean volume due to ssh temporal evolution
79         DO jk = 1, jpk
80            cvol(:,:,jk) = e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
81         END DO
82         IF ( ll_trcstat .OR. kt == nitrst .OR. ( ln_check_mass .AND. kt == nitend )   &
83            & .OR. iom_use( "pno3tot" ) .OR. iom_use( "ppo4tot" ) .OR. iom_use( "psiltot" )   &
84            & .OR. iom_use( "palktot" ) .OR. iom_use( "pfertot" ) )                           &
85            &     areatot = glob_sum( 'trcstp', cvol(:,:,:) )
86      ENDIF
87      !
88      IF( l_trcdm2dc )   CALL trc_mean_qsr( kt )
89      !   
90      !
91      IF(sn_cfctl%l_prttrc) THEN
92         WRITE(charout,FMT="('kt =', I4,'  d/m/y =',I2,I2,I4)") kt, nday, nmonth, nyear
93         CALL prt_ctl_trc_info(charout)
94      ENDIF
95      !
96      tr(:,:,:,:,Krhs) = 0._wp
97      !
98      CALL trc_rst_opn  ( kt )                            ! Open tracer restart file
99      IF( lrst_trc )  CALL trc_rst_cal  ( kt, 'WRITE' )   ! calendar
100      CALL trc_wri      ( kt,      Kmm            )       ! output of passive tracers with iom I/O manager
101      CALL trc_sms      ( kt, Kbb, Kmm, Krhs      )       ! tracers: sinks and sources
102      CALL trc_trp      ( kt, Kbb, Kmm, Krhs, Kaa )       ! transport of passive tracers
103           !
104           ! Note passive tracers have been time-filtered in trc_trp but the time level
105           ! indices will not be swapped until after tra_atf/dyn_atf/ssh_atf in stp. Subsequent calls here
106           ! anticipate this update which will be: Nrhs= Nbb ; Nbb = Nnn ; Nnn = Naa ; Naa = Nrhs
107           ! and use the filtered levels explicitly.
108           !
109      IF( kt == nittrc000 ) THEN
110         CALL iom_close( numrtr )                         ! close input tracer restart file
111         IF(lwm) CALL FLUSH( numont )                     ! flush namelist output
112      ENDIF
113      IF( lrst_trc )            CALL trc_rst_wri  ( kt, Kmm, Kaa, Kbb  )       ! write tracer restart file
114      IF( lk_trdmxl_trc  )      CALL trd_mxl_trc  ( kt,      Kaa       )       ! trends: Mixed-layer
115      !
116      IF( ln_top_euler ) THEN 
117         ! For Euler timestepping for TOP we need to copy the "after" to the "now" fields
118         ! here then after the (leapfrog) swapping of the time-level indices in OCE/step.F90 we have
119         ! "before" fields = "now" fields.
120         tr(:,:,:,:,Kmm) = tr(:,:,:,:,Kaa)
121      ENDIF
122      !
123      IF (ll_trcstat) THEN
124         ztrai = 0._wp                                                   !  content of all tracers
125         DO jn = 1, jptra
126            ztrai = ztrai + glob_sum( 'trcstp', tr(:,:,:,jn,Kaa) * cvol(:,:,:)   )
127         END DO
128         IF( lwm ) WRITE(numstr,9300) kt,  ztrai / areatot
129      ENDIF
1309300  FORMAT(i10,D23.16)
131      !
132      IF( ln_timing )   CALL timing_stop('trc_stp')
133      !
134   END SUBROUTINE trc_stp
135
136
137   SUBROUTINE trc_stp_ctl
138      !!----------------------------------------------------------------------
139      !!                     ***  ROUTINE trc_stp_ctl  ***
140      !! ** Purpose :        Control  + ocean volume
141      !!----------------------------------------------------------------------
142      !
143      ! Define logical parameter ton control dirunal cycle in TOP
144      l_trcdm2dc = ln_dm2dc .OR. ( ln_cpl .AND. ncpl_qsr_freq /= 1 )
145      l_trcdm2dc = l_trcdm2dc  .AND. .NOT. l_offline
146      IF( l_trcdm2dc .AND. lwp )   CALL ctl_warn( 'Coupling with passive tracers and used of diurnal cycle.',   &
147         &                           'Computation of a daily mean shortwave for some biogeochemical models ' )
148      !
149   END SUBROUTINE trc_stp_ctl
150
151
152   SUBROUTINE trc_mean_qsr( kt )
153      !!----------------------------------------------------------------------
154      !!             ***  ROUTINE trc_mean_qsr  ***
155      !!
156      !! ** Purpose :  Compute daily mean qsr for biogeochemical model in case
157      !!               of diurnal cycle
158      !!
159      !! ** Method  : store in TOP the qsr every hour ( or every time-step if the latter
160      !!              is greater than 1 hour ) and then, compute the  mean with
161      !!              a moving average over 24 hours.
162      !!              In coupled mode, the sampling is done at every coupling frequency
163      !!----------------------------------------------------------------------
164      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
165      !
166      INTEGER  ::   jn   ! dummy loop indices
167      REAL(wp) ::   zkt, zrec     ! local scalars
168      CHARACTER(len=1) ::   cl1   ! 1 character
169      CHARACTER(len=2) ::   cl2   ! 2 characters
170      !!----------------------------------------------------------------------
171      !
172      IF( ln_timing )   CALL timing_start('trc_mean_qsr')
173      !
174      IF( kt == nittrc000 ) THEN
175         IF( ln_cpl )  THEN 
176            rdt_sampl = rday / ncpl_qsr_freq
177            nb_rec_per_day = ncpl_qsr_freq
178         ELSE 
179            rdt_sampl = MAX( 3600., rdttrc )
180            nb_rec_per_day = INT( rday / rdt_sampl )
181         ENDIF
182         !
183         IF(lwp) THEN
184            WRITE(numout,*) 
185            WRITE(numout,*) ' Sampling frequency dt = ', rdt_sampl, 's','   Number of sampling per day  nrec = ', nb_rec_per_day
186            WRITE(numout,*) 
187         ENDIF
188         !
189         ALLOCATE( qsr_arr(jpi,jpj,nb_rec_per_day ) )
190         !
191         !                                            !* Restart: read in restart file
192         IF( ln_rsttr .AND. nn_rsttr /= 0 .AND. iom_varid( numrtr, 'qsr_mean' , ldstop = .FALSE. ) > 0  &
193           &                              .AND. iom_varid( numrtr, 'qsr_arr_1', ldstop = .FALSE. ) > 0  &
194           &                              .AND. iom_varid( numrtr, 'ktdcy'    , ldstop = .FALSE. ) > 0  &
195           &                              .AND. iom_varid( numrtr, 'nrdcy'    , ldstop = .FALSE. ) > 0  ) THEN
196
197            CALL iom_get( numrtr, 'ktdcy', zkt ) 
198            rsecfst = INT( zkt ) * rdttrc
199            IF(lwp) WRITE(numout,*) 'trc_qsr_mean:   qsr_mean read in the restart file at time-step rsecfst =', rsecfst, ' s '
200            CALL iom_get( numrtr, jpdom_autoglo, 'qsr_mean', qsr_mean )   !  A mean of qsr
201            CALL iom_get( numrtr, 'nrdcy', zrec )   !  Number of record per days
202            IF( INT( zrec ) == nb_rec_per_day ) THEN
203               DO jn = 1, nb_rec_per_day 
204                  IF( jn <= 9 )  THEN
205                    WRITE(cl1,'(i1)') jn
206                    CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl1, qsr_arr(:,:,jn) )   !  A mean of qsr
207                  ELSE
208                    WRITE(cl2,'(i2.2)') jn
209                    CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl2, qsr_arr(:,:,jn) )   !  A mean of qsr
210                  ENDIF
211              END DO
212            ELSE
213               DO jn = 1, nb_rec_per_day
214                  qsr_arr(:,:,jn) = qsr_mean(:,:)
215               ENDDO
216            ENDIF
217         ELSE                                         !* no restart: set from nit000 values
218            IF(lwp) WRITE(numout,*) 'trc_qsr_mean:   qsr_mean set to nit000 values'
219            rsecfst  = kt * rdttrc
220            !
221            qsr_mean(:,:) = qsr(:,:)
222            DO jn = 1, nb_rec_per_day
223               qsr_arr(:,:,jn) = qsr_mean(:,:)
224            END DO
225         ENDIF
226         !
227      ENDIF
228      !
229      rseclast = kt * rdttrc
230      !
231      llnew   = ( rseclast - rsecfst ) .ge.  rdt_sampl    !   new shortwave to store
232      IF( llnew ) THEN
233          ktdcy = kt
234          IF( lwp .AND. kt < nittrc000 + 100 ) WRITE(numout,*) ' New shortwave to sample for TOP at time kt = ', ktdcy, &
235             &                      ' time = ', rseclast/3600.,'hours '
236          rsecfst = rseclast
237          DO jn = 1, nb_rec_per_day - 1
238             qsr_arr(:,:,jn) = qsr_arr(:,:,jn+1)
239          ENDDO
240          qsr_arr (:,:,nb_rec_per_day) = qsr(:,:)
241          qsr_mean(:,:                ) = SUM( qsr_arr(:,:,:), 3 ) / nb_rec_per_day
242      ENDIF
243      !
244      IF( lrst_trc ) THEN    !* Write the mean of qsr in restart file
245         IF(lwp) WRITE(numout,*)
246         IF(lwp) WRITE(numout,*) 'trc_mean_qsr : write qsr_mean in restart file  kt =', kt
247         IF(lwp) WRITE(numout,*) '~~~~~~~'
248         zkt  = REAL( ktdcy, wp )
249         zrec = REAL( nb_rec_per_day, wp )
250         CALL iom_rstput( kt, nitrst, numrtw, 'ktdcy', zkt  )
251         CALL iom_rstput( kt, nitrst, numrtw, 'nrdcy', zrec )
252          DO jn = 1, nb_rec_per_day 
253             IF( jn <= 9 )  THEN
254               WRITE(cl1,'(i1)') jn
255               CALL iom_rstput( kt, nitrst, numrtw, 'qsr_arr_'//cl1, qsr_arr(:,:,jn) )
256             ELSE
257               WRITE(cl2,'(i2.2)') jn
258               CALL iom_rstput( kt, nitrst, numrtw, 'qsr_arr_'//cl2, qsr_arr(:,:,jn) )
259             ENDIF
260         END DO
261         CALL iom_rstput( kt, nitrst, numrtw, 'qsr_mean', qsr_mean(:,:) )
262      ENDIF
263      !
264      IF( ln_timing )   CALL timing_stop('trc_mean_qsr')
265      !
266   END SUBROUTINE trc_mean_qsr
267
268#else
269   !!----------------------------------------------------------------------
270   !!   Default key                                     NO passive tracers
271   !!----------------------------------------------------------------------
272CONTAINS
273   SUBROUTINE trc_stp( kt )        ! Empty routine
274      WRITE(*,*) 'trc_stp: You should not have seen this print! error?', kt
275   END SUBROUTINE trc_stp
276#endif
277
278   !!======================================================================
279END MODULE trcstp
Note: See TracBrowser for help on using the repository browser.