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/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/sbcdcy.F90 @ 3586

Last change on this file since 3586 was 3586, checked in by cbricaud, 11 years ago

add modification from dev_r3342_MERCATOR7_SST in dev_MERCATOR_2012_rev3555
  • Property svn:keywords set to Id
File size: 11.4 KB
Line 
1MODULE sbcdcy
2   !!======================================================================
3   !!                    ***  MODULE  sbcdcy  ***
4   !! Ocean forcing:  compute the diurnal cycle
5   !!======================================================================
6   !! History : OPA  !  2005-02  (D. Bernie)  Original code
7   !!   NEMO    2.0  !  2006-02  (S. Masson, G. Madec)  adaptation to NEMO
8   !!           3.1  !  2009-07  (J.M. Molines)  adaptation to v3.1
9   !!----------------------------------------------------------------------
10
11   !!----------------------------------------------------------------------
12   !!  sbc_dcy : solar flux at kt from daily mean, taking diurnal cycle into account
13   !!----------------------------------------------------------------------
14   USE oce              ! ocean dynamics and tracers
15   USE phycst           ! ocean physics
16   USE dom_oce          ! ocean space and time domain
17   USE sbc_oce          ! Surface boundary condition: ocean fields
18   USE in_out_manager   ! I/O manager
19   USE lib_mpp          ! MPP library
20   USE timing           ! Timing
21
22   IMPLICIT NONE
23   PRIVATE
24   
25   INTEGER, PUBLIC ::   nday_qsr   !: day when parameters were computed
26   
27   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   raa , rbb  , rcc  , rab     ! diurnal cycle parameters
28   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rtmd, rdawn, rdusk, rscal   !    -      -       -
29 
30   PUBLIC   sbc_dcy        ! routine called by sbc
31
32   !!----------------------------------------------------------------------
33   !! NEMO/OPA 3.3 , NEMO-consortium (2010)
34   !! $Id$
35   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
36   !!----------------------------------------------------------------------
37CONTAINS
38
39      INTEGER FUNCTION sbc_dcy_alloc()
40         !!----------------------------------------------------------------------
41         !!                ***  FUNCTION sbc_dcy_alloc  ***
42         !!----------------------------------------------------------------------
43         ALLOCATE( raa (jpi,jpj) , rbb  (jpi,jpj) , rcc  (jpi,jpj) , rab  (jpi,jpj) ,     &
44            &      rtmd(jpi,jpj) , rdawn(jpi,jpj) , rdusk(jpi,jpj) , rscal(jpi,jpj) , STAT=sbc_dcy_alloc )
45            !
46         IF( lk_mpp             )   CALL mpp_sum ( sbc_dcy_alloc )
47         IF( sbc_dcy_alloc /= 0 )   CALL ctl_warn('sbc_dcy_alloc: failed to allocate arrays')
48      END FUNCTION sbc_dcy_alloc
49
50
51   FUNCTION sbc_dcy( pqsrin, l_mask ) RESULT( zqsrout )
52      !!----------------------------------------------------------------------
53      !!                  ***  ROUTINE sbc_dcy  ***
54      !!
55      !! ** Purpose : introduce a diurnal cycle of qsr from daily values
56      !!
57      !! ** Method  : see Appendix A of Bernie et al. 2007.
58      !!
59      !! ** Action  : redistribute daily QSR on each time step following the diurnal cycle
60      !!
61      !! reference  : Bernie, DJ, E Guilyardi, G Madec, JM Slingo, and SJ Woolnough, 2007
62      !!              Impact of resolving the diurnal cycle in an ocean--atmosphere GCM.
63      !!              Part 1: a diurnally forced OGCM. Climate Dynamics 29:6, 575-590.
64      !!----------------------------------------------------------------------
65      LOGICAL, OPTIONAL, INTENT(in) :: l_mask ! use the routine for night mask computation
66      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqsrin    ! input daily QSR flux
67      !!
68      INTEGER  ::   ji, jj                                       ! dummy loop indices
69      INTEGER, DIMENSION(jpi,jpj) :: imask_night ! night mask
70      REAL(wp) ::   ztwopi, zinvtwopi, zconvrad 
71      REAL(wp) ::   zlo, zup, zlousd, zupusd
72      REAL(wp) ::   zdsws, zdecrad, ztx, zsin, zcos
73      REAL(wp) ::   ztmp, ztmp1, ztmp2, ztest
74      REAL(wp) ::   ztmpm, ztmpm1, ztmpm2
75      REAL(wp), DIMENSION(jpi,jpj) ::   zqsrout                  ! output QSR flux with diurnal cycle
76      !---------------------------statement functions------------------------
77      REAL(wp) ::   fintegral, pt1, pt2, paaa, pbbb, pccc        ! dummy statement function arguments
78      fintegral( pt1, pt2, paaa, pbbb, pccc ) =                         &
79         &   paaa * pt2 + zinvtwopi * pbbb * SIN(pccc + ztwopi * pt2)   &
80         & - paaa * pt1 - zinvtwopi * pbbb * SIN(pccc + ztwopi * pt1)
81      !!---------------------------------------------------------------------
82      !
83      IF( nn_timing == 1 )  CALL timing_start('sbc_dcy')
84      !
85      ! Initialization
86      ! --------------
87      ztwopi    = 2._wp * rpi
88      zinvtwopi = 1._wp / ztwopi
89      zconvrad  = ztwopi / 360._wp
90
91      ! When are we during the day (from 0 to 1)
92      zlo = ( REAL(nsec_day, wp) - 0.5_wp * rdttra(1) ) / rday
93      zup = zlo + ( REAL(nn_fsbc, wp)     * rdttra(1) ) / rday
94      !                                         
95      IF( nday_qsr == -1 ) THEN       ! first time step only 
96         IF(lwp) THEN
97            WRITE(numout,*)
98            WRITE(numout,*) 'sbc_dcy : introduce diurnal cycle from daily mean qsr'
99            WRITE(numout,*) '~~~~~~~'
100            WRITE(numout,*)
101         ENDIF
102         ! allocate sbcdcy arrays
103         IF( sbc_dcy_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_dcy_alloc : unable to allocate arrays' )
104         ! Compute rcc needed to compute the time integral of the diurnal cycle
105         rcc(:,:) = zconvrad * glamt(:,:) - rpi
106         ! time of midday
107         rtmd(:,:) = 0.5 - glamt(:,:) / 360.
108         rtmd(:,:) = MOD( (rtmd(:,:) + 1.), 1. )
109      ENDIF
110
111      ! If this is a new day, we have to update the dawn, dusk and scaling function 
112      !----------------------
113   
114      !     2.1 dawn and dusk 
115
116      ! nday is the number of days since the beginning of the current month
117      IF( nday_qsr /= nday ) THEN 
118         ! save the day of the year and the daily mean of qsr
119         nday_qsr = nday 
120         ! number of days since the previous winter solstice (supposed to be always 21 December)         
121         zdsws = REAL(11 + nday_year, wp)
122         ! declination of the earths orbit
123         zdecrad = (-23.5 * zconvrad) * COS( zdsws * ztwopi / REAL(nyear_len(1),wp) )
124         ! Compute A and B needed to compute the time integral of the diurnal cycle
125
126         zsin = SIN( zdecrad )   ;   zcos = COS( zdecrad )
127         DO jj = 1, jpj
128            DO ji = 1, jpi
129               ztmp = zconvrad * gphit(ji,jj)
130               raa(ji,jj) = SIN( ztmp ) * zsin
131               rbb(ji,jj) = COS( ztmp ) * zcos
132            END DO 
133         END DO 
134         ! Compute the time of dawn and dusk
135
136         ! rab to test if the day time is equal to 0, less than 24h of full day       
137         rab(:,:) = -raa(:,:) / rbb(:,:)
138         DO jj = 1, jpj
139            DO ji = 1, jpi
140               IF ( ABS(rab(ji,jj)) < 1 ) THEN         ! day duration is less than 24h
141         ! When is it night?
142                  ztx = zinvtwopi * (ACOS(rab(ji,jj)) - rcc(ji,jj))
143                  ztest = -rbb(ji,jj) * SIN( rcc(ji,jj) + ztwopi * ztx )
144         ! is it dawn or dusk?
145                  IF ( ztest > 0 ) THEN
146                     rdawn(ji,jj) = ztx
147                     rdusk(ji,jj) = rtmd(ji,jj) + ( rtmd(ji,jj) - rdawn(ji,jj) )
148                  ELSE
149                     rdusk(ji,jj) = ztx
150                     rdawn(ji,jj) = rtmd(ji,jj) - ( rdusk(ji,jj) - rtmd(ji,jj) )
151                  ENDIF
152               ELSE
153                  rdawn(ji,jj) = rtmd(ji,jj) + 0.5
154                  rdusk(ji,jj) = rdawn(ji,jj)
155               ENDIF
156             END DO 
157         END DO 
158         rdawn(:,:) = MOD( (rdawn(:,:) + 1._wp), 1._wp )
159         rdusk(:,:) = MOD( (rdusk(:,:) + 1._wp), 1._wp )
160         !     2.2 Compute the scalling function:
161         !         S* = the inverse of the time integral of the diurnal cycle from dawm to dusk
162         DO jj = 1, jpj
163            DO ji = 1, jpi
164               IF ( ABS(rab(ji,jj)) < 1 ) THEN         ! day duration is less than 24h
165                  IF ( rdawn(ji,jj) < rdusk(ji,jj) ) THEN      ! day time in one part
166                     rscal(ji,jj) = fintegral(rdawn(ji,jj), rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
167                     rscal(ji,jj) = 1. / rscal(ji,jj)
168                  ELSE                                         ! day time in two parts
169                     rscal(ji,jj) = fintegral(0., rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))   &
170                        &         + fintegral(rdawn(ji,jj), 1., raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
171                     rscal(ji,jj) = 1. / rscal(ji,jj)
172                  ENDIF
173               ELSE
174                  IF ( raa(ji,jj) > rbb(ji,jj) ) THEN         ! 24h day
175                     rscal(ji,jj) = fintegral(0., 1., raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
176                     rscal(ji,jj) = 1. / rscal(ji,jj)
177                  ELSE                                          ! No day
178                     rscal(ji,jj) = 0.e0
179                  ENDIF
180               ENDIF
181            END DO 
182         END DO 
183         !
184         ztmp = rday / ( rdttra(1) * REAL(nn_fsbc, wp) )
185         rscal(:,:) = rscal(:,:) * ztmp
186         !
187      ENDIF 
188         !     3. update qsr with the diurnal cycle
189         !     ------------------------------------
190
191      imask_night(:,:) = 0
192      DO jj = 1, jpj
193         DO ji = 1, jpi
194            ztmpm = 0.0
195            IF( ABS(rab(ji,jj)) < 1 ) THEN         ! day duration is less than 24h
196               !
197               IF( rdawn(ji,jj) < rdusk(ji,jj) ) THEN       ! day time in one part
198                  zlousd = MAX(zlo, rdawn(ji,jj))
199                  zlousd = MIN(zlousd, zup)
200                  zupusd = MIN(zup, rdusk(ji,jj))
201                  zupusd = MAX(zupusd, zlo)
202                  ztmp = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
203                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
204                  ztmpm = zupusd - zlousd
205                  IF ( ztmpm .EQ. 0 ) imask_night(ji,jj) = 1
206                  !
207               ELSE                                         ! day time in two parts
208                  zlousd = MIN(zlo, rdusk(ji,jj))
209                  zupusd = MIN(zup, rdusk(ji,jj))
210                  ztmp1 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
211                  ztmpm1=zupusd-zlousd
212                  zlousd = MAX(zlo, rdawn(ji,jj))
213                  zupusd = MAX(zup, rdawn(ji,jj))
214                  ztmp2 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
215                  ztmpm2 =zupusd-zlousd
216                  ztmp = ztmp1 + ztmp2
217                  ztmpm = ztmpm1 + ztmpm2
218                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
219                  IF (ztmpm .EQ. 0.) imask_night(ji,jj) = 1
220               ENDIF
221            ELSE                                   ! 24h light or 24h night
222               !
223               IF( raa(ji,jj) > rbb(ji,jj) ) THEN           ! 24h day
224                  ztmp = fintegral(zlo, zup, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
225                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
226                  imask_night(ji,jj) = 0
227                  !
228               ELSE                                         ! No day
229                  zqsrout(ji,jj) = 0.e0
230                  imask_night(ji,jj) = 1
231               ENDIF
232            ENDIF
233         END DO 
234      END DO 
235      !
236      IF ( PRESENT(l_mask) .AND. l_mask ) THEN
237         zqsrout(:,:) = float(imask_night(:,:))
238      ENDIF
239      !
240      IF( nn_timing == 1 )  CALL timing_stop('sbc_dcy')
241      !
242   END FUNCTION sbc_dcy
243
244   !!======================================================================
245END MODULE sbcdcy
Note: See TracBrowser for help on using the repository browser.