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.
sbcdcy.F90 in branches/dev_r2174_DCY/NEMO/OPA_SRC/SBC – NEMO

source: branches/dev_r2174_DCY/NEMO/OPA_SRC/SBC/sbcdcy.F90 @ 2187

Last change on this file since 2187 was 2187, checked in by smasson, 14 years ago

add diurna cycle in dev_r2174_DCY, see ticket:730
File size: 9.5 KB
Line 
1MODULE sbcdcy
2   !!======================================================================
3   !!                    ***  MODULE  sbcdcy  ***
4   !! Ocean forcing:  compute the diurnal cycle
5   !!======================================================================
6   !! History : 8.2  !  2005-02  (D. Bernie)  Original code
7   !!           9.0  !  2006-02  (S. Masson, G. Madec)  adaptation to OPA9
8   !!           3.1  !  2009-07  (J.M. Molines)  adaptation to nemo3.1
9   !!----------------------------------------------------------------------
10
11   !!----------------------------------------------------------------------
12   !!  sbc_dcy : compute solar flux at kt from daily mean, taking
13   !!            diurnal cycle into account
14   !!----------------------------------------------------------------------
15   USE oce              ! ocean dynamics and tracers
16   USE phycst           ! ocean physics
17   USE dom_oce          ! ocean space and time domain
18   USE in_out_manager   ! I/O manager
19
20   IMPLICIT NONE
21   PRIVATE
22   INTEGER                      ::   idayqsr                                            ! day when parameters were computed
23   REAL(wp), DIMENSION(jpi,jpj) ::   zaaa, zbbb, zccc, zab, ztmd, zdawn, zdusk, zscal   ! parameters to compute the diurnal cycle
24   REAL(wp), DIMENSION(jpi,jpj) ::   qsr_daily                                          ! to hold daily mean QSR
25 
26   PUBLIC sbc_dcy       ! routine called by sbc
27
28   !!----------------------------------------------------------------------
29   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)
30   !! $Id$
31   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
32   !!----------------------------------------------------------------------
33
34CONTAINS
35
36      SUBROUTINE sbc_dcy( kt, pqsr )
37      !!----------------------------------------------------------------------
38      !!                  ***  ROUTINE sbc_dcy  ***
39      !!
40      !! ** Purpose : introduce a diurnal cycle of qsr from daily values
41      !!
42      !! ** Method  : see Appendix A of
43      !!              Bernie, DJ, Guilyardi, E, Madec, G, Slingo, JM and Woolnough, SJ
44      !!              Impact of resolving the diurnal cycle in an ocean--atmosphere GCM. Part 1: a diurnally forced OGCM
45      !!              Climate Dynamics 29:6, 575-590 (2007)
46      !!
47      !! ** Action  : redistribute daily QSR on each time step following the diurnal cycle
48      !!----------------------------------------------------------------------
49      INTEGER,                      INTENT( in    ) ::   kt     ! ocean time-step index
50      REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) ::   pqsr   ! QSR flux with diurnal cycle
51      !!
52      INTEGER  ::   ji, jj                                      ! dummy loop indices
53      REAL(wp) ::   fintegral, pt1, pt2, paaa, pbbb, pccc       !
54      REAL(wp) ::   ztwopi, zinvtwopi, zconvrad 
55      REAL(wp) ::   zlo, zup, zlousd, zupusd
56      REAL(wp) ::   zdsws, zdecrad, ztx
57      REAL(wp) ::   ztmp, ztmp1, ztmp2, ztest
58      !!---------------------------------------------------------------------
59
60      !----------------------------------------------------------------------
61      ! statement functions
62
63      fintegral(pt1, pt2, paaa, pbbb, pccc) =                           &
64         &   paaa * pt2 + zinvtwopi * pbbb * SIN(pccc + ztwopi * pt2)   &
65         & - paaa * pt1 - zinvtwopi * pbbb * SIN(pccc + ztwopi * pt1)
66      !----------------------------------------------------------------------
67
68      ! Initialization
69      ! --------------
70      ztwopi    = 2. * rpi
71      zinvtwopi = 1. / ztwopi
72      zconvrad  = ztwopi / 360.
73
74      ! When are we during the day (from 0 to 1)
75      zlo = MOD( rdt / rday * REAL( kt-nit000, wp ), 1.)
76      zup = zlo + rdt / rday
77
78      !                                         
79      IF (kt == nit000) THEN                   
80         !
81         IF(lwp) THEN
82            WRITE(numout,*)
83            WRITE(numout,*) 'sbc_dcy : introduce diurnal cycle from daily mean qsr'
84            WRITE(numout,*) '~~~~~~~'
85            WRITE(numout,*)
86         ENDIF
87         idayqsr = 0
88         ! Compute C needed to compute the time integral of the diurnal cycle
89         zccc(:,:) = zconvrad * glamt(:,:) - rpi
90         ! time of midday
91         ztmd(:,:) = 0.5 - glamt(:,:) / 360.
92         ztmd(:,:) = MOD((ztmd(:,:) + 1.), 1.)
93      ENDIF
94
95      ! If this is a new day, we have to update the dawn, dusk and scaling function 
96      !----------------------
97   
98      !     2.1 dawn and dusk 
99
100      ! nday is the number of days since the beginning of the current month
101      IF( idayqsr /= nday ) THEN 
102         ! save the day of the year and the daily mean of qsr
103         idayqsr = nday 
104         ! number of days since the previous winter solstice (supposed to be always 21 December)         
105         zdsws = 11 + nday_year
106         ! declination of the earths orbit
107         zdecrad = (-23.5 * zconvrad) * COS( zdsws * ztwopi / REAL(nyear_len(1),wp) )
108         ! save the daily QSR for nest hours of the day
109         qsr_daily(:,:) = pqsr(:,:)
110         ! Compute A and B needed to compute the time integral of the diurnal cycle
111       
112         DO jj = 1, jpj
113            DO ji = 1, jpi
114               ztmp = zconvrad * gphit(ji,jj)
115               zaaa(ji,jj) = SIN(ztmp) * SIN(zdecrad)
116               zbbb(ji,jj) = COS(ztmp) * COS(zdecrad)
117            END DO 
118         END DO 
119
120         ! Compute the time of dawn and dusk
121
122         ! zab to test if the day time is equal to 0, less than 24h of full day       
123         zab(:,:) = -zaaa(:,:) / zbbb(:,:)
124         DO jj = 1, jpj
125            DO ji = 1, jpi
126               IF ( ABS(zab(ji,jj)) < 1 ) THEN
127         ! day duration is less than 24h
128         ! When is it night?
129                  ztx = zinvtwopi * (ACOS(zab(ji,jj)) - zccc(ji,jj))
130                  ztest = -zbbb(ji,jj) * SIN( zccc(ji,jj) + ztwopi * ztx )
131         ! is it dawn or dusk?
132                  IF ( ztest > 0 ) THEN
133                     zdawn(ji,jj) = ztx
134                     zdusk(ji,jj) = ztmd(ji,jj) + ( ztmd(ji,jj) - zdawn(ji,jj) )
135                  ELSE
136                     zdusk(ji,jj) = ztx
137                     zdawn(ji,jj) = ztmd(ji,jj) - ( zdusk(ji,jj) - ztmd(ji,jj) )
138                  ENDIF
139               ELSE
140                  zdawn(ji,jj) = ztmd(ji,jj) + 0.5
141                  zdusk(ji,jj) = zdawn(ji,jj)
142               ENDIF
143             END DO 
144         END DO 
145         zdawn(:,:) = MOD((zdawn(:,:) + 1.), 1.)
146         zdusk(:,:) = MOD((zdusk(:,:) + 1.), 1.)
147
148
149         !     2.2 Compute the scalling function:
150         !         S* = the inverse of the time integral of the diurnal cycle from dawm to dusk
151
152         DO jj = 1, jpj
153            DO ji = 1, jpi
154               IF ( ABS(zab(ji,jj)) < 1 ) THEN
155         ! day duration is less than 24h
156                  IF ( zdawn(ji,jj) < zdusk(ji,jj) ) THEN
157         ! day time in one part
158                     zscal(ji,jj) = fintegral(zdawn(ji,jj), zdusk(ji,jj), zaaa(ji,jj), zbbb(ji,jj), zccc(ji,jj)) 
159                     zscal(ji,jj) = 1. / zscal(ji,jj)
160                  ELSE
161         ! day time in two parts
162                     zscal(ji,jj) = fintegral(0., zdusk(ji,jj), zaaa(ji,jj), zbbb(ji,jj), zccc(ji,jj))   &
163                        &         + fintegral(zdawn(ji,jj), 1., zaaa(ji,jj), zbbb(ji,jj), zccc(ji,jj)) 
164                     zscal(ji,jj) = 1. / zscal(ji,jj)
165                  ENDIF
166               ELSE
167                  IF ( zaaa(ji,jj) > zbbb(ji,jj) ) THEN
168         ! 24h day
169                     zscal(ji,jj) = fintegral(0., 1., zaaa(ji,jj), zbbb(ji,jj), zccc(ji,jj)) 
170                     zscal(ji,jj) = 1. / zscal(ji,jj)
171                  ELSE
172         ! No day
173                     zscal(ji,jj) = 0.
174                  ENDIF
175               ENDIF
176            END DO 
177         END DO 
178         !
179         ztmp = rday / rdt
180         zscal(:,:) = zscal(:,:) * ztmp
181
182      ENDIF 
183
184         !     3. compute qsr with the diurnal cycle
185         !     -----------------------
186
187      DO jj = 1, jpj
188         DO ji = 1, jpi
189            IF ( ABS(zab(ji,jj)) < 1 ) THEN
190         ! day duration is less than 24h
191                  IF ( zdawn(ji,jj) < zdusk(ji,jj) ) THEN
192         ! day time in one part
193                     zlousd = MAX(zlo, zdawn(ji,jj))
194                     zlousd = MIN(zlousd, zup)
195                     zupusd = MIN(zup, zdusk(ji,jj))
196                     zupusd = MAX(zupusd, zlo)
197                     ztmp = fintegral(zlousd, zupusd, zaaa(ji,jj), zbbb(ji,jj), zccc(ji,jj)) 
198                     pqsr(ji,jj) = qsr_daily(ji,jj) * ztmp * zscal(ji,jj)
199                  ELSE
200         ! day time in two parts
201                     zlousd = MIN(zlo, zdusk(ji,jj))
202                     zupusd = MIN(zup, zdusk(ji,jj))
203                     ztmp1 = fintegral(zlousd, zupusd, zaaa(ji,jj), zbbb(ji,jj), zccc(ji,jj)) 
204                     zlousd = MAX(zlo, zdawn(ji,jj))
205                     zupusd = MAX(zup, zdawn(ji,jj))
206                     ztmp2 = fintegral(zlousd, zupusd, zaaa(ji,jj), zbbb(ji,jj), zccc(ji,jj)) 
207                     ztmp = ztmp1 + ztmp2
208                     pqsr(ji,jj) = qsr_daily(ji,jj) * ztmp * zscal(ji,jj)
209                  ENDIF
210            ELSE
211                  IF ( zaaa(ji,jj) > zbbb(ji,jj) ) THEN
212         ! 24h day
213                     ztmp = fintegral(zlo, zup, zaaa(ji,jj), zbbb(ji,jj), zccc(ji,jj)) 
214                     pqsr(ji,jj) = qsr_daily(ji,jj) * ztmp * zscal(ji,jj)
215                  ELSE
216         ! No day
217                     pqsr(ji,jj) = 0.
218                  ENDIF
219            ENDIF
220         END DO 
221      END DO 
222
223   END SUBROUTINE sbc_dcy
224
225   !!======================================================================
226END MODULE sbcdcy
Note: See TracBrowser for help on using the repository browser.