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.
daymod.F90 on Ticket #1037 – Attachment – NEMO

Ticket #1037: daymod.F90

File daymod.F90, 17.5 KB (added by poddo, 11 years ago)

new version of daymod, more general and flexible

Line 
1MODULE daymod
2   !!======================================================================
3   !!                       ***  MODULE  daymod  ***
4   !! Ocean        :  calendar
5   !!=====================================================================
6   !! History :  OPA  ! 1994-09  (M. Pontaud M. Imbard)  Original code
7   !!                 ! 1997-03  (O. Marti)
8   !!                 ! 1997-05  (G. Madec)
9   !!                 ! 1997-08  (M. Imbard)
10   !!   NEMO     1.0  ! 2003-09  (G. Madec)  F90 + nyear, nmonth, nday
11   !!                 ! 2004-01  (A.M. Treguier) new calculation based on adatrj
12   !!                 ! 2006-08  (G. Madec)  surface module major update
13   !!                 ! 2012-12  (P. Oddo) generalizzation of the algorithm
14   !!                       it is now possible to start  at any time of the day
15   !!----------------------------------------------------------------------     
16
17   !!----------------------------------------------------------------------
18   !!   day        : calendar
19   !! 
20   !!           -------------------------------
21   !!           ----------- WARNING -----------
22   !!
23   !!   we suppose that the time step is deviding the number of second of in a day
24   !!             ---> MOD( rday, rdttra(1) ) == 0
25   !!
26   !!           ----------- WARNING -----------
27   !!           -------------------------------
28   !! 
29   !!----------------------------------------------------------------------
30   USE dom_oce         ! ocean space and time domain
31   USE phycst          ! physical constants
32   USE in_out_manager  ! I/O manager
33   USE iom             !
34   USE ioipsl, ONLY :   ymds2ju,ju2ymds   ! for calendar
35   USE prtctl          ! Print control
36   USE restart         !
37   USE trc_oce, ONLY : lk_offline ! offline flag
38   USE timing          ! Timing
39
40   IMPLICIT NONE
41   PRIVATE
42
43   PUBLIC   day        ! called by step.F90
44   PUBLIC   day_init   ! called by istate.F90
45
46   INTEGER ::   nsecd, nsecd05, ndt, ndt05
47
48   !!----------------------------------------------------------------------
49   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
50   !! $Id$
51   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
52   !!----------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE day_init
56      !!----------------------------------------------------------------------
57      !!                   ***  ROUTINE day_init  ***
58      !!
59      !! ** Purpose :   Initialization of the calendar values to their values 1 time step before nit000
60      !!                because day will be called at the beginning of step
61      !!
62      !! ** Action  : - nyear        : current year
63      !!              - nmonth       : current month of the year nyear
64      !!              - nday         : current day of the month nmonth
65      !!              - nday_year    : current day of the year nyear
66      !!              - nsec_year    : current time step counted in second since 00h jan 1st of the current year
67      !!              - nsec_month   : current time step counted in second since 00h 1st day of the current month
68      !!              - nsec_day     : current time step counted in second since 00h of the current day
69      !!              - nsec1jan000  : second since Jan. 1st 00h of nit000 year and Jan. 1st 00h of the current year
70      !!              - nmonth_len, nyear_len, nmonth_half, nmonth_end through day_mth
71      !!----------------------------------------------------------------------
72      INTEGER  ::   inbday, idweek
73      REAL(wp) ::   zjul,fracday
74      !!----------------------------------------------------------------------
75      !
76      ! all calendar staff is based on the fact that MOD( rday, rdttra(1) ) == 0
77      IF( MOD( rday     , rdttra(1) ) /= 0. )   CALL ctl_stop( 'the time step must devide the number of second of in a day' )
78      IF( MOD( rday     , 2.        ) /= 0. )   CALL ctl_stop( 'the number of second of in a day must be an even number'    )
79      IF( MOD( rdttra(1), 2.        ) /= 0. )   CALL ctl_stop( 'the time step (in second) must be an even number'           )
80      nsecd   = NINT(rday           )
81      nsecd05 = NINT(0.5 * rday     )
82      ndt     = NINT(      rdttra(1))
83      ndt05   = NINT(0.5 * rdttra(1))
84
85      IF( .NOT. lk_offline ) CALL day_rst( nit000, 'READ' ) 
86
87      ! set the calandar from ndastp (read in restart file and namelist)
88
89      nyear   =   ndastp / 10000
90      nmonth  = ( ndastp - (nyear * 10000) ) / 100
91      nday    =   ndastp - (nyear * 10000) - ( nmonth * 100 ) 
92 
93      ! round adatrj to the 6 digits potential problem with very small time-step
94      adatrj= NINT(adatrj*1000000.0_wp)/1000000.0_wp
95      nsec_day= INT((adatrj-INT(adatrj))*rday)
96      ! compute the julian
97      CALL ymds2ju( nyear, nmonth, nday, REAL(nsec_day), fjulday ) 
98
99      ! Compute date (year, month, day, second) 1/2 time-step before the integration start
100      CALL ju2ymds(fjulday-(ndt05/rday),nyear,nmonth,nday,fracday)
101      nsec_day=INT(fracday)
102
103      ! Re-Compute the julian day at  1/2 time-step before the integration start
104      CALL ymds2ju( nyear, nmonth, nday, REAL(nsec_day), fjulday ) 
105
106      ! Maybe the following statement can be removed
107      IF( ABS(fjulday - REAL(NINT(fjulday),wp)) < 0.1 / rday )   fjulday = REAL(NINT(fjulday),wp)   ! avoid truncation error
108
109      nsec1jan000 = 0
110      CALL day_mth
111
112      IF (nmonth == 12) THEN
113         nsec1jan000 = nsec1jan000 - nsecd * nyear_len(0)
114      ENDIF 
115
116      ! day since january 1st
117      nday_year = nday + SUM( nmonth_len(1:nmonth - 1) )
118
119      !compute number of days between last monday and today     
120      CALL ymds2ju( 1900, 01, 01, 0.0, zjul )  ! compute julian day value of 01.01.1900 (our reference that was a Monday)
121      inbday = NINT(INT(fjulday) - zjul)            ! compute nb day between  01.01.1900 and current day 
122      idweek = MOD(inbday, 7)                  ! compute nb day between last monday and current day 
123
124      ! number of seconds since the beginning of current year/month/week/day at the middle of the time-step
125      nsec_week  = idweek    * nsecd + nsec_day 
126      nsec_month = (nday-1)  * nsecd + nsec_day 
127      nsec_year  = (nday_year-1) * nsecd + nsec_day 
128
129      ! control print
130      IF(lwp) WRITE(numout,'(a,i6,a,i2,a,i2,a,i6)')' ==============>> 1/2 time step before the start of the run DATE Y/M/D = ',   &
131           &                   nyear, '/', nmonth, '/', nday, '  nsec_day:', nsec_day, '  nsec_week:', nsec_week
132
133      ! Up to now, calendar parameters are related to the end of previous run (nit000-1)
134      ! call day to set the calendar parameters at the begining of the current simulaton. needed by iom_init
135      CALL day( nit000 )
136      !
137   END SUBROUTINE day_init
138
139
140   SUBROUTINE day_mth
141      !!----------------------------------------------------------------------
142      !!                   ***  ROUTINE day_init  ***
143      !!
144      !! ** Purpose :   calendar values related to the months
145      !!
146      !! ** Action  : - nmonth_len    : length in days of the months of the current year
147      !!              - nyear_len     : length in days of the previous/current year
148      !!              - nmonth_half   : second since the beginning of the year and the halft of the months
149      !!              - nmonth_end    : second since the beginning of the year and the end of the months
150      !!----------------------------------------------------------------------
151      INTEGER  ::   jm               ! dummy loop indice
152      !!----------------------------------------------------------------------
153
154      ! length of the month of the current year (from nleapy, read in namelist)
155      IF ( nleapy < 2 ) THEN
156         nmonth_len(:) = (/ 31, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31, 31 /)
157         nyear_len(:) = 365
158         IF ( nleapy == 1 ) THEN   ! we are using calandar with leap years
159            IF ( MOD(nyear-1, 4) == 0 .AND. ( MOD(nyear-1, 400) == 0 .OR. MOD(nyear-1, 100) /= 0 ) ) THEN
160               nyear_len(0) = 366
161            ENDIF
162            IF ( MOD(nyear, 4) == 0 .AND. ( MOD(nyear, 400) == 0 .OR. MOD(nyear, 100) /= 0 ) ) THEN
163               nmonth_len(2) = 29
164               nyear_len(1) = 366
165            ENDIF
166         ENDIF
167      ELSE
168         nmonth_len(:) = nleapy   ! all months with nleapy days per year
169         nyear_len(:) = 12 * nleapy
170      ENDIF
171
172      ! half month in second since the begining of the year:
173      ! time since Jan 1st   0     1     2    ...    11    12    13
174      !          ---------*--|--*--|--*--| ... |--*--|--*--|--*--|--------------------------------------
175      !                 <---> <---> <--->  ...  <---> <---> <--->       
176      ! month number      0     1     2    ...    11    12    13
177      !
178      ! nmonth_half(jm) = rday * REAL( 0.5 * nmonth_len(jm) + SUM(nmonth_len(1:jm-1)) )
179      nmonth_half(0) = - nsecd05 * nmonth_len(0)
180      DO jm = 1, 13
181         nmonth_half(jm) = nmonth_half(jm-1) + nsecd05 * ( nmonth_len(jm-1) + nmonth_len(jm) )
182      END DO
183
184      nmonth_end(0) = 0
185      DO jm = 1, 13
186         nmonth_end(jm) = nmonth_end(jm-1) + nsecd * nmonth_len(jm)
187      END DO
188      !           
189   END SUBROUTINE
190
191
192   SUBROUTINE day( kt )
193      !!----------------------------------------------------------------------
194      !!                      ***  ROUTINE day  ***
195      !!
196      !! ** Purpose :   Compute the date with a day iteration IF necessary.
197      !!
198      !! ** Method  : - ???
199      !!
200      !! ** Action  : - nyear     : current year
201      !!              - nmonth    : current month of the year nyear
202      !!              - nday      : current day of the month nmonth
203      !!              - nday_year : current day of the year nyear
204      !!              - ndastp    : = nyear*10000 + nmonth*100 + nday
205      !!              - adatrj    : date in days since the beginning of the run
206      !!              - nsec_year : current time of the year (in second since 00h, jan 1st)
207      !!----------------------------------------------------------------------     
208      INTEGER, INTENT(in) ::   kt        ! ocean time-step indices
209      !
210      CHARACTER (len=25) ::   charout
211      REAL(wp)           ::   zprec      ! fraction of day corresponding to 0.1 second
212      !!----------------------------------------------------------------------
213      !
214      IF( nn_timing == 1 )  CALL timing_start('day')
215      !
216      zprec = 0.1 / rday
217      !                                                 ! New time-step
218      nsec_year  = nsec_year  + ndt 
219      nsec_month = nsec_month + ndt                 
220      nsec_week  = nsec_week  + ndt
221      nsec_day   = nsec_day   + ndt               
222      adatrj  = adatrj  + rdttra(1) / rday
223      fjulday = fjulday + rdttra(1) / rday
224      IF( ABS(fjulday - REAL(NINT(fjulday),wp)) < zprec )   fjulday = REAL(NINT(fjulday),wp)   ! avoid truncation error
225      IF( ABS(adatrj  - REAL(NINT(adatrj ),wp)) < zprec )   adatrj  = REAL(NINT(adatrj ),wp)   ! avoid truncation error
226     
227      IF( nsec_day > nsecd ) THEN                       ! New day
228         !
229         nday      = nday + 1
230         nday_year = nday_year + 1
231         nsec_day  = ndt05
232         !
233         IF( nday == nmonth_len(nmonth) + 1 ) THEN      ! New month
234            nday   = 1
235            nmonth = nmonth + 1
236            nsec_month = ndt05
237            IF( nmonth == 13 ) THEN                     ! New year
238               nyear     = nyear + 1
239               nmonth    = 1
240               nday_year = 1
241               nsec_year = ndt05
242               nsec1jan000 = nsec1jan000 + nsecd * nyear_len(1)
243               IF( nleapy == 1 )   CALL day_mth
244            ENDIF
245         ENDIF
246         !
247         ndastp = nyear * 10000 + nmonth * 100 + nday   ! New date
248         !
249         !compute first day of the year in julian days
250         CALL ymds2ju( nyear, 01, 01, 0.0, fjulstartyear )
251         !
252         IF(lwp) WRITE(numout,'(a,i8,a,i4.4,a,i2.2,a,i2.2,a,i3.3)') '======>> time-step =', kt,   &
253              &   '      New day, DATE Y/M/D = ', nyear, '/', nmonth, '/', nday, '      nday_year = ', nday_year
254         IF(lwp) WRITE(numout,'(a,i8,a,i7,a,i5)') '         nsec_year = ', nsec_year,   &
255              &   '   nsec_month = ', nsec_month, '   nsec_day = ', nsec_day, '   nsec_week = ', nsec_week
256      ENDIF
257
258      IF( nsec_week > 7*nsecd )   nsec_week = ndt05     ! New week
259     
260      IF(ln_ctl) THEN
261         WRITE(charout,FMT="('kt =', I4,'  d/m/y =',I2,I2,I4)") kt, nday, nmonth, nyear
262         CALL prt_ctl_info(charout)
263      ENDIF
264
265      IF( .NOT. lk_offline ) CALL rst_opn( kt )               ! Open the restart file if needed and control lrst_oce
266      IF( lrst_oce         ) CALL day_rst( kt, 'WRITE' )      ! write day restart information
267      !
268      IF( nn_timing == 1 )  CALL timing_stop('day')
269      !
270   END SUBROUTINE day
271
272
273   SUBROUTINE day_rst( kt, cdrw )
274      !!---------------------------------------------------------------------
275      !!                   ***  ROUTINE ts_rst  ***
276      !!
277      !!  ** Purpose : Read or write calendar in restart file:
278      !!
279      !!  WRITE(READ) mode:
280      !!       kt        : number of time step since the begining of the experiment at the
281      !!                   end of the current(previous) run
282      !!       adatrj(0) : number of elapsed days since the begining of the experiment at the
283      !!                   end of the current(previous) run (REAL -> keep fractions of day)
284      !!       ndastp    : date at the end of the current(previous) run (coded as yyyymmdd integer)
285      !!
286      !!   According to namelist parameter nrstdt,
287      !!       nrstdt = 0  no control on the date (nit000 is  arbitrary).
288      !!       nrstdt = 1  we verify that nit000 is equal to the last
289      !!                   time step of previous run + 1.
290      !!       In both those options, the  exact duration of the experiment
291      !!       since the beginning (cumulated duration of all previous restart runs)
292      !!       is not stored in the restart and is assumed to be (nit000-1)*rdt.
293      !!       This is valid is the time step has remained constant.
294      !!
295      !!       nrstdt = 2  the duration of the experiment in days (adatrj)
296      !!                    has been stored in the restart file.
297      !!----------------------------------------------------------------------
298      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
299      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
300      !
301      REAL(wp) ::   zkt, zndastp
302      !!----------------------------------------------------------------------
303     
304      IF( TRIM(cdrw) == 'READ' ) THEN
305
306         IF( iom_varid( numror, 'kt', ldstop = .FALSE. ) > 0 ) THEN
307            ! Get Calendar informations
308            CALL iom_get( numror, 'kt', zkt )   ! last time-step of previous run
309            IF(lwp) THEN
310               WRITE(numout,*) ' *** Info read in restart : '
311               WRITE(numout,*) '   previous time-step                               : ', NINT( zkt )
312               WRITE(numout,*) ' *** restart option'
313               SELECT CASE ( nrstdt )
314               CASE ( 0 )   ;   WRITE(numout,*) ' nrstdt = 0 : no control of nit000'
315               CASE ( 1 )   ;   WRITE(numout,*) ' nrstdt = 1 : no control the date at nit000 (use ndate0 read in the namelist)'
316               CASE ( 2 )   ;   WRITE(numout,*) ' nrstdt = 2 : calendar parameters read in restart'
317               END SELECT
318               WRITE(numout,*)
319            ENDIF
320            ! Control of date
321            IF( nit000 - NINT( zkt ) /= 1 .AND. nrstdt /= 0 )                                         & 
322                 &   CALL ctl_stop( ' ===>>>> : problem with nit000 for the restart',                 & 
323                 &                  ' verify the restart file or rerun with nrstdt = 0 (namelist)' )
324            ! define ndastp and adatrj
325            IF ( nrstdt == 2 ) THEN 
326               ! read the parameters correspondting to nit000 - 1 (last time step of previous run)
327               CALL iom_get( numror, 'ndastp', zndastp )
328               ndastp = NINT( zndastp )
329               CALL iom_get( numror, 'adatrj', adatrj  )
330            ELSE 
331               ! parameters correspondting to nit000
332               ndastp = ndate0 
333               adatrj = ( REAL( nit000 -1, wp ) * rdttra(1) ) / rday 
334               ! note this is wrong if time step has changed during run
335            ENDIF
336         ELSE
337            ! parameters correspondting to nit000
338            ndastp = ndate0 
339            adatrj = ( REAL( nit000 -1, wp ) * rdttra(1) ) / rday 
340         ENDIF
341         IF( ABS(adatrj  - REAL(NINT(adatrj),wp)) < 0.1 / rday )   adatrj = REAL(NINT(adatrj),wp)   ! avoid truncation error
342         !
343         IF(lwp) THEN
344            WRITE(numout,*) ' *** Info used values : '
345            WRITE(numout,*) '   date ndastp                                      : ', ndastp
346            WRITE(numout,*) '   number of elapsed days since the begining of run : ', adatrj
347            WRITE(numout,*)
348         ENDIF
349         !
350      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
351         !
352         IF( kt == nitrst ) THEN
353            IF(lwp) WRITE(numout,*)
354            IF(lwp) WRITE(numout,*) 'rst_write : write oce restart file  kt =', kt
355            IF(lwp) WRITE(numout,*) '~~~~~~~'         
356         ENDIF
357         ! calendar control
358         CALL iom_rstput( kt, nitrst, numrow, 'kt'     , REAL( kt    , wp) )   ! time-step
359         CALL iom_rstput( kt, nitrst, numrow, 'ndastp' , REAL( ndastp, wp) )   ! date
360         CALL iom_rstput( kt, nitrst, numrow, 'adatrj' , adatrj            )   ! number of elapsed days since
361         !                                                                     ! the begining of the run [s]
362      ENDIF
363      !
364   END SUBROUTINE day_rst
365
366   !!======================================================================
367END MODULE daymod