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

source: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcdcy.F90 @ 2777

Last change on this file since 2777 was 2715, checked in by rblod, 13 years ago

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