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 NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC – NEMO

source: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcdcy.F90 @ 12063

Last change on this file since 12063 was 12015, checked in by gsamson, 4 years ago

dev_ASINTER-01-05_merged: merge dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk branch @rev11988 with dev_r11265_ASINTER-01_Guillaume_ABL1D branch @rev11937 (tickets #2159 and #2131); ORCA2_ICE(_ABL) reproductibility OK

  • Property svn:keywords set to Id
File size: 12.8 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   !!           4.*  !  2019-10  (L. Brodeau)  nothing really new, but the routine
10   !!                ! "sbc_dcy_param" has been extracted from old function "sbc_dcy"
11   !!                ! => this allows the warm-layer param of COARE3* to know the time
12   !!                ! of dawn and dusk even if "ln_dm2dc=.false." (rdawn_dcy & rdusk_dcy
13   !!                ! are now public)
14   !!----------------------------------------------------------------------
15
16   !!----------------------------------------------------------------------
17   !!  sbc_dcy : solar flux at kt from daily mean, taking diurnal cycle into account
18   !!----------------------------------------------------------------------
19   USE oce              ! ocean dynamics and tracers
20   USE phycst           ! ocean physics
21   USE dom_oce          ! ocean space and time domain
22   USE sbc_oce          ! Surface boundary condition: ocean fields
23   !
24   USE in_out_manager   ! I/O manager
25   USE lib_mpp          ! MPP library
26
27   IMPLICIT NONE
28   PRIVATE
29
30   INTEGER, PUBLIC ::   nday_qsr   !: day when parameters were computed
31
32   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   raa , rbb  , rcc  , rab     ! diurnal cycle parameters
33   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rtmd, rscal   !    -      -       -
34   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: rdawn_dcy, rdusk_dcy   !    -      -       -
35
36   PUBLIC   sbc_dcy        ! routine called by sbc
37   PUBLIC   sbc_dcy_param  ! routine used here and called by warm-layer parameterization (sbcblk_skin_coare*)
38
39   !!----------------------------------------------------------------------
40   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
41   !! $Id$
42   !! Software governed by the CeCILL license (see ./LICENSE)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   INTEGER FUNCTION sbc_dcy_alloc()
47      !!----------------------------------------------------------------------
48      !!                ***  FUNCTION sbc_dcy_alloc  ***
49      !!----------------------------------------------------------------------
50      ALLOCATE( raa (jpi,jpj) , rbb  (jpi,jpj) , rcc  (jpi,jpj) , rab  (jpi,jpj) ,     &
51         &      rtmd(jpi,jpj) , rdawn_dcy(jpi,jpj) , rdusk_dcy(jpi,jpj) , rscal(jpi,jpj) , STAT=sbc_dcy_alloc )
52      !
53      CALL mpp_sum ( 'sbcdcy', sbc_dcy_alloc )
54      IF( sbc_dcy_alloc /= 0 )   CALL ctl_stop( 'STOP', 'sbc_dcy_alloc: failed to allocate arrays' )
55   END FUNCTION sbc_dcy_alloc
56
57
58   FUNCTION sbc_dcy( pqsrin, l_mask ) RESULT( zqsrout )
59      !!----------------------------------------------------------------------
60      !!                  ***  ROUTINE sbc_dcy  ***
61      !!
62      !! ** Purpose : introduce a diurnal cycle of qsr from daily values
63      !!
64      !! ** Method  : see Appendix A of Bernie et al. 2007.
65      !!
66      !! ** Action  : redistribute daily QSR on each time step following the diurnal cycle
67      !!
68      !! reference  : Bernie, DJ, E Guilyardi, G Madec, JM Slingo, and SJ Woolnough, 2007
69      !!              Impact of resolving the diurnal cycle in an ocean--atmosphere GCM.
70      !!              Part 1: a diurnally forced OGCM. Climate Dynamics 29:6, 575-590.
71      !!----------------------------------------------------------------------
72      LOGICAL , OPTIONAL          , INTENT(in) ::   l_mask    ! use the routine for night mask computation
73      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqsrin    ! input daily QSR flux
74      REAL(wp), DIMENSION(jpi,jpj)             ::   zqsrout   ! output QSR flux with diurnal cycle
75      !!
76      INTEGER  ::   ji, jj                                       ! dummy loop indices
77      INTEGER, DIMENSION(jpi,jpj) :: imask_night ! night mask
78      REAL(wp) ::   zlo, zup, zlousd, zupusd
79      REAL(wp) ::   ztmp, ztmp1, ztmp2
80      REAL(wp) ::   ztmpm, ztmpm1, ztmpm2
81      !!---------------------------------------------------------------------
82      !
83      ! Initialization
84      ! --------------
85      ! When are we during the day (from 0 to 1)
86      zlo = ( REAL(nsec_day, wp) - 0.5_wp * rdt ) / rday
87      zup = zlo + ( REAL(nn_fsbc, wp)     * rdt ) / rday
88      !
89      IF( nday_qsr == -1 ) THEN       ! first time step only
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
96      ENDIF
97
98      ! Setting parameters for each new day:
99      CALL sbc_dcy_param()
100
101      !CALL iom_put( "rdusk_dcy", rdusk_dcy(:,:)*tmask(:,:,1) ) !LB
102      !CALL iom_put( "rdawn_dcy", rdawn_dcy(:,:)*tmask(:,:,1) ) !LB
103      !CALL iom_put( "rscal_dcy", rscal(:,:)*tmask(:,:,1) ) !LB
104
105
106      !     3. update qsr with the diurnal cycle
107      !     ------------------------------------
108
109      imask_night(:,:) = 0
110      DO jj = 1, jpj
111         DO ji = 1, jpi
112            ztmpm = 0._wp
113            IF( ABS(rab(ji,jj)) < 1. ) THEN         ! day duration is less than 24h
114               !
115               IF( rdawn_dcy(ji,jj) < rdusk_dcy(ji,jj) ) THEN       ! day time in one part
116                  zlousd = MAX(zlo, rdawn_dcy(ji,jj))
117                  zlousd = MIN(zlousd, zup)
118                  zupusd = MIN(zup, rdusk_dcy(ji,jj))
119                  zupusd = MAX(zupusd, zlo)
120                  ztmp = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))
121                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
122                  ztmpm = zupusd - zlousd
123                  IF( ztmpm .EQ. 0 ) imask_night(ji,jj) = 1
124                  !
125               ELSE                                         ! day time in two parts
126                  zlousd = MIN(zlo, rdusk_dcy(ji,jj))
127                  zupusd = MIN(zup, rdusk_dcy(ji,jj))
128                  ztmp1 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))
129                  ztmpm1=zupusd-zlousd
130                  zlousd = MAX(zlo, rdawn_dcy(ji,jj))
131                  zupusd = MAX(zup, rdawn_dcy(ji,jj))
132                  ztmp2 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))
133                  ztmpm2 =zupusd-zlousd
134                  ztmp = ztmp1 + ztmp2
135                  ztmpm = ztmpm1 + ztmpm2
136                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
137                  IF(ztmpm .EQ. 0.) imask_night(ji,jj) = 1
138               ENDIF
139            ELSE                                   ! 24h light or 24h night
140               !
141               IF( raa(ji,jj) > rbb(ji,jj) ) THEN           ! 24h day
142                  ztmp = fintegral(zlo, zup, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))
143                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
144                  imask_night(ji,jj) = 0
145                  !
146               ELSE                                         ! No day
147                  zqsrout(ji,jj) = 0.0_wp
148                  imask_night(ji,jj) = 1
149               ENDIF
150            ENDIF
151         END DO
152      END DO
153      !
154      IF( PRESENT(l_mask) .AND. l_mask ) THEN
155         zqsrout(:,:) = float(imask_night(:,:))
156      ENDIF
157      !
158   END FUNCTION sbc_dcy
159
160
161   SUBROUTINE sbc_dcy_param( )
162      !!
163      INTEGER  ::   ji, jj                                       ! dummy loop indices
164      !INTEGER, DIMENSION(jpi,jpj) :: imask_night ! night mask
165      REAL(wp) ::   zdsws, zdecrad, ztx, zsin, zcos
166      REAL(wp) ::   ztmp, ztest
167      !---------------------------statement functions------------------------
168      !
169      IF( nday_qsr == -1 ) THEN       ! first time step only
170         ! allocate sbcdcy arrays
171         IF( sbc_dcy_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_dcy_alloc : unable to allocate arrays' )
172         ! Compute rcc needed to compute the time integral of the diurnal cycle
173         rcc(:,:) = rad * glamt(:,:) - rpi
174         ! time of midday
175         rtmd(:,:) = 0.5_wp - glamt(:,:) / 360._wp
176         rtmd(:,:) = MOD( (rtmd(:,:) + 1._wp) , 1._wp)
177      ENDIF
178
179      ! If this is a new day, we have to update the dawn, dusk and scaling function
180      !----------------------
181
182      !     2.1 dawn and dusk
183
184      ! nday is the number of days since the beginning of the current month
185      IF( nday_qsr /= nday ) THEN
186         ! save the day of the year and the daily mean of qsr
187         nday_qsr = nday
188         ! number of days since the previous winter solstice (supposed to be always 21 December)
189         zdsws = REAL(11 + nday_year, wp)
190         ! declination of the earths orbit
191         zdecrad = (-23.5_wp * rad) * COS( zdsws * 2._wp*rpi / REAL(nyear_len(1),wp) )
192         ! Compute A and B needed to compute the time integral of the diurnal cycle
193
194         zsin = SIN( zdecrad )   ;   zcos = COS( zdecrad )
195         DO jj = 1, jpj
196            DO ji = 1, jpi
197               ztmp = rad * gphit(ji,jj)
198               raa(ji,jj) = SIN( ztmp ) * zsin
199               rbb(ji,jj) = COS( ztmp ) * zcos
200            END DO
201         END DO
202         ! Compute the time of dawn and dusk
203
204         ! rab to test if the day time is equal to 0, less than 24h of full day
205         rab(:,:) = -raa(:,:) / rbb(:,:)
206         DO jj = 1, jpj
207            DO ji = 1, jpi
208               IF( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h
209                  ! When is it night?
210                  ztx = 1._wp/(2._wp*rpi) * (ACOS(rab(ji,jj)) - rcc(ji,jj))
211                  ztest = -rbb(ji,jj) * SIN( rcc(ji,jj) + 2._wp*rpi * ztx )
212                  ! is it dawn or dusk?
213                  IF( ztest > 0._wp ) THEN
214                     rdawn_dcy(ji,jj) = ztx
215                     rdusk_dcy(ji,jj) = rtmd(ji,jj) + ( rtmd(ji,jj) - rdawn_dcy(ji,jj) )
216                  ELSE
217                     rdusk_dcy(ji,jj) = ztx
218                     rdawn_dcy(ji,jj) = rtmd(ji,jj) - ( rdusk_dcy(ji,jj) - rtmd(ji,jj) )
219                  ENDIF
220               ELSE
221                  rdawn_dcy(ji,jj) = rtmd(ji,jj) + 0.5_wp
222                  rdusk_dcy(ji,jj) = rdawn_dcy(ji,jj)
223               ENDIF
224            END DO
225         END DO
226         rdawn_dcy(:,:) = MOD( (rdawn_dcy(:,:) + 1._wp), 1._wp )
227         rdusk_dcy(:,:) = MOD( (rdusk_dcy(:,:) + 1._wp), 1._wp )
228         !     2.2 Compute the scaling function:
229         !         S* = the inverse of the time integral of the diurnal cycle from dawn to dusk
230         !         Avoid possible infinite scaling factor, associated with very short daylight
231         !         periods, by ignoring periods less than 1/1000th of a day (ticket #1040)
232         DO jj = 1, jpj
233            DO ji = 1, jpi
234               IF( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h
235                  rscal(ji,jj) = 0.0_wp
236                  IF( rdawn_dcy(ji,jj) < rdusk_dcy(ji,jj) ) THEN      ! day time in one part
237                     IF( (rdusk_dcy(ji,jj) - rdawn_dcy(ji,jj) ) .ge. 0.001_wp ) THEN
238                        rscal(ji,jj) = fintegral(rdawn_dcy(ji,jj), rdusk_dcy(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))
239                        rscal(ji,jj) = 1._wp / rscal(ji,jj)
240                     ENDIF
241                  ELSE                                         ! day time in two parts
242                     IF( (rdusk_dcy(ji,jj) + (1._wp - rdawn_dcy(ji,jj)) ) .ge. 0.001_wp ) THEN
243                        rscal(ji,jj) = fintegral(0._wp, rdusk_dcy(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))   &
244                           &         + fintegral(rdawn_dcy(ji,jj), 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))
245                        rscal(ji,jj) = 1. / rscal(ji,jj)
246                     ENDIF
247                  ENDIF
248               ELSE
249                  IF( raa(ji,jj) > rbb(ji,jj) ) THEN         ! 24h day
250                     rscal(ji,jj) = fintegral(0._wp, 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))
251                     rscal(ji,jj) = 1._wp / rscal(ji,jj)
252                  ELSE                                          ! No day
253                     rscal(ji,jj) = 0.0_wp
254                  ENDIF
255               ENDIF
256            END DO
257         END DO
258         !
259         ztmp = rday / ( rdt * REAL(nn_fsbc, wp) )
260         rscal(:,:) = rscal(:,:) * ztmp
261         !
262      ENDIF !IF( nday_qsr /= nday )
263      !
264   END SUBROUTINE sbc_dcy_param
265
266
267   FUNCTION fintegral( pt1, pt2, paaa, pbbb, pccc )
268      REAL(wp), INTENT(in) :: pt1, pt2, paaa, pbbb, pccc
269      REAL(wp) :: fintegral
270      fintegral =   paaa * pt2 + 1._wp/(2._wp*rpi) * pbbb * SIN(pccc + 2._wp*rpi*pt2)   &
271         &        - paaa * pt1 - 1._wp/(2._wp*rpi) * pbbb * SIN(pccc + 2._wp*rpi*pt1)
272   END FUNCTION fintegral
273
274   !!======================================================================
275END MODULE sbcdcy
Note: See TracBrowser for help on using the repository browser.