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_new on Ticket #632 – Attachment – NEMO

Ticket #632: daymod.F90_new

File daymod.F90_new, 18.2 KB (added by cbricaud, 14 years ago)
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   !!----------------------------------------------------------------------     
14
15   !!----------------------------------------------------------------------
16   !!   day        : calendar
17   !! 
18   !!           -------------------------------
19   !!           ----------- WARNING -----------
20   !!
21   !!   we suppose that the time step is deviding the number of second of in a day
22   !!             ---> MOD( rday, rdttra(1) ) == 0
23   !!
24   !!           ----------- WARNING -----------
25   !!           -------------------------------
26   !! 
27   !!----------------------------------------------------------------------
28   USE dom_oce         ! ocean space and time domain
29   USE phycst          ! physical constants
30   USE in_out_manager  ! I/O manager
31   USE iom             !
32   USE ioipsl, ONLY :   ymds2ju   ! for calendar
33   USE prtctl          ! Print control
34   USE restart         !
35
36   IMPLICIT NONE
37   PRIVATE
38
39   PUBLIC   day        ! called by step.F90
40   PUBLIC   day_init   ! called by istate.F90
41
42   INTEGER ::   nsecd, nsecd05, ndt, ndt05
43   INTEGER , PUBLIC ::   njuloffset  !:
44   REAL(wp), PUBLIC ::   fjulstartyear !: firt day of the current year in julian days
45   REAL(wp), PUBLIC ::   fjulstartweek !: firt day of the current week in julian days
46   INTEGER , PUBLIC ::   nsec_week   !: current time step counted in second since 00h 1st day of the current week
47   CHARACTER(len=12),PARAMETER                ::   cfirstdayweek  = 'thursday'
48   CHARACTER(len=15),DIMENSION(7)             ::   cl_week
49
50   !!----------------------------------------------------------------------
51   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
52   !! $Id: daymod.F90 1730 2009-11-16 14:34:19Z smasson $
53   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
54   !!----------------------------------------------------------------------
55
56CONTAINS
57
58   SUBROUTINE day_init
59      !!----------------------------------------------------------------------
60      !!                   ***  ROUTINE day_init  ***
61      !!
62      !! ** Purpose :   Initialization of the calendar values to their values 1 time step before nit000
63      !!                because day will be called at the beginning of step
64      !!
65      !! ** Action  : - nyear        : current year
66      !!              - nmonth       : current month of the year nyear
67      !!              - nday         : current day of the month nmonth
68      !!              - nday_year    : current day of the year nyear
69      !!              - nsec_year    : current time step counted in second since 00h jan 1st of the current year
70      !!              - nsec_month   : current time step counted in second since 00h 1st day of the current month
71      !!              - nsec_day     : current time step counted in second since 00h of the current day
72      !!              - nsec1jan000  : second since Jan. 1st 00h of nit000 year and Jan. 1st 00h of the current year
73      !!              - nmonth_len, nyear_len, nmonth_half, nmonth_end through day_mth
74      !!----------------------------------------------------------------------
75      REAL(wp)::   rjulstartweek !: firt day of the current week in julian days
76
77      ! all calendar staff is based on the fact that MOD( rday, rdttra(1) ) == 0
78      IF( MOD( rday     , rdttra(1) ) /= 0. )   CALL ctl_stop( 'the time step must devide the number of second of in a day' )
79      IF( MOD( rday     , 2.        ) /= 0. )   CALL ctl_stop( 'the number of second of in a day must be an even number'    )
80      IF( MOD( rdttra(1), 2.        ) /= 0. )   CALL ctl_stop( 'the time step (in second) must be an even number'           )
81      nsecd   = NINT(rday           )
82      nsecd05 = NINT(0.5 * rday     )
83      ndt     = NINT(      rdttra(1))
84      ndt05   = NINT(0.5 * rdttra(1))
85
86      CALL day_rst( nit000, 'READ' )
87
88      ! set the calandar from ndastp (read in restart file and namelist)
89
90      nyear   =   ndastp / 10000
91      nmonth  = ( ndastp - (nyear * 10000) ) / 100
92      nday    =   ndastp - (nyear * 10000) - ( nmonth * 100 )
93
94      !compute njuloffset: it is the julian day that corresponds to the first day of the first week in the calendar
95      IF( nleapy == 0 )THEN ; cl_week=(/"wednesday","thursday ","friday   ","saturday ","sunday   ","monday   ","tuesday  "/)
96      ELSE                  ; cl_week=(/"thursday ","friday   ","saturday ","sunday   ","monday   ","tuesday  ","wednesday"/)
97      ENDIF
98      DO njuloffset=1,7
99         IF(  cl_week(njuloffset)==cfirstdayweek ) EXIT
100      ENDDO
101      IF( njuloffset .GT. 7 ) CALL ctl_stop( 'day_init: weekly file: wrong day for cfirstdayweek:  ',TRIM(cfirstdayweek) )
102      njuloffset=njuloffset-1
103      CALL ymds2ju( nyear, nmonth, nday, 0.0, fjulday )                                     ! we assume that we start run at 00:00
104      IF( ABS(fjulday - REAL(NINT(fjulday),wp)) < 0.1 / rday )  fjulday = REAL(NINT(fjulday),wp)   ! avoid truncation error
105      rjulstartweek=INT( (fjulday-njuloffset)/7 )*7 + njuloffset                            !compute fjulstart for nsec_week
106
107      fjulday = fjulday + 1.                             ! move back to the day at nit000 (and not at nit000 - 1)
108
109      nsec1jan000 = 0
110      CALL day_mth
111     
112      IF ( nday == 0 ) THEN     !   for ex if ndastp = ndate0 - 1
113         nmonth = nmonth - 1 
114         nday = nmonth_len(nmonth)
115      ENDIF
116      IF ( nmonth == 0 ) THEN   ! go at the end of previous year
117         nmonth = 12
118         nyear = nyear - 1
119         nsec1jan000 = nsec1jan000 - nsecd * nyear_len(0)
120         IF( nleapy == 1 )   CALL day_mth
121      ENDIF
122     
123      ! day since january 1st
124      nday_year = nday + SUM( nmonth_len(1:nmonth - 1) )
125     
126      ! number of seconds since the beginning of current year/month at the middle of the time-step
127      nsec_year  = nday_year * nsecd - ndt05   ! 1 time step before the middle of the first time step
128      nsec_month = nday      * nsecd - ndt05   ! because day will be called at the beginning of step
129      nsec_week  = REAL( fjulday   - rjulstartweek  ) *rday - 0.5 * rdttra(1)
130      nsec_day   =             nsecd - ndt05
131
132      ! control print
133      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 = ',   &
134           &                   nyear, '/', nmonth, '/', nday, '  nsec_day:', nsec_day
135
136      ! Up to now, calendar parameters are related to the end of previous run (nit000-1)
137      ! call day to set the calendar parameters at the begining of the current simulaton. needed by iom_init
138      CALL day( nit000 )
139
140     
141   END SUBROUTINE day_init
142
143
144   SUBROUTINE day_mth
145      !!----------------------------------------------------------------------
146      !!                   ***  ROUTINE day_init  ***
147      !!
148      !! ** Purpose :   calendar values related to the months
149      !!
150      !! ** Action  : - nmonth_len    : length in days of the months of the current year
151      !!              - nyear_len     : length in days of the previous/current year
152      !!              - nmonth_half   : second since the beginning of the year and the halft of the months
153      !!              - nmonth_end    : second since the beginning of the year and the end of the months
154      !!----------------------------------------------------------------------
155      INTEGER  ::   jm               ! dummy loop indice
156      !!----------------------------------------------------------------------
157
158      ! length of the month of the current year (from nleapy, read in namelist)
159      IF ( nleapy < 2 ) THEN
160         nmonth_len(:) = (/ 31, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31, 31 /)
161         nyear_len(:) = 365
162         IF ( nleapy == 1 ) THEN   ! we are using calandar with leap years
163            IF ( MOD(nyear-1, 4) == 0 .AND. ( MOD(nyear-1, 400) == 0 .OR. MOD(nyear-1, 100) /= 0 ) ) THEN
164               nyear_len(0) = 366
165            ENDIF
166            IF ( MOD(nyear, 4) == 0 .AND. ( MOD(nyear, 400) == 0 .OR. MOD(nyear, 100) /= 0 ) ) THEN
167               nmonth_len(2) = 29
168               nyear_len(1) = 366
169            ENDIF
170         ENDIF
171      ELSE
172         nmonth_len(:) = nleapy   ! all months with nleapy days per year
173         nyear_len(:) = 12 * nleapy
174      ENDIF
175
176      ! half month in second since the begining of the year:
177      ! time since Jan 1st   0     1     2    ...    11    12    13
178      !          ---------*--|--*--|--*--| ... |--*--|--*--|--*--|--------------------------------------
179      !                 <---> <---> <--->  ...  <---> <---> <--->       
180      ! month number      0     1     2    ...    11    12    13
181      !
182      ! nmonth_half(jm) = rday * REAL( 0.5 * nmonth_len(jm) + SUM(nmonth_len(1:jm-1)) )
183      nmonth_half(0) = - nsecd05 * nmonth_len(0)
184      DO jm = 1, 13
185         nmonth_half(jm) = nmonth_half(jm-1) + nsecd05 * ( nmonth_len(jm-1) + nmonth_len(jm) )
186      END DO
187
188      nmonth_end(0) = 0
189      DO jm = 1, 13
190         nmonth_end(jm) = nmonth_end(jm-1) + nsecd * nmonth_len(jm)
191      END DO
192      !           
193   END SUBROUTINE
194
195
196   SUBROUTINE day( kt )
197      !!----------------------------------------------------------------------
198      !!                      ***  ROUTINE day  ***
199      !!
200      !! ** Purpose :   Compute the date with a day iteration IF necessary.
201      !!
202      !! ** Method  : - ???
203      !!
204      !! ** Action  : - nyear     : current year
205      !!              - nmonth    : current month of the year nyear
206      !!              - nday      : current day of the month nmonth
207      !!              - nday_year : current day of the year nyear
208      !!              - ndastp    : = nyear*10000 + nmonth*100 + nday
209      !!              - adatrj    : date in days since the beginning of the run
210      !!              - nsec_year : current time of the year (in second since 00h, jan 1st)
211      !!----------------------------------------------------------------------     
212      INTEGER, INTENT(in) ::   kt        ! ocean time-step indices
213      !
214      CHARACTER (len=25) ::   charout
215      REAL(wp)           ::   zprec      ! fraction of day corresponding to 0.1 second
216      REAL(wp)           ::   rsec       !temp variable
217      !!----------------------------------------------------------------------
218      zprec = 0.1 / rday
219      !                                                 ! New time-step
220      nsec_year  = nsec_year  + ndt
221      nsec_month = nsec_month + ndt                 
222      nsec_week  = nsec_week  + ndt
223      nsec_day   = nsec_day   + ndt               
224      adatrj  = adatrj  + rdttra(1) / rday
225      fjulday = fjulday + rdttra(1) / rday
226      IF( ABS(fjulday - REAL(NINT(fjulday),wp)) < zprec )   fjulday = REAL(NINT(fjulday),wp)   ! avoid truncation error
227      IF( ABS(adatrj  - REAL(NINT(adatrj ),wp)) < zprec )   adatrj  = REAL(NINT(adatrj ),wp)   ! avoid truncation error
228     
229      IF( nsec_day > nsecd ) THEN                        ! NEW day
230         !
231         nday      = nday + 1
232         nday_year = nday_year + 1
233         nsec_day  = ndt05
234         !
235         IF(nsec_week > 86400*7)nsec_week=0.5 * ndt
236         IF( nday == nmonth_len(nmonth) + 1 ) THEN      ! NEW month
237            nday   = 1
238            nmonth = nmonth + 1
239            nsec_month = ndt05
240            IF( nmonth == 13 ) THEN                     ! NEW year
241               nyear     = nyear + 1
242               nmonth    = 1
243               nday_year = 1
244               nsec_year = ndt05
245               nsec1jan000 = nsec1jan000 + nsecd * nyear_len(1)
246               IF( nleapy == 1 )   CALL day_mth
247            ENDIF
248         ENDIF
249         !
250         ndastp = nyear * 10000 + nmonth * 100 + nday   ! NEW date
251         !
252         !compute first day of the year in julian days
253         CALL ymds2ju( nyear, 01, 01, 0.0, fjulstartyear )
254         !compute first day of the week  in julian days
255         fjulstartweek=INT( (fjulday-njuloffset)/7 )*7 + njuloffset !compute fjulstartweek
256
257         IF(lwp) WRITE(numout,'(a,i8,a,i4.4,a,i2.2,a,i2.2,a,i3.3)') '======>> time-step =', kt,   &
258              &   '      New day, DATE Y/M/D = ', nyear, '/', nmonth, '/', nday, '      nday_year = ', nday_year
259         IF(lwp) WRITE(numout,'(a,i8,a,i7,a,i7,a,i5)') '         nsec_year = ', nsec_year,   &
260              &   '   nsec_month = ', nsec_month, '   nsec_week = ', nsec_week, '   nsec_day = ', nsec_day
261      ENDIF
262     
263      IF(ln_ctl) THEN
264         WRITE(charout,FMT="('kt =', I4,'  d/m/y =',I2,I2,I4)") kt, nday, nmonth, nyear
265         CALL prt_ctl_info(charout)
266      ENDIF
267
268      IF( lrst_oce )   CALL day_rst( kt, 'WRITE' )
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 - 1 (as we start the step loop with a call to day)
332               ndastp = ndate0 - 1     ! ndate0 read in the namelist in dom_nam, we assume that we start run at 00:00
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 - 1 (as we start the step loop with a call to day)
338            ndastp = ndate0 - 1        ! ndate0 read in the namelist in dom_nam, we assume that we start run at 00:00
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