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 @ 2216

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

diurnal cycle in coupled mode in dev_r2174_DCY, see ticket:730
File size: 9.7 KB
RevLine 
[2187]1MODULE sbcdcy
2   !!======================================================================
3   !!                    ***  MODULE  sbcdcy  ***
4   !! Ocean forcing:  compute the diurnal cycle
5   !!======================================================================
[2188]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
[2187]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
[2210]18   USE sbc_oce          ! Surface boundary condition: ocean fields
[2187]19   USE in_out_manager   ! I/O manager
20
21   IMPLICIT NONE
22   PRIVATE
[2216]23   INTEGER, PUBLIC              ::   nday_qsr                    ! day when parameters were computed
[2188]24   REAL(wp), DIMENSION(jpi,jpj) ::   raa , rbb  , rcc  , rab     ! parameters used to compute the diurnal cycle
25   REAL(wp), DIMENSION(jpi,jpj) ::   rtmd, rdawn, rdusk, rscal   !     -       -         -           -      -
[2187]26 
[2188]27   PUBLIC   sbc_dcy     ! routine called by sbc
[2187]28
29   !!----------------------------------------------------------------------
[2188]30   !! NEMO/OPA 3.3 , NEMO-consortium (2010)
[2187]31   !! $Id$
32   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
33   !!----------------------------------------------------------------------
34CONTAINS
35
[2216]36      FUNCTION sbc_dcy( pqsrin ) RESULT( zqsrout )
[2187]37      !!----------------------------------------------------------------------
38      !!                  ***  ROUTINE sbc_dcy  ***
39      !!
40      !! ** Purpose : introduce a diurnal cycle of qsr from daily values
41      !!
[2188]42      !! ** Method  : see Appendix A of Bernie et al. 2007.
[2187]43      !!
44      !! ** Action  : redistribute daily QSR on each time step following the diurnal cycle
[2188]45      !!
46      !! reference  : Bernie, DJ, E Guilyardi, G Madec, JM Slingo, and SJ Woolnough, 2007
47      !!              Impact of resolving the diurnal cycle in an ocean--atmosphere GCM.
48      !!              Part 1: a diurnally forced OGCM. Climate Dynamics 29:6, 575-590.
[2187]49      !!----------------------------------------------------------------------
[2210]50      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pqsrin    ! input daily QSR flux
[2187]51      !!
[2210]52      INTEGER  ::   ji, jj                                       ! dummy loop indices
[2187]53      REAL(wp) ::   ztwopi, zinvtwopi, zconvrad 
54      REAL(wp) ::   zlo, zup, zlousd, zupusd
[2210]55      REAL(wp) ::   zdsws, zdecrad, ztx, zsin, zcos
[2187]56      REAL(wp) ::   ztmp, ztmp1, ztmp2, ztest
[2210]57      REAL(wp), DIMENSION(jpi,jpj) ::   zqsrout                  ! output QSR flux with diurnal cycle
[2188]58      !---------------------------statement functions------------------------
[2210]59      REAL(wp) ::   fintegral, pt1, pt2, paaa, pbbb, pccc        ! dummy statement function arguments
[2188]60      fintegral( pt1, pt2, paaa, pbbb, pccc ) =                         &
[2187]61         &   paaa * pt2 + zinvtwopi * pbbb * SIN(pccc + ztwopi * pt2)   &
62         & - paaa * pt1 - zinvtwopi * pbbb * SIN(pccc + ztwopi * pt1)
[2188]63      !!---------------------------------------------------------------------
[2187]64
65      ! Initialization
66      ! --------------
67      ztwopi    = 2. * rpi
68      zinvtwopi = 1. / ztwopi
69      zconvrad  = ztwopi / 360.
70
71      ! When are we during the day (from 0 to 1)
[2210]72      zlo = ( REAL(nsec_day, wp) - 0.5   * rdttra(1) ) / rday
73      zup = zlo + ( REAL(nn_fsbc, wp) * rdttra(1) ) / rday
[2187]74
75      !                                         
[2216]76      IF( nday_qsr == -1 ) THEN       ! first time step only               
[2187]77         IF(lwp) THEN
78            WRITE(numout,*)
79            WRITE(numout,*) 'sbc_dcy : introduce diurnal cycle from daily mean qsr'
80            WRITE(numout,*) '~~~~~~~'
81            WRITE(numout,*)
82         ENDIF
[2188]83         ! Compute rcc needed to compute the time integral of the diurnal cycle
84         rcc(:,:) = zconvrad * glamt(:,:) - rpi
[2187]85         ! time of midday
[2188]86         rtmd(:,:) = 0.5 - glamt(:,:) / 360.
87         rtmd(:,:) = MOD( (rtmd(:,:) + 1.), 1. )
[2187]88      ENDIF
89
90      ! If this is a new day, we have to update the dawn, dusk and scaling function 
91      !----------------------
92   
93      !     2.1 dawn and dusk 
94
95      ! nday is the number of days since the beginning of the current month
[2188]96      IF( nday_qsr /= nday ) THEN 
[2187]97         ! save the day of the year and the daily mean of qsr
[2188]98         nday_qsr = nday 
[2187]99         ! number of days since the previous winter solstice (supposed to be always 21 December)         
[2210]100         zdsws = REAL(11 + nday_year, wp)
[2187]101         ! declination of the earths orbit
102         zdecrad = (-23.5 * zconvrad) * COS( zdsws * ztwopi / REAL(nyear_len(1),wp) )
103         ! Compute A and B needed to compute the time integral of the diurnal cycle
104       
[2210]105         zsin = SIN( zdecrad )   ;   zcos = COS( zdecrad )
[2187]106         DO jj = 1, jpj
107            DO ji = 1, jpi
108               ztmp = zconvrad * gphit(ji,jj)
[2210]109               raa(ji,jj) = SIN( ztmp ) * zsin
110               rbb(ji,jj) = COS( ztmp ) * zcos
[2187]111            END DO 
112         END DO 
113
114         ! Compute the time of dawn and dusk
115
[2188]116         ! rab to test if the day time is equal to 0, less than 24h of full day       
117         rab(:,:) = -raa(:,:) / rbb(:,:)
[2187]118         DO jj = 1, jpj
119            DO ji = 1, jpi
[2188]120               IF ( ABS(rab(ji,jj)) < 1 ) THEN         ! day duration is less than 24h
[2187]121         ! When is it night?
[2188]122                  ztx = zinvtwopi * (ACOS(rab(ji,jj)) - rcc(ji,jj))
123                  ztest = -rbb(ji,jj) * SIN( rcc(ji,jj) + ztwopi * ztx )
[2187]124         ! is it dawn or dusk?
125                  IF ( ztest > 0 ) THEN
[2188]126                     rdawn(ji,jj) = ztx
127                     rdusk(ji,jj) = rtmd(ji,jj) + ( rtmd(ji,jj) - rdawn(ji,jj) )
[2187]128                  ELSE
[2188]129                     rdusk(ji,jj) = ztx
130                     rdawn(ji,jj) = rtmd(ji,jj) - ( rdusk(ji,jj) - rtmd(ji,jj) )
[2187]131                  ENDIF
132               ELSE
[2188]133                  rdawn(ji,jj) = rtmd(ji,jj) + 0.5
134                  rdusk(ji,jj) = rdawn(ji,jj)
[2187]135               ENDIF
136             END DO 
137         END DO 
[2188]138         rdawn(:,:) = MOD((rdawn(:,:) + 1.), 1.)
139         rdusk(:,:) = MOD((rdusk(:,:) + 1.), 1.)
[2187]140
141         !     2.2 Compute the scalling function:
142         !         S* = the inverse of the time integral of the diurnal cycle from dawm to dusk
143         DO jj = 1, jpj
144            DO ji = 1, jpi
[2188]145               IF ( ABS(rab(ji,jj)) < 1 ) THEN         ! day duration is less than 24h
146                  IF ( rdawn(ji,jj) < rdusk(ji,jj) ) THEN      ! day time in one part
147                     rscal(ji,jj) = fintegral(rdawn(ji,jj), rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
148                     rscal(ji,jj) = 1. / rscal(ji,jj)
149                  ELSE                                         ! day time in two parts
150                     rscal(ji,jj) = fintegral(0., rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))   &
151                        &         + fintegral(rdawn(ji,jj), 1., raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
152                     rscal(ji,jj) = 1. / rscal(ji,jj)
[2187]153                  ENDIF
154               ELSE
[2188]155                  IF ( raa(ji,jj) > rbb(ji,jj) ) THEN         ! 24h day
156                     rscal(ji,jj) = fintegral(0., 1., raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
157                     rscal(ji,jj) = 1. / rscal(ji,jj)
158                  ELSE                                          ! No day
159                     rscal(ji,jj) = 0.e0
[2187]160                  ENDIF
161               ENDIF
162            END DO 
163         END DO 
164         !
[2210]165         ztmp = rday / ( rdttra(1) * REAL(nn_fsbc, wp) )
[2188]166         rscal(:,:) = rscal(:,:) * ztmp
[2187]167
168      ENDIF 
169
[2188]170         !     3. update qsr with the diurnal cycle
171         !     ------------------------------------
[2187]172
173      DO jj = 1, jpj
174         DO ji = 1, jpi
[2188]175            IF( ABS(rab(ji,jj)) < 1 ) THEN         ! day duration is less than 24h
176               !
177               IF( rdawn(ji,jj) < rdusk(ji,jj) ) THEN       ! day time in one part
178                  zlousd = MAX(zlo, rdawn(ji,jj))
179                  zlousd = MIN(zlousd, zup)
180                  zupusd = MIN(zup, rdusk(ji,jj))
181                  zupusd = MAX(zupusd, zlo)
182                  ztmp = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
[2210]183                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
[2188]184                  !
185               ELSE                                         ! day time in two parts
186                  zlousd = MIN(zlo, rdusk(ji,jj))
187                  zupusd = MIN(zup, rdusk(ji,jj))
188                  ztmp1 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
189                  zlousd = MAX(zlo, rdawn(ji,jj))
190                  zupusd = MAX(zup, rdawn(ji,jj))
191                  ztmp2 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
192                  ztmp = ztmp1 + ztmp2
[2210]193                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
[2188]194               ENDIF
195            ELSE                                   ! 24h light or 24h night
196               !
[2210]197               IF( raa(ji,jj) > rbb(ji,jj) ) THEN           ! 24h day
[2188]198                  ztmp = fintegral(zlo, zup, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
[2210]199                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
[2188]200                  !
201               ELSE                                         ! No day
[2210]202                  zqsrout(ji,jj) = 0.e0
[2188]203               ENDIF
[2187]204            ENDIF
205         END DO 
206      END DO 
[2188]207      !
[2210]208   END FUNCTION sbc_dcy
[2187]209
210   !!======================================================================
211END MODULE sbcdcy
Note: See TracBrowser for help on using the repository browser.