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

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcdcy.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.
File size: 11.9 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            IF(lflush) CALL flush(numout)
102         ENDIF
103         ! allocate sbcdcy arrays
104         IF( sbc_dcy_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_dcy_alloc : unable to allocate arrays' )
105         ! Compute rcc needed to compute the time integral of the diurnal cycle
106         rcc(:,:) = zconvrad * glamt(:,:) - rpi
107         ! time of midday
108         rtmd(:,:) = 0.5_wp - glamt(:,:) / 360._wp
109         rtmd(:,:) = MOD( (rtmd(:,:) + 1._wp) , 1._wp)
110      ENDIF
111
112      ! If this is a new day, we have to update the dawn, dusk and scaling function 
113      !----------------------
114   
115      !     2.1 dawn and dusk 
116
117      ! nday is the number of days since the beginning of the current month
118      IF( nday_qsr /= nday ) THEN 
119         ! save the day of the year and the daily mean of qsr
120         nday_qsr = nday 
121         ! number of days since the previous winter solstice (supposed to be always 21 December)         
122         zdsws = REAL(11 + nday_year, wp)
123         ! declination of the earths orbit
124         zdecrad = (-23.5_wp * zconvrad) * COS( zdsws * ztwopi / REAL(nyear_len(1),wp) )
125         ! Compute A and B needed to compute the time integral of the diurnal cycle
126
127         zsin = SIN( zdecrad )   ;   zcos = COS( zdecrad )
128         DO jj = 1, jpj
129            DO ji = 1, jpi
130               ztmp = zconvrad * gphit(ji,jj)
131               raa(ji,jj) = SIN( ztmp ) * zsin
132               rbb(ji,jj) = COS( ztmp ) * zcos
133            END DO 
134         END DO 
135         ! Compute the time of dawn and dusk
136
137         ! rab to test if the day time is equal to 0, less than 24h of full day       
138         rab(:,:) = -raa(:,:) / rbb(:,:)
139         DO jj = 1, jpj
140            DO ji = 1, jpi
141               IF ( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h
142         ! When is it night?
143                  ztx = zinvtwopi * (ACOS(rab(ji,jj)) - rcc(ji,jj))
144                  ztest = -rbb(ji,jj) * SIN( rcc(ji,jj) + ztwopi * ztx )
145         ! is it dawn or dusk?
146                  IF ( ztest > 0._wp ) THEN
147                     rdawn(ji,jj) = ztx
148                     rdusk(ji,jj) = rtmd(ji,jj) + ( rtmd(ji,jj) - rdawn(ji,jj) )
149                  ELSE
150                     rdusk(ji,jj) = ztx
151                     rdawn(ji,jj) = rtmd(ji,jj) - ( rdusk(ji,jj) - rtmd(ji,jj) )
152                  ENDIF
153               ELSE
154                  rdawn(ji,jj) = rtmd(ji,jj) + 0.5_wp
155                  rdusk(ji,jj) = rdawn(ji,jj)
156               ENDIF
157             END DO 
158         END DO 
159         rdawn(:,:) = MOD( (rdawn(:,:) + 1._wp), 1._wp )
160         rdusk(:,:) = MOD( (rdusk(:,:) + 1._wp), 1._wp )
161         !     2.2 Compute the scaling function:
162         !         S* = the inverse of the time integral of the diurnal cycle from dawn to dusk
163         !         Avoid possible infinite scaling factor, associated with very short daylight
164         !         periods, by ignoring periods less than 1/1000th of a day (ticket #1040)
165         DO jj = 1, jpj
166            DO ji = 1, jpi
167               IF ( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h
168                  rscal(ji,jj) = 0.0_wp
169                  IF ( rdawn(ji,jj) < rdusk(ji,jj) ) THEN      ! day time in one part
170                     IF( (rdusk(ji,jj) - rdawn(ji,jj) ) .ge. 0.001_wp ) THEN
171                       rscal(ji,jj) = fintegral(rdawn(ji,jj), rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
172                       rscal(ji,jj) = 1._wp / rscal(ji,jj)
173                     ENDIF
174                  ELSE                                         ! day time in two parts
175                     IF( (rdusk(ji,jj) + (1._wp - rdawn(ji,jj)) ) .ge. 0.001_wp ) THEN
176                       rscal(ji,jj) = fintegral(0._wp, rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))   &
177                          &         + fintegral(rdawn(ji,jj), 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
178                       rscal(ji,jj) = 1. / rscal(ji,jj)
179                     ENDIF
180                  ENDIF
181               ELSE
182                  IF ( raa(ji,jj) > rbb(ji,jj) ) THEN         ! 24h day
183                     rscal(ji,jj) = fintegral(0._wp, 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
184                     rscal(ji,jj) = 1._wp / rscal(ji,jj)
185                  ELSE                                          ! No day
186                     rscal(ji,jj) = 0.0_wp
187                  ENDIF
188               ENDIF
189            END DO 
190         END DO 
191         !
192         ztmp = rday / ( rdttra(1) * REAL(nn_fsbc, wp) )
193         rscal(:,:) = rscal(:,:) * ztmp
194         !
195      ENDIF 
196         !     3. update qsr with the diurnal cycle
197         !     ------------------------------------
198
199      imask_night(:,:) = 0
200      DO jj = 1, jpj
201         DO ji = 1, jpi
202            ztmpm = 0.0
203            IF( ABS(rab(ji,jj)) < 1. ) THEN         ! day duration is less than 24h
204               !
205               IF( rdawn(ji,jj) < rdusk(ji,jj) ) THEN       ! day time in one part
206                  zlousd = MAX(zlo, rdawn(ji,jj))
207                  zlousd = MIN(zlousd, zup)
208                  zupusd = MIN(zup, rdusk(ji,jj))
209                  zupusd = MAX(zupusd, zlo)
210                  ztmp = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
211                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
212                  ztmpm = zupusd - zlousd
213                  IF ( ztmpm .EQ. 0 ) imask_night(ji,jj) = 1
214                  !
215               ELSE                                         ! day time in two parts
216                  zlousd = MIN(zlo, rdusk(ji,jj))
217                  zupusd = MIN(zup, rdusk(ji,jj))
218                  ztmp1 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
219                  ztmpm1=zupusd-zlousd
220                  zlousd = MAX(zlo, rdawn(ji,jj))
221                  zupusd = MAX(zup, rdawn(ji,jj))
222                  ztmp2 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
223                  ztmpm2 =zupusd-zlousd
224                  ztmp = ztmp1 + ztmp2
225                  ztmpm = ztmpm1 + ztmpm2
226                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
227                  IF (ztmpm .EQ. 0.) imask_night(ji,jj) = 1
228               ENDIF
229            ELSE                                   ! 24h light or 24h night
230               !
231               IF( raa(ji,jj) > rbb(ji,jj) ) THEN           ! 24h day
232                  ztmp = fintegral(zlo, zup, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
233                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
234                  imask_night(ji,jj) = 0
235                  !
236               ELSE                                         ! No day
237                  zqsrout(ji,jj) = 0.0_wp
238                  imask_night(ji,jj) = 1
239               ENDIF
240            ENDIF
241         END DO 
242      END DO 
243      !
244      IF ( PRESENT(l_mask) .AND. l_mask ) THEN
245         zqsrout(:,:) = float(imask_night(:,:))
246      ENDIF
247      !
248      IF( nn_timing == 1 )  CALL timing_stop('sbc_dcy')
249      !
250   END FUNCTION sbc_dcy
251
252   !!======================================================================
253END MODULE sbcdcy
Note: See TracBrowser for help on using the repository browser.