source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/SBC/sbcdcy.F90 @ 10297

Last change on this file since 10297 was 10297, checked in by smasson, 2 years ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: action 2a: add report calls of mppmin/max/sum, see #2133

  • Property svn:keywords set to Id
File size: 11.7 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
[9124]18   !
[2198]19   USE in_out_manager   ! I/O manager
[2715]20   USE lib_mpp          ! MPP library
[2198]21
22   IMPLICIT NONE
23   PRIVATE
[2715]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   !    -      -       -
[2198]29 
[2715]30   PUBLIC   sbc_dcy        ! routine called by sbc
[2198]31
32   !!----------------------------------------------------------------------
[10068]33   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[2198]34   !! $Id$
[10068]35   !! Software governed by the CeCILL license (see ./LICENSE)
[2198]36   !!----------------------------------------------------------------------
37CONTAINS
38
[2715]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            !
[10297]46         IF( lk_mpp             )   CALL mpp_sum ( 'sbcdcy', sbc_dcy_alloc )
[2715]47         IF( sbc_dcy_alloc /= 0 )   CALL ctl_warn('sbc_dcy_alloc: failed to allocate arrays')
48      END FUNCTION sbc_dcy_alloc
49
50
[3651]51   FUNCTION sbc_dcy( pqsrin, l_mask ) RESULT( zqsrout )
[2198]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      !!----------------------------------------------------------------------
[9124]65      LOGICAL , OPTIONAL          , INTENT(in) ::   l_mask    ! use the routine for night mask computation
[2228]66      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqsrin    ! input daily QSR flux
[9124]67      REAL(wp), DIMENSION(jpi,jpj)             ::   zqsrout   ! output QSR flux with diurnal cycle
[2198]68      !!
[2228]69      INTEGER  ::   ji, jj                                       ! dummy loop indices
[3651]70      INTEGER, DIMENSION(jpi,jpj) :: imask_night ! night mask
[2198]71      REAL(wp) ::   ztwopi, zinvtwopi, zconvrad 
72      REAL(wp) ::   zlo, zup, zlousd, zupusd
[2228]73      REAL(wp) ::   zdsws, zdecrad, ztx, zsin, zcos
[2198]74      REAL(wp) ::   ztmp, ztmp1, ztmp2, ztest
[3651]75      REAL(wp) ::   ztmpm, ztmpm1, ztmpm2
[2198]76      !---------------------------statement functions------------------------
[2228]77      REAL(wp) ::   fintegral, pt1, pt2, paaa, pbbb, pccc        ! dummy statement function arguments
[2198]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      !!---------------------------------------------------------------------
[3294]82      !
[2198]83      ! Initialization
84      ! --------------
[2715]85      ztwopi    = 2._wp * rpi
86      zinvtwopi = 1._wp / ztwopi
87      zconvrad  = ztwopi / 360._wp
[2198]88
89      ! When are we during the day (from 0 to 1)
[6140]90      zlo = ( REAL(nsec_day, wp) - 0.5_wp * rdt ) / rday
91      zup = zlo + ( REAL(nn_fsbc, wp)     * rdt ) / rday
[2198]92      !                                         
[3651]93      IF( nday_qsr == -1 ) THEN       ! first time step only 
[2198]94         IF(lwp) THEN
95            WRITE(numout,*)
96            WRITE(numout,*) 'sbc_dcy : introduce diurnal cycle from daily mean qsr'
97            WRITE(numout,*) '~~~~~~~'
98            WRITE(numout,*)
99         ENDIF
[2715]100         ! allocate sbcdcy arrays
101         IF( sbc_dcy_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_dcy_alloc : unable to allocate arrays' )
[2198]102         ! Compute rcc needed to compute the time integral of the diurnal cycle
103         rcc(:,:) = zconvrad * glamt(:,:) - rpi
104         ! time of midday
[3764]105         rtmd(:,:) = 0.5_wp - glamt(:,:) / 360._wp
106         rtmd(:,:) = MOD( (rtmd(:,:) + 1._wp) , 1._wp)
[2198]107      ENDIF
108
109      ! If this is a new day, we have to update the dawn, dusk and scaling function 
110      !----------------------
111   
112      !     2.1 dawn and dusk 
113
114      ! nday is the number of days since the beginning of the current month
115      IF( nday_qsr /= nday ) THEN 
116         ! save the day of the year and the daily mean of qsr
117         nday_qsr = nday 
118         ! number of days since the previous winter solstice (supposed to be always 21 December)         
[2228]119         zdsws = REAL(11 + nday_year, wp)
[2198]120         ! declination of the earths orbit
[3764]121         zdecrad = (-23.5_wp * zconvrad) * COS( zdsws * ztwopi / REAL(nyear_len(1),wp) )
[2198]122         ! Compute A and B needed to compute the time integral of the diurnal cycle
[3651]123
[2228]124         zsin = SIN( zdecrad )   ;   zcos = COS( zdecrad )
[2198]125         DO jj = 1, jpj
126            DO ji = 1, jpi
127               ztmp = zconvrad * gphit(ji,jj)
[2228]128               raa(ji,jj) = SIN( ztmp ) * zsin
129               rbb(ji,jj) = COS( ztmp ) * zcos
[2198]130            END DO 
131         END DO 
132         ! Compute the time of dawn and dusk
133
134         ! rab to test if the day time is equal to 0, less than 24h of full day       
135         rab(:,:) = -raa(:,:) / rbb(:,:)
136         DO jj = 1, jpj
137            DO ji = 1, jpi
[3764]138               IF ( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h
[2198]139         ! When is it night?
140                  ztx = zinvtwopi * (ACOS(rab(ji,jj)) - rcc(ji,jj))
141                  ztest = -rbb(ji,jj) * SIN( rcc(ji,jj) + ztwopi * ztx )
142         ! is it dawn or dusk?
[3764]143                  IF ( ztest > 0._wp ) THEN
[2198]144                     rdawn(ji,jj) = ztx
145                     rdusk(ji,jj) = rtmd(ji,jj) + ( rtmd(ji,jj) - rdawn(ji,jj) )
146                  ELSE
147                     rdusk(ji,jj) = ztx
148                     rdawn(ji,jj) = rtmd(ji,jj) - ( rdusk(ji,jj) - rtmd(ji,jj) )
149                  ENDIF
150               ELSE
[3764]151                  rdawn(ji,jj) = rtmd(ji,jj) + 0.5_wp
[2198]152                  rdusk(ji,jj) = rdawn(ji,jj)
153               ENDIF
154             END DO 
155         END DO 
[2715]156         rdawn(:,:) = MOD( (rdawn(:,:) + 1._wp), 1._wp )
157         rdusk(:,:) = MOD( (rdusk(:,:) + 1._wp), 1._wp )
[3764]158         !     2.2 Compute the scaling function:
159         !         S* = the inverse of the time integral of the diurnal cycle from dawn to dusk
160         !         Avoid possible infinite scaling factor, associated with very short daylight
161         !         periods, by ignoring periods less than 1/1000th of a day (ticket #1040)
[2198]162         DO jj = 1, jpj
163            DO ji = 1, jpi
[3764]164               IF ( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h
165                  rscal(ji,jj) = 0.0_wp
[2198]166                  IF ( rdawn(ji,jj) < rdusk(ji,jj) ) THEN      ! day time in one part
[3764]167                     IF( (rdusk(ji,jj) - rdawn(ji,jj) ) .ge. 0.001_wp ) THEN
168                       rscal(ji,jj) = fintegral(rdawn(ji,jj), rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
169                       rscal(ji,jj) = 1._wp / rscal(ji,jj)
170                     ENDIF
[2198]171                  ELSE                                         ! day time in two parts
[3764]172                     IF( (rdusk(ji,jj) + (1._wp - rdawn(ji,jj)) ) .ge. 0.001_wp ) THEN
173                       rscal(ji,jj) = fintegral(0._wp, rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))   &
174                          &         + fintegral(rdawn(ji,jj), 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
175                       rscal(ji,jj) = 1. / rscal(ji,jj)
176                     ENDIF
[2198]177                  ENDIF
178               ELSE
179                  IF ( raa(ji,jj) > rbb(ji,jj) ) THEN         ! 24h day
[3764]180                     rscal(ji,jj) = fintegral(0._wp, 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
181                     rscal(ji,jj) = 1._wp / rscal(ji,jj)
[2198]182                  ELSE                                          ! No day
[3764]183                     rscal(ji,jj) = 0.0_wp
[2198]184                  ENDIF
185               ENDIF
186            END DO 
187         END DO 
188         !
[6140]189         ztmp = rday / ( rdt * REAL(nn_fsbc, wp) )
[2198]190         rscal(:,:) = rscal(:,:) * ztmp
[2715]191         !
[2198]192      ENDIF 
193         !     3. update qsr with the diurnal cycle
194         !     ------------------------------------
195
[3651]196      imask_night(:,:) = 0
[2198]197      DO jj = 1, jpj
198         DO ji = 1, jpi
[9124]199            ztmpm = 0._wp
[3764]200            IF( ABS(rab(ji,jj)) < 1. ) THEN         ! day duration is less than 24h
[2198]201               !
202               IF( rdawn(ji,jj) < rdusk(ji,jj) ) THEN       ! day time in one part
203                  zlousd = MAX(zlo, rdawn(ji,jj))
204                  zlousd = MIN(zlousd, zup)
205                  zupusd = MIN(zup, rdusk(ji,jj))
206                  zupusd = MAX(zupusd, zlo)
207                  ztmp = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
[2228]208                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
[3651]209                  ztmpm = zupusd - zlousd
210                  IF ( ztmpm .EQ. 0 ) imask_night(ji,jj) = 1
[2198]211                  !
212               ELSE                                         ! day time in two parts
213                  zlousd = MIN(zlo, rdusk(ji,jj))
214                  zupusd = MIN(zup, rdusk(ji,jj))
215                  ztmp1 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
[3651]216                  ztmpm1=zupusd-zlousd
[2198]217                  zlousd = MAX(zlo, rdawn(ji,jj))
218                  zupusd = MAX(zup, rdawn(ji,jj))
219                  ztmp2 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
[3651]220                  ztmpm2 =zupusd-zlousd
[2198]221                  ztmp = ztmp1 + ztmp2
[3651]222                  ztmpm = ztmpm1 + ztmpm2
[2228]223                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
[3651]224                  IF (ztmpm .EQ. 0.) imask_night(ji,jj) = 1
[2198]225               ENDIF
226            ELSE                                   ! 24h light or 24h night
227               !
[2228]228               IF( raa(ji,jj) > rbb(ji,jj) ) THEN           ! 24h day
[2198]229                  ztmp = fintegral(zlo, zup, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
[2228]230                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
[3651]231                  imask_night(ji,jj) = 0
[2198]232                  !
233               ELSE                                         ! No day
[3764]234                  zqsrout(ji,jj) = 0.0_wp
[3651]235                  imask_night(ji,jj) = 1
[2198]236               ENDIF
237            ENDIF
238         END DO 
239      END DO 
240      !
[9124]241      IF( PRESENT(l_mask) .AND. l_mask ) THEN
[3651]242         zqsrout(:,:) = float(imask_night(:,:))
243      ENDIF
244      !
[2228]245   END FUNCTION sbc_dcy
[2198]246
247   !!======================================================================
248END MODULE sbcdcy
Note: See TracBrowser for help on using the repository browser.