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.
trcstp.F90 in NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/TOP – NEMO

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, 4 years 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
RevLine 
[1457]1MODULE trcstp
2   !!======================================================================
3   !!                       ***  MODULE trcstp  ***
4   !! Time-stepping    : time loop of opa for passive tracer
5   !!======================================================================
[2528]6   !! History :  1.0  !  2004-03  (C. Ethe)  Original
[12377]7   !!            4.1  !  2019-08  (A. Coward, D. Storkey) rewrite in preparation for new timestepping scheme
[2528]8   !!----------------------------------------------------------------------
[1457]9#if defined key_top
10   !!----------------------------------------------------------------------
[9019]11   !!   trc_stp       : passive tracer system time-stepping
[1457]12   !!----------------------------------------------------------------------
[9019]13   USE oce_trc        ! ocean dynamics and active tracers variables
[4306]14   USE sbc_oce
[2528]15   USE trc
[9019]16   USE trctrp         ! passive tracers transport
17   USE trcsms         ! passive tracers sources and sinks
[1457]18   USE trcwri
19   USE trcrst
[4990]20   USE trdtrc_oce
21   USE trdmxl_trc
[10425]22   USE sms_pisces,  ONLY : ln_check_mass
[9019]23   !
24   USE prtctl_trc     ! Print control for debbuging
25   USE iom            !
26   USE in_out_manager !
[1457]27
28   IMPLICIT NONE
29   PRIVATE
30
[2528]31   PUBLIC   trc_stp    ! called by step
[3294]32
[9019]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
[5385]38
[1457]39   !!----------------------------------------------------------------------
[10067]40   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
[7753]41   !! $Id$
[10068]42   !! Software governed by the CeCILL license (see ./LICENSE)
[1457]43   !!----------------------------------------------------------------------
44CONTAINS
45
[12377]46   SUBROUTINE trc_stp( kt, Kbb, Kmm, Krhs, Kaa )
[1457]47      !!-------------------------------------------------------------------
48      !!                     ***  ROUTINE trc_stp  ***
49      !!                     
[9019]50      !! ** Purpose :   Time loop of opa for passive tracer
[1457]51      !!
[9019]52      !! ** Method  :   Compute the passive tracers trends
53      !!                Update the passive tracers
[1457]54      !!-------------------------------------------------------------------
[12377]55      INTEGER, INTENT( in ) :: kt                  ! ocean time-step index
56      INTEGER, INTENT( in ) :: Kbb, Kmm, Krhs, Kaa ! time level indices
[9019]57      !
58      INTEGER ::   jk, jn   ! dummy loop indices
59      REAL(wp)::   ztrai    ! local scalar
[10570]60      LOGICAL ::   ll_trcstat ! local logical
[9019]61      CHARACTER (len=25) ::   charout   !
[2528]62      !!-------------------------------------------------------------------
[3294]63      !
[9124]64      IF( ln_timing )   CALL timing_start('trc_stp')
[3294]65      !
[12397]66      IF( l_1st_euler .OR. ln_top_euler ) THEN     ! at nittrc000
[7646]67         r2dttrc =  rdttrc           ! = rdttrc (use or restarting with Euler time stepping)
[12377]68      ELSEIF( kt <= nittrc000 + 1 ) THEN                                     ! at nittrc000 or nittrc000+1
69         r2dttrc = 2. * rdttrc       ! = 2 rdttrc (leapfrog)
[7646]70      ENDIF
71      !
[12377]72      ll_trcstat  = ( sn_cfctl%l_trcstat ) .AND. &
[10570]73     &              ( ( MOD( kt, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend ) )
[12136]74
75      IF( kt == nittrc000 )                      CALL trc_stp_ctl   ! control
[4990]76      IF( kt == nittrc000 .AND. lk_trdmxl_trc )  CALL trd_mxl_trc_init    ! trends: Mixed-layer
[3294]77      !
[6140]78      IF( .NOT.ln_linssh ) THEN                                           ! update ocean volume due to ssh temporal evolution
[3294]79         DO jk = 1, jpk
[12377]80            cvol(:,:,jk) = e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
[3294]81         END DO
[12377]82         IF ( ll_trcstat .OR. kt == nitrst .OR. ( ln_check_mass .AND. kt == nitend )   &
[10425]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(:,:,:) )
[3294]86      ENDIF
[5385]87      !
88      IF( l_trcdm2dc )   CALL trc_mean_qsr( kt )
[3294]89      !   
[12377]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)
[1457]94      ENDIF
[3294]95      !
[12377]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      !
[10570]123      IF (ll_trcstat) THEN
[10425]124         ztrai = 0._wp                                                   !  content of all tracers
125         DO jn = 1, jptra
[12377]126            ztrai = ztrai + glob_sum( 'trcstp', tr(:,:,:,jn,Kaa) * cvol(:,:,:)   )
[10425]127         END DO
128         IF( lwm ) WRITE(numstr,9300) kt,  ztrai / areatot
129      ENDIF
[6942]1309300  FORMAT(i10,D23.16)
[3294]131      !
[9124]132      IF( ln_timing )   CALL timing_stop('trc_stp')
[3294]133      !
[1457]134   END SUBROUTINE trc_stp
135
[12377]136
[12136]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
[9019]150
[12136]151
[5385]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      !!
[6942]159      !! ** Method  : store in TOP the qsr every hour ( or every time-step if the latter
[5385]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      !!----------------------------------------------------------------------
[9019]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      !
[9124]172      IF( ln_timing )   CALL timing_start('trc_mean_qsr')
173      !
[5385]174      IF( kt == nittrc000 ) THEN
[5407]175         IF( ln_cpl )  THEN 
[6968]176            rdt_sampl = rday / ncpl_qsr_freq
[6942]177            nb_rec_per_day = ncpl_qsr_freq
[5385]178         ELSE 
[6981]179            rdt_sampl = MAX( 3600., rdttrc )
[6968]180            nb_rec_per_day = INT( rday / rdt_sampl )
[5385]181         ENDIF
182         !
[9019]183         IF(lwp) THEN
[5385]184            WRITE(numout,*) 
[6942]185            WRITE(numout,*) ' Sampling frequency dt = ', rdt_sampl, 's','   Number of sampling per day  nrec = ', nb_rec_per_day
[5385]186            WRITE(numout,*) 
187         ENDIF
188         !
[6942]189         ALLOCATE( qsr_arr(jpi,jpj,nb_rec_per_day ) )
[5385]190         !
[6942]191         !                                            !* Restart: read in restart file
[7787]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 ) 
[7812]198            rsecfst = INT( zkt ) * rdttrc
[6968]199            IF(lwp) WRITE(numout,*) 'trc_qsr_mean:   qsr_mean read in the restart file at time-step rsecfst =', rsecfst, ' s '
[6942]200            CALL iom_get( numrtr, jpdom_autoglo, 'qsr_mean', qsr_mean )   !  A mean of qsr
[7787]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
[9019]211              END DO
[7787]212            ELSE
213               DO jn = 1, nb_rec_per_day
214                  qsr_arr(:,:,jn) = qsr_mean(:,:)
215               ENDDO
216            ENDIF
[6942]217         ELSE                                         !* no restart: set from nit000 values
218            IF(lwp) WRITE(numout,*) 'trc_qsr_mean:   qsr_mean set to nit000 values'
[6981]219            rsecfst  = kt * rdttrc
[6942]220            !
221            qsr_mean(:,:) = qsr(:,:)
222            DO jn = 1, nb_rec_per_day
223               qsr_arr(:,:,jn) = qsr_mean(:,:)
[9019]224            END DO
[6942]225         ENDIF
[5385]226         !
227      ENDIF
228      !
[6981]229      rseclast = kt * rdttrc
[6942]230      !
[6968]231      llnew   = ( rseclast - rsecfst ) .ge.  rdt_sampl    !   new shortwave to store
[6942]232      IF( llnew ) THEN
[7787]233          ktdcy = kt
234          IF( lwp .AND. kt < nittrc000 + 100 ) WRITE(numout,*) ' New shortwave to sample for TOP at time kt = ', ktdcy, &
[6968]235             &                      ' time = ', rseclast/3600.,'hours '
236          rsecfst = rseclast
[6942]237          DO jn = 1, nb_rec_per_day - 1
[5385]238             qsr_arr(:,:,jn) = qsr_arr(:,:,jn+1)
[6942]239          ENDDO
240          qsr_arr (:,:,nb_rec_per_day) = qsr(:,:)
241          qsr_mean(:,:                ) = SUM( qsr_arr(:,:,:), 3 ) / nb_rec_per_day
[5385]242      ENDIF
243      !
[6942]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,*) '~~~~~~~'
[7787]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 )
[6942]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
[9019]260         END DO
[6942]261         CALL iom_rstput( kt, nitrst, numrtw, 'qsr_mean', qsr_mean(:,:) )
262      ENDIF
263      !
[9124]264      IF( ln_timing )   CALL timing_stop('trc_mean_qsr')
265      !
[5385]266   END SUBROUTINE trc_mean_qsr
267
[1457]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.