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

Last change on this file 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
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
[3294]20   USE timing           ! Timing
[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   !!----------------------------------------------------------------------
33   !! NEMO/OPA 3.3 , NEMO-consortium (2010)
[6486]34   !! $Id$
[2715]35   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[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            !
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
[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      !!----------------------------------------------------------------------
[3651]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
[2198]67      !!
[2228]68      INTEGER  ::   ji, jj                                       ! dummy loop indices
[3651]69      INTEGER, DIMENSION(jpi,jpj) :: imask_night ! night mask
[2198]70      REAL(wp) ::   ztwopi, zinvtwopi, zconvrad 
71      REAL(wp) ::   zlo, zup, zlousd, zupusd
[2228]72      REAL(wp) ::   zdsws, zdecrad, ztx, zsin, zcos
[2198]73      REAL(wp) ::   ztmp, ztmp1, ztmp2, ztest
[3651]74      REAL(wp) ::   ztmpm, ztmpm1, ztmpm2
[2228]75      REAL(wp), DIMENSION(jpi,jpj) ::   zqsrout                  ! output QSR flux with diurnal cycle
[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      !
83      IF( nn_timing == 1 )  CALL timing_start('sbc_dcy')
84      !
[2198]85      ! Initialization
86      ! --------------
[2715]87      ztwopi    = 2._wp * rpi
88      zinvtwopi = 1._wp / ztwopi
89      zconvrad  = ztwopi / 360._wp
[2198]90
91      ! When are we during the day (from 0 to 1)
[2715]92      zlo = ( REAL(nsec_day, wp) - 0.5_wp * rdttra(1) ) / rday
93      zup = zlo + ( REAL(nn_fsbc, wp)     * rdttra(1) ) / rday
[2198]94      !                                         
[3651]95      IF( nday_qsr == -1 ) THEN       ! first time step only 
[2198]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,*)
[11101]101            IF(lflush) CALL flush(numout)
[2198]102         ENDIF
[2715]103         ! allocate sbcdcy arrays
104         IF( sbc_dcy_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_dcy_alloc : unable to allocate arrays' )
[2198]105         ! Compute rcc needed to compute the time integral of the diurnal cycle
106         rcc(:,:) = zconvrad * glamt(:,:) - rpi
107         ! time of midday
[3764]108         rtmd(:,:) = 0.5_wp - glamt(:,:) / 360._wp
109         rtmd(:,:) = MOD( (rtmd(:,:) + 1._wp) , 1._wp)
[2198]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)         
[2228]122         zdsws = REAL(11 + nday_year, wp)
[2198]123         ! declination of the earths orbit
[3764]124         zdecrad = (-23.5_wp * zconvrad) * COS( zdsws * ztwopi / REAL(nyear_len(1),wp) )
[2198]125         ! Compute A and B needed to compute the time integral of the diurnal cycle
[3651]126
[2228]127         zsin = SIN( zdecrad )   ;   zcos = COS( zdecrad )
[2198]128         DO jj = 1, jpj
129            DO ji = 1, jpi
130               ztmp = zconvrad * gphit(ji,jj)
[2228]131               raa(ji,jj) = SIN( ztmp ) * zsin
132               rbb(ji,jj) = COS( ztmp ) * zcos
[2198]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
[3764]141               IF ( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h
[2198]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?
[3764]146                  IF ( ztest > 0._wp ) THEN
[2198]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
[3764]154                  rdawn(ji,jj) = rtmd(ji,jj) + 0.5_wp
[2198]155                  rdusk(ji,jj) = rdawn(ji,jj)
156               ENDIF
157             END DO 
158         END DO 
[2715]159         rdawn(:,:) = MOD( (rdawn(:,:) + 1._wp), 1._wp )
160         rdusk(:,:) = MOD( (rdusk(:,:) + 1._wp), 1._wp )
[3764]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)
[2198]165         DO jj = 1, jpj
166            DO ji = 1, jpi
[3764]167               IF ( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h
168                  rscal(ji,jj) = 0.0_wp
[2198]169                  IF ( rdawn(ji,jj) < rdusk(ji,jj) ) THEN      ! day time in one part
[3764]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
[2198]174                  ELSE                                         ! day time in two parts
[3764]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
[2198]180                  ENDIF
181               ELSE
182                  IF ( raa(ji,jj) > rbb(ji,jj) ) THEN         ! 24h day
[3764]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)
[2198]185                  ELSE                                          ! No day
[3764]186                     rscal(ji,jj) = 0.0_wp
[2198]187                  ENDIF
188               ENDIF
189            END DO 
190         END DO 
191         !
[2228]192         ztmp = rday / ( rdttra(1) * REAL(nn_fsbc, wp) )
[2198]193         rscal(:,:) = rscal(:,:) * ztmp
[2715]194         !
[2198]195      ENDIF 
196         !     3. update qsr with the diurnal cycle
197         !     ------------------------------------
198
[3651]199      imask_night(:,:) = 0
[2198]200      DO jj = 1, jpj
201         DO ji = 1, jpi
[3651]202            ztmpm = 0.0
[3764]203            IF( ABS(rab(ji,jj)) < 1. ) THEN         ! day duration is less than 24h
[2198]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)) 
[2228]211                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
[3651]212                  ztmpm = zupusd - zlousd
213                  IF ( ztmpm .EQ. 0 ) imask_night(ji,jj) = 1
[2198]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)) 
[3651]219                  ztmpm1=zupusd-zlousd
[2198]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)) 
[3651]223                  ztmpm2 =zupusd-zlousd
[2198]224                  ztmp = ztmp1 + ztmp2
[3651]225                  ztmpm = ztmpm1 + ztmpm2
[2228]226                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
[3651]227                  IF (ztmpm .EQ. 0.) imask_night(ji,jj) = 1
[2198]228               ENDIF
229            ELSE                                   ! 24h light or 24h night
230               !
[2228]231               IF( raa(ji,jj) > rbb(ji,jj) ) THEN           ! 24h day
[2198]232                  ztmp = fintegral(zlo, zup, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
[2228]233                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj)
[3651]234                  imask_night(ji,jj) = 0
[2198]235                  !
236               ELSE                                         ! No day
[3764]237                  zqsrout(ji,jj) = 0.0_wp
[3651]238                  imask_night(ji,jj) = 1
[2198]239               ENDIF
240            ENDIF
241         END DO 
242      END DO 
243      !
[3651]244      IF ( PRESENT(l_mask) .AND. l_mask ) THEN
245         zqsrout(:,:) = float(imask_night(:,:))
246      ENDIF
247      !
[3294]248      IF( nn_timing == 1 )  CALL timing_stop('sbc_dcy')
249      !
[2228]250   END FUNCTION sbc_dcy
[2198]251
252   !!======================================================================
253END MODULE sbcdcy
Note: See TracBrowser for help on using the repository browser.