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/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcdcy.F90 @ 3270

Last change on this file since 3270 was 3211, checked in by spickles2, 13 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 10.7 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
21   IMPLICIT NONE
22   PRIVATE
23   
24   INTEGER, PUBLIC ::   nday_qsr   !: day when parameters were computed
25   
26   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   raa , rbb  , rcc  , rab     ! diurnal cycle parameters
27   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rtmd, rdawn, rdusk, rscal   !    -      -       -
28 
29   PUBLIC   sbc_dcy        ! routine called by sbc
30
31   !! * Control permutation of array indices
32#  include "oce_ftrans.h90"
33#  include "dom_oce_ftrans.h90"
34#  include "sbc_oce_ftrans.h90"
35
36   !!----------------------------------------------------------------------
37   !! NEMO/OPA 3.3 , NEMO-consortium (2010)
38   !! $Id$
39   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
40   !!----------------------------------------------------------------------
41CONTAINS
42
43      INTEGER FUNCTION sbc_dcy_alloc()
44         !!----------------------------------------------------------------------
45         !!                ***  FUNCTION sbc_dcy_alloc  ***
46         !!----------------------------------------------------------------------
47         ALLOCATE( raa (jpi,jpj) , rbb  (jpi,jpj) , rcc  (jpi,jpj) , rab  (jpi,jpj) ,     &
48            &      rtmd(jpi,jpj) , rdawn(jpi,jpj) , rdusk(jpi,jpj) , rscal(jpi,jpj) , STAT=sbc_dcy_alloc )
49            !
50         IF( lk_mpp             )   CALL mpp_sum ( sbc_dcy_alloc )
51         IF( sbc_dcy_alloc /= 0 )   CALL ctl_warn('sbc_dcy_alloc: failed to allocate arrays')
52      END FUNCTION sbc_dcy_alloc
53
54
55   FUNCTION sbc_dcy( pqsrin ) RESULT( zqsrout )
56      !!----------------------------------------------------------------------
57      !!                  ***  ROUTINE sbc_dcy  ***
58      !!
59      !! ** Purpose : introduce a diurnal cycle of qsr from daily values
60      !!
61      !! ** Method  : see Appendix A of Bernie et al. 2007.
62      !!
63      !! ** Action  : redistribute daily QSR on each time step following the diurnal cycle
64      !!
65      !! reference  : Bernie, DJ, E Guilyardi, G Madec, JM Slingo, and SJ Woolnough, 2007
66      !!              Impact of resolving the diurnal cycle in an ocean--atmosphere GCM.
67      !!              Part 1: a diurnally forced OGCM. Climate Dynamics 29:6, 575-590.
68      !!----------------------------------------------------------------------
69      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqsrin    ! input daily QSR flux
70      !!
71      INTEGER  ::   ji, jj                                       ! dummy loop indices
72      REAL(wp) ::   ztwopi, zinvtwopi, zconvrad 
73      REAL(wp) ::   zlo, zup, zlousd, zupusd
74      REAL(wp) ::   zdsws, zdecrad, ztx, zsin, zcos
75      REAL(wp) ::   ztmp, ztmp1, ztmp2, ztest
76      REAL(wp), DIMENSION(jpi,jpj) ::   zqsrout                  ! output QSR flux with diurnal cycle
77      !---------------------------statement functions------------------------
78      REAL(wp) ::   fintegral, pt1, pt2, paaa, pbbb, pccc        ! dummy statement function arguments
79      fintegral( pt1, pt2, paaa, pbbb, pccc ) =                         &
80         &   paaa * pt2 + zinvtwopi * pbbb * SIN(pccc + ztwopi * pt2)   &
81         & - paaa * pt1 - zinvtwopi * pbbb * SIN(pccc + ztwopi * pt1)
82      !!---------------------------------------------------------------------
83
84      ! Initialization
85      ! --------------
86      ztwopi    = 2._wp * rpi
87      zinvtwopi = 1._wp / ztwopi
88      zconvrad  = ztwopi / 360._wp
89
90      ! When are we during the day (from 0 to 1)
91      zlo = ( REAL(nsec_day, wp) - 0.5_wp * rdttra(1) ) / rday
92      zup = zlo + ( REAL(nn_fsbc, wp)     * rdttra(1) ) / rday
93      !                                         
94      IF( nday_qsr == -1 ) THEN       ! first time step only               
95         IF(lwp) THEN
96            WRITE(numout,*)
97            WRITE(numout,*) 'sbc_dcy : introduce diurnal cycle from daily mean qsr'
98            WRITE(numout,*) '~~~~~~~'
99            WRITE(numout,*)
100         ENDIF
101         ! allocate sbcdcy arrays
102         IF( sbc_dcy_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_dcy_alloc : unable to allocate arrays' )
103         ! Compute rcc needed to compute the time integral of the diurnal cycle
104         rcc(:,:) = zconvrad * glamt(:,:) - rpi
105         ! time of midday
106         rtmd(:,:) = 0.5 - glamt(:,:) / 360.
107         rtmd(:,:) = MOD( (rtmd(:,:) + 1.), 1. )
108      ENDIF
109
110      ! If this is a new day, we have to update the dawn, dusk and scaling function 
111      !----------------------
112   
113      !     2.1 dawn and dusk 
114
115      ! nday is the number of days since the beginning of the current month
116      IF( nday_qsr /= nday ) THEN 
117         ! save the day of the year and the daily mean of qsr
118         nday_qsr = nday 
119         ! number of days since the previous winter solstice (supposed to be always 21 December)         
120         zdsws = REAL(11 + nday_year, wp)
121         ! declination of the earths orbit
122         zdecrad = (-23.5 * zconvrad) * COS( zdsws * ztwopi / REAL(nyear_len(1),wp) )
123         ! Compute A and B needed to compute the time integral of the diurnal cycle
124       
125         zsin = SIN( zdecrad )   ;   zcos = COS( zdecrad )
126         DO jj = 1, jpj
127            DO ji = 1, jpi
128               ztmp = zconvrad * gphit(ji,jj)
129               raa(ji,jj) = SIN( ztmp ) * zsin
130               rbb(ji,jj) = COS( ztmp ) * zcos
131            END DO 
132         END DO 
133
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
161         !     2.2 Compute the scalling function:
162         !         S* = the inverse of the time integral of the diurnal cycle from dawm to dusk
163         DO jj = 1, jpj
164            DO ji = 1, jpi
165               IF ( ABS(rab(ji,jj)) < 1 ) THEN         ! day duration is less than 24h
166                  IF ( rdawn(ji,jj) < rdusk(ji,jj) ) THEN      ! day time in one part
167                     rscal(ji,jj) = fintegral(rdawn(ji,jj), rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
168                     rscal(ji,jj) = 1. / rscal(ji,jj)
169                  ELSE                                         ! day time in two parts
170                     rscal(ji,jj) = fintegral(0., rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))   &
171                        &         + fintegral(rdawn(ji,jj), 1., raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
172                     rscal(ji,jj) = 1. / rscal(ji,jj)
173                  ENDIF
174               ELSE
175                  IF ( raa(ji,jj) > rbb(ji,jj) ) THEN         ! 24h day
176                     rscal(ji,jj) = fintegral(0., 1., raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
177                     rscal(ji,jj) = 1. / rscal(ji,jj)
178                  ELSE                                          ! No day
179                     rscal(ji,jj) = 0.e0
180                  ENDIF
181               ENDIF
182            END DO 
183         END DO 
184         !
185         ztmp = rday / ( rdttra(1) * REAL(nn_fsbc, wp) )
186         rscal(:,:) = rscal(:,:) * ztmp
187         !
188      ENDIF 
189
190         !     3. update qsr with the diurnal cycle
191         !     ------------------------------------
192
193      DO jj = 1, jpj
194         DO ji = 1, jpi
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                  !
205               ELSE                                         ! day time in two parts
206                  zlousd = MIN(zlo, rdusk(ji,jj))
207                  zupusd = MIN(zup, rdusk(ji,jj))
208                  ztmp1 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
209                  zlousd = MAX(zlo, rdawn(ji,jj))
210                  zupusd = MAX(zup, rdawn(ji,jj))
211                  ztmp2 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
212                  ztmp = ztmp1 + ztmp2
213                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
214               ENDIF
215            ELSE                                   ! 24h light or 24h night
216               !
217               IF( raa(ji,jj) > rbb(ji,jj) ) THEN           ! 24h day
218                  ztmp = fintegral(zlo, zup, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
219                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
220                  !
221               ELSE                                         ! No day
222                  zqsrout(ji,jj) = 0.e0
223               ENDIF
224            ENDIF
225         END DO 
226      END DO 
227      !
228   END FUNCTION sbc_dcy
229
230   !!======================================================================
231END MODULE sbcdcy
Note: See TracBrowser for help on using the repository browser.