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.
sbcblk_clio.F90 in trunk/NEMO/OPA_SRC/SBC – NEMO

source: trunk/NEMO/OPA_SRC/SBC/sbcblk_clio.F90 @ 892

Last change on this file since 892 was 892, checked in by rblod, 16 years ago

Correct light bugs after the introduction of the new surface module, see ticket #113

File size: 53.2 KB
Line 
1MODULE sbcblk_clio
2   !!======================================================================
3   !!                   ***  MODULE  sbcblk_clio  ***
4   !! Ocean forcing:  bulk thermohaline forcing of the ocean (or ice)
5   !!=====================================================================
6   !! History :  OPA  !  1997-06 (Louvain-La-Neuve)  Original code
7   !!                 !  2001-04 (C. Ethe) add flx_blk_declin
8   !!   NEMO     2.0  !  2002-08 (C. Ethe, G. Madec) F90: Free form and module
9   !!            3.0  !  2008-03 (C. Talandier, G. Madec) surface module + LIM3
10   !!----------------------------------------------------------------------
11   !!   sbc_blk_clio   : CLIO bulk formulation: read and update required input fields
12   !!   blk_clio_oce   : ocean CLIO bulk formulea: compute momentum, heat and freswater fluxes for the ocean
13   !!   blk_ice_clio   : ice   CLIO bulk formulea: compute momentum, heat and freswater fluxes for the sea-ice
14   !!   blk_clio_qsr_oce : shortwave radiation for ocean computed from the cloud cover
15   !!   blk_clio_qsr_ice : shortwave radiation for ice   computed from the cloud cover
16   !!   flx_blk_declin : solar declinaison
17   !!----------------------------------------------------------------------
18   USE oce             ! ocean dynamics and tracers
19   USE dom_oce         ! ocean space and time domain
20   USE phycst          ! physical constants
21   USE daymod          ! calendar
22   USE ocfzpt          ! ocean freezing point
23   USE fldread         ! read input fields
24   USE sbc_oce         ! Surface boundary condition: ocean fields
25   USE iom             ! I/O manager library
26   USE in_out_manager  ! I/O manager
27   USE lib_mpp         ! distribued memory computing library
28   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
29
30   USE albedo
31   USE prtctl          ! Print control
32#if defined key_lim3
33   USE par_ice
34   USE ice
35   USE ice_oce         ! For ice surface temperature
36#elif defined key_lim2
37   USE par_ice_2
38   USE ice_2
39#endif
40   IMPLICIT NONE
41   PRIVATE
42
43   PUBLIC sbc_blk_clio        ! routine called by sbcmod.F90
44   PUBLIC blk_ice_clio        ! routine called by sbcice_lim.F90
45
46   INTEGER , PARAMETER ::   jpfld   = 7           ! maximum number of files to read
47   INTEGER , PARAMETER ::   jp_utau = 1           ! index of wind stress (i-component)      (N/m2)    at U-point
48   INTEGER , PARAMETER ::   jp_vtau = 2           ! index of wind stress (j-component)      (N/m2)    at V-point
49   INTEGER , PARAMETER ::   jp_wndm = 3           ! index of 10m wind module                 (m/s)    at T-point
50   INTEGER , PARAMETER ::   jp_humi = 4           ! index of specific humidity               ( % )
51   INTEGER , PARAMETER ::   jp_ccov = 5           ! index of cloud cover                     ( % )
52   INTEGER , PARAMETER ::   jp_tair = 6           ! index of 10m air temperature             (Kelvin)
53   INTEGER , PARAMETER ::   jp_prec = 7           ! index of total precipitation (rain+snow) (Kg/m2/s)
54   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read)
55
56   INTEGER, PARAMETER  ::   jpintsr = 24          ! number of time step between sunrise and sunset
57   !                                              ! uses for heat flux computation
58   LOGICAL ::   lbulk_init = .TRUE.               ! flag, bulk initialization done or not)
59
60   REAL(wp) ::   cai = 1.40e-3 ! best estimate of atm drag in order to get correct FS export in ORCA2-LIM
61   REAL(wp) ::   cao = 1.00e-3 ! chosen by default  ==> should depends on many things...  !!gmto be updated
62
63   REAL(wp) ::   yearday     !: number of days per year   
64   REAL(wp) ::   rdtbs2      !: number of days per year   
65   
66   REAL(wp), DIMENSION(19)  ::  budyko            ! BUDYKO's coefficient (cloudiness effect on LW radiation)
67   DATA budyko / 1.00, 0.98, 0.95, 0.92, 0.89, 0.86, 0.83, 0.80, 0.78, 0.75,   &
68      &          0.72, 0.69, 0.67, 0.64, 0.61, 0.58, 0.56, 0.53, 0.50 /
69   REAL(wp), DIMENSION(20)  :: tauco              ! cloud optical depth coefficient
70   DATA tauco / 6.6, 6.6, 7.0, 7.2, 7.1, 6.8, 6.5, 6.6, 7.1, 7.6,   &
71      &         6.6, 6.1, 5.6, 5.5, 5.8, 5.8, 5.6, 5.6, 5.6, 5.6 /
72   !!
73   REAL(wp), DIMENSION(jpi,jpj) ::   sbudyko      ! cloudiness effect on LW radiation
74   REAL(wp), DIMENSION(jpi,jpj) ::   stauc        ! cloud optical depth
75   
76   REAL(wp)  ::   zeps    = 1.e-20                ! constant values
77   REAL(wp)  ::   zeps0   = 1.e-13 
78
79#  include "vectopt_loop_substitute.h90"
80   !!----------------------------------------------------------------------
81   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)
82   !! $Id:$
83   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
84   !!----------------------------------------------------------------------
85
86CONTAINS
87
88
89   SUBROUTINE sbc_blk_clio( kt )
90      !!---------------------------------------------------------------------
91      !!                    ***  ROUTINE sbc_blk_clio  ***
92      !!                   
93      !! ** Purpose :   provide at each time step the surface ocean fluxes
94      !!      (momentum, heat, freshwater and runoff)
95      !!
96      !! ** Method  :   READ each fluxes in NetCDF files
97      !!      The i-component of the stress                utau   (N/m2)
98      !!      The j-component of the stress                vtau   (N/m2)
99      !!      the net downward heat flux                   qtot   (watt/m2)
100      !!      the net downward radiative flux              qsr    (watt/m2)
101      !!      the net upward water (evapo - precip)        emp    (kg/m2/s)
102      !!                Assumptions made:
103      !!       - each file content an entire year (read record, not the time axis)
104      !!       - first and last record are part of the previous and next year
105      !!         (useful for time interpolation)
106      !!       - the number of records is 2 + 365*24 / freqh(jf)
107      !!         or 366 in leap year case
108      !!
109      !!      C A U T I O N : never mask the surface stress fields
110      !!                      the stress is assumed to be in the mesh referential
111      !!                      i.e. the (i,j) referential
112      !!
113      !! ** Action  :   defined at each time-step at the air-sea interface
114      !!              - utau  &  vtau   : stress components in geographical ref.
115      !!              - qns   &  qsr    : non solar and solar heat fluxes
116      !!              - emp             : evap - precip (volume flux)
117      !!              - emps            : evap - precip (concentration/dillution)
118      !!----------------------------------------------------------------------
119      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
120      !!
121      INTEGER  ::   jf, ifpr, jfpr     ! dummy indices
122      INTEGER  ::   ierror             ! return error code
123      !!
124      CHARACTER(len=100) ::  cn_dir                            !   Root directory for location of CLIO files
125      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read
126      TYPE(FLD_N) ::   sn_utau, sn_vtau, sn_wndm, sn_tair      ! informations about the fields to be read
127      TYPE(FLD_N) ::   sn_humi, sn_ccov, sn_prec               !   "                                 "
128      !!
129      NAMELIST/namsbc_clio/ cn_dir, sn_utau, sn_vtau, sn_wndm, sn_humi,   &
130         &                          sn_ccov, sn_tair, sn_prec
131      !!---------------------------------------------------------------------
132
133      !                                         ! ====================== !
134      IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
135         !                                      ! ====================== !
136         ! set file information (default values)
137         cn_dir = './'       ! directory in which the model is executed
138
139         ! (NB: frequency positive => hours, negative => months)
140         !            !    file     ! frequency !  variable  ! time intep !  clim  ! starting !
141         !            !    name     !  (hours)  !   name     !   (T/F)    !  (0/1) !  record  !
142         sn_utau = FLD_N( 'utau'    ,    24.    ,  'utau'    ,  .true.    ,    0   ,     0    ) 
143         sn_vtau = FLD_N( 'vtau'    ,    24.    ,  'vtau'    ,  .true.    ,    0   ,     0    ) 
144         sn_wndm = FLD_N( 'mwnd10m' ,    24.    ,  'm_10'    ,  .true.    ,    0   ,     0    ) 
145         sn_tair = FLD_N( 'tair10m' ,    24.    ,  't_10'    ,  .FALSE.   ,    0   ,     0    ) 
146         sn_humi = FLD_N( 'humi10m' ,    24.    ,  'q_10'    ,  .FALSE.   ,    0   ,     0    ) 
147         sn_ccov = FLD_N( 'ccover'  ,   -12.    ,  'cloud'   ,  .TRUE.    ,    0   ,     0    ) 
148         sn_prec = FLD_N( 'precip'  ,   -12.    ,  'precip'  ,  .TRUE.    ,    0   ,     0    ) 
149
150         REWIND( numnam )                    ! ... read in namlist namsbc_clio
151         READ  ( numnam, namsbc_clio )
152
153         ! store namelist information in an array
154         slf_i(jp_utau) = sn_utau   ;   slf_i(jp_vtau) = sn_vtau   ;   slf_i(jp_wndm) = sn_wndm
155         slf_i(jp_tair) = sn_tair   ;   slf_i(jp_humi) = sn_humi
156         slf_i(jp_ccov) = sn_ccov   ;   slf_i(jp_prec) = sn_prec
157         
158         ! set sf structure
159         ALLOCATE( sf(jpfld), STAT=ierror )
160         IF( ierror > 0 ) THEN
161            CALL ctl_stop( 'sbc_blk_clio: unable to allocate sf structure' )   ;   RETURN
162         ENDIF
163         !
164         DO jf = 1, jpfld
165            WRITE(sf(jf)%clrootname,'(a,a)' )   TRIM( cn_dir ), TRIM( slf_i(jf)%clname )
166            sf(jf)%freqh   = slf_i(jf)%freqh
167            sf(jf)%clvar   = slf_i(jf)%clvar
168            sf(jf)%ln_tint = slf_i(jf)%ln_tint
169            sf(jf)%nclim   = slf_i(jf)%nclim
170            sf(jf)%nstrec  = slf_i(jf)%nstrec
171         END DO
172
173         IF(lwp) THEN      ! control print
174            WRITE(numout,*)           
175            WRITE(numout,*) 'sbc_blk_clio : flux formulattion for ocean surface boundary condition'
176            WRITE(numout,*) '~~~~~~~~~~~~ '
177            WRITE(numout,*) '          namsbc_clio Namelist'
178            WRITE(numout,*) '          list of files and frequency (>0: in hours ; <0 in months)'
179            DO jf = 1, jpfld
180                WRITE(numout,*) '               file root name: ' , TRIM( sf(jf)%clrootname ),   &
181                   &                          ' variable name: '  , TRIM( sf(jf)%clvar      )
182                WRITE(numout,*) '               frequency: '      ,       sf(jf)%freqh       ,   &
183                   &                          ' time interp: '    ,       sf(jf)%ln_tint     ,   &
184                   &                          ' climatology: '    ,       sf(jf)%nclim       ,   &
185                   &                          ' starting record: ',       sf(jf)%nstrec
186            END DO
187         ENDIF
188         !
189      ENDIF
190      !                                         ! ====================== !
191      !                                         !    At each time-step   !
192      !                                         ! ====================== !
193      !
194      CALL fld_read( kt, nn_fsbc, sf )                ! input fields provided at the current time-step
195      !
196#if defined key_lim3     
197      tatm_ice(:,:) = sf(jp_tair)%fnow(:,:)     !RB ugly patch
198#endif
199      !
200      IF(lwp .AND. nitend-nit000 <= 100 ) THEN
201         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
202            WRITE(numout,*)
203            WRITE(numout,*) ' read monthly CLIO fluxes: ok, kt: ', kt
204            WRITE(numout,*)
205            ifpr = INT(jpi/8)      ;      jfpr = INT(jpj/10)
206            WRITE(numout,*) TRIM(sf(jp_utau)%clvar),' day: ',ndastp
207            CALL prihre( sf(jp_utau)%fnow,jpi,jpj,1,jpi,ifpr,1,jpj,jfpr,0.,numout )
208            WRITE(numout,*)
209            WRITE(numout,*) TRIM(sf(jp_vtau)%clvar),' day: ',ndastp
210            CALL prihre( sf(jp_vtau)%fnow,jpi,jpj,1,jpi,ifpr,1,jpj,jfpr,0.,numout )
211            WRITE(numout,*)
212            WRITE(numout,*) TRIM(sf(jp_humi)%clvar),' day: ',ndastp
213            CALL prihre( sf(jp_humi)%fnow,jpi,jpj,1,jpi,ifpr,1,jpj,jfpr,0.,numout )
214            WRITE(numout,*)
215            WRITE(numout,*) TRIM(sf(jp_wndm)%clvar),' day: ',ndastp
216            CALL prihre( sf(jp_wndm)%fnow,jpi,jpj,1,jpi,ifpr,1,jpj,jfpr,0.,numout )
217            WRITE(numout,*)
218            WRITE(numout,*) TRIM(sf(jp_ccov)%clvar),' day: ',ndastp
219            CALL prihre( sf(jp_ccov)%fnow,jpi,jpj,1,jpi,ifpr,1,jpj,jfpr,0.,numout )
220            WRITE(numout,*)
221            WRITE(numout,*) TRIM(sf(jp_prec)%clvar),' day: ',ndastp
222            CALL prihre( sf(jp_prec)%fnow,jpi,jpj,1,jpi,ifpr,1,jpj,jfpr,0.,numout )
223            WRITE(numout,*)
224            WRITE(numout,*) TRIM(sf(jp_tair)%clvar),' day: ',ndastp
225            CALL prihre( sf(jp_tair)%fnow,jpi,jpj,1,jpi,ifpr,1,jpj,jfpr,0.,numout )
226            WRITE(numout,*)
227         ENDIF
228      ENDIF
229
230      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
231          CALL blk_oce_clio( sst_m, ssu_m, ssv_m )    ! compute the surface ocean fluxes using CLIO bulk formulea
232      ENDIF                                           !
233      !
234   END SUBROUTINE sbc_blk_clio
235
236
237   SUBROUTINE blk_oce_clio( pst, pu, pv )
238      !!---------------------------------------------------------------------------
239      !!                     ***  ROUTINE blk_oce_clio  ***
240      !!                 
241      !!  ** Purpose :   Compute momentum, heat and freshwater fluxes at ocean surface
242      !!               using CLIO bulk formulea
243      !!         
244      !!  ** Method  :   The flux of heat at the ocean surfaces are derived
245      !!       from semi-empirical ( or bulk ) formulae which relate the flux to
246      !!       the properties of the surface and of the lower atmosphere. Here, we
247      !!       follow the work of Oberhuber, 1988   
248      !!               - momentum flux (stresses) directly read in files at U- and V-points
249      !!               - compute ocean/ice albedos (call albedo_oce/albedo_ice) 
250      !!               - compute shortwave radiation for ocean (call blk_clio_qsr_oce)
251      !!               - compute long-wave radiation for the ocean
252      !!               - compute the turbulent heat fluxes over the ocean
253      !!               - deduce the evaporation over the ocean
254      !!  ** Action  :   Fluxes over the ocean:
255      !!               - utau, vtau  i- and j-component of the wind stress
256      !!               - qns, qsr    non-slor and solar heat flux
257      !!               - emp, emps   evaporation minus precipitation
258      !!----------------------------------------------------------------------
259      REAL(wp), INTENT(in), DIMENSION(jpi,jpj) ::   pst   ! surface temperature                      [Celcius]
260      REAL(wp), INTENT(in), DIMENSION(jpi,jpj) ::   pu    ! surface current at U-point (i-component) [m/s]
261      REAL(wp), INTENT(in), DIMENSION(jpi,jpj) ::   pv    ! surface current at V-point (j-component) [m/s]
262      !!
263      INTEGER  ::   ji, jj   ! dummy loop indices
264      !!
265      REAL(wp) ::   zrhova, zcsho, zcleo, zcldeff               ! temporary scalars
266      REAL(wp) ::   zqsato, zdteta, zdeltaq, ztvmoy, zobouks    !    -         -
267      REAL(wp) ::   zpsims, zpsihs, zpsils, zobouku, zxins, zpsimu   !    -         -
268      REAL(wp) ::   zpsihu, zpsilu, zstab,zpsim, zpsih, zpsil   !    -         -
269      REAL(wp) ::   zvatmg, zcmn, zchn, zcln, zcmcmn, zdenum    !    -         -
270      REAL(wp) ::   zdtetar, ztvmoyr, zlxins, zchcm, zclcm      !    -         -
271      REAL(wp) ::   zmt1, zmt2, zmt3, ztatm3, ztamr, ztaevbk    !    -         -
272      REAL(wp) ::   zsst, ztatm, zcco1, zpatm, zinda            !    -         -
273      REAL(wp) ::   zrhoa, zev, zes, zeso, zqatm, zevsqr        !    -         -
274      !!
275      REAL(wp), DIMENSION(jpi,jpj) ::   zqlw        ! long-wave heat flux over ocean
276      REAL(wp), DIMENSION(jpi,jpj) ::   zqla        ! latent heat flux over ocean
277      REAL(wp), DIMENSION(jpi,jpj) ::   zqsb        ! sensible heat flux over ocean
278      !!---------------------------------------------------------------------
279
280      zpatm = 101000.      ! atmospheric pressure  (assumed constant here)
281
282      !------------------------------------!
283      !   momentum fluxes  (utau, vtau )   !
284      !------------------------------------!
285!CDIR COLLAPSE
286      DO jj = 1 , jpj
287         DO ji = 1, jpi
288            utau(ji,jj) = sf(jp_utau)%fnow(ji,jj)
289            vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj)
290         END DO
291      END DO
292
293      !------------------------------------------------!
294      !   Shortwave radiation for ocean and snow/ice   !
295      !------------------------------------------------!
296     
297      CALL blk_clio_qsr_oce( qsr )
298
299      ! CAUTION: ocean shortwave radiation sets to zero if more than 50% of sea-ice !!gm to be removed
300      DO jj = 1, jpj
301         DO ji = 1, jpi
302            zinda    = MAX(  0.e0, SIGN(  1.e0, -( -1.5 - freeze(ji,jj) )  )  )
303            qsr(ji,jj) = zinda * qsr(ji,jj)
304         END DO
305      END DO
306
307
308      !------------------------!
309      !   Other ocean fluxes   !
310      !------------------------!
311!CDIR NOVERRCHK
312!CDIR COLLAPSE
313      DO jj = 1, jpj
314!CDIR NOVERRCHK
315         DO ji = 1, jpi
316            !
317            zsst  = pst(ji,jj)              + rt0           ! converte Celcius to Kelvin the SST and Tair
318            ztatm = sf(jp_tair)%fnow(ji,jj) + rt0           ! and set minimum value far above 0 K (=rt0 over land)
319            zcco1 = 1.0 - sf(jp_ccov)%fnow(ji,jj)           ! fraction of clear sky ( 1 - cloud cover)
320            zrhoa = zpatm / ( 287.04 * ztatm )              ! air density (equation of state for dry air)
321            ztamr = ztatm - rtt                             ! Saturation water vapour
322            zmt1  = SIGN( 17.269,  ztamr )                  !           ||
323            zmt2  = SIGN( 21.875,  ztamr )                  !          \  /
324            zmt3  = SIGN( 28.200, -ztamr )                  !           \/
325            zes   = 611.0 * EXP(  ABS( ztamr ) * MIN ( zmt1, zmt2 ) / ( ztatm - 35.86  + MAX( 0.e0, zmt3 ) )  )
326            zev    = sf(jp_humi)%fnow(ji,jj) * zes          ! vapour pressure 
327            zevsqr = SQRT( zev * 0.01 )                     ! square-root of vapour pressure
328            zqatm = 0.622 * zev / ( zpatm - 0.378 * zev )   ! specific humidity
329
330            !--------------------------------------!
331            !  long-wave radiation over the ocean  !  ( Berliand 1952 ; all latitudes )
332            !--------------------------------------!
333            ztatm3  = ztatm * ztatm * ztatm
334            zcldeff = 1.0 - sbudyko(ji,jj) * sf(jp_ccov)%fnow(ji,jj) * sf(jp_ccov)%fnow(ji,jj)   
335            ztaevbk = ztatm * ztatm3 * zcldeff * ( 0.39 - 0.05 * zevsqr ) 
336            !
337            zqlw(ji,jj) = - emic * stefan * ( ztaevbk + 4. * ztatm3 * ( zsst - ztatm ) ) 
338
339            !--------------------------------------------------
340            !  Latent and sensible heat fluxes over the ocean
341            !--------------------------------------------------
342            !                                                          ! vapour pressure at saturation of ocean
343            zeso =  611.0 * EXP ( 17.2693884 * ( zsst - rtt ) * tmask(ji,jj,1) / ( zsst - 35.86 ) )
344
345            zqsato = ( 0.622 * zeso ) / ( zpatm - 0.378 * zeso )       ! humidity close to the ocean surface (at saturation)
346
347            ! Drag coefficients from Large and Pond (1981,1982)
348            !                                                          ! Stability parameters
349            zdteta  = zsst - ztatm
350            zdeltaq = zqatm - zqsato
351            ztvmoy  = ztatm * ( 1. + 2.2e-3 * ztatm * zqatm )
352            zdenum  = MAX( sf(jp_wndm)%fnow(ji,jj) * sf(jp_wndm)%fnow(ji,jj) * ztvmoy, zeps )
353            zdtetar = zdteta / zdenum
354            ztvmoyr = ztvmoy * ztvmoy * zdeltaq / zdenum
355            !                                                          ! case of stable atmospheric conditions
356            zobouks = -70.0 * 10. * ( zdtetar + 3.2e-3 * ztvmoyr )
357            zobouks = MAX( 0.e0, zobouks )
358            zpsims = -7.0 * zobouks
359            zpsihs =  zpsims
360            zpsils =  zpsims
361            !                                                          ! case of unstable atmospheric conditions
362            zobouku = MIN(  0.e0, -100.0 * 10.0 * ( zdtetar + 2.2e-3 * ztvmoyr )  )
363            zxins   = ( 1. - 16. * zobouku )**0.25
364            zlxins  = LOG( ( 1. + zxins * zxins ) / 2. )
365            zpsimu  = 2. * LOG( ( 1 + zxins ) * 0.5 )  + zlxins - 2. * ATAN( zxins ) + rpi * 0.5
366            zpsihu  = 2. * zlxins
367            zpsilu  = zpsihu
368            !                                                          ! intermediate values
369            zstab   = MAX( 0.e0, SIGN( 1.e0, zdteta ) )
370            zpsim   = zstab * zpsimu + ( 1.0 - zstab ) * zpsims
371            zpsih   = zstab * zpsihu + ( 1.0 - zstab ) * zpsihs
372            zpsil   = zpsih
373           
374            zvatmg         = MAX( 0.032 * 1.5e-3 * sf(jp_wndm)%fnow(ji,jj) * sf(jp_wndm)%fnow(ji,jj) / grav, zeps )
375            zcmn           = vkarmn / LOG ( 10. / zvatmg )
376            zchn           = 0.0327 * zcmn
377            zcln           = 0.0346 * zcmn
378            zcmcmn         = 1. / ( 1. - zcmn * zpsim / vkarmn )
379            zchcm          = zcmcmn / ( 1. - zchn * zpsih / ( vkarmn * zcmn ) )
380            zclcm          = zchcm
381            !                                                          ! transfert coef. (Large and Pond 1981,1982)
382            zcsho          = zchn * zchcm                               
383            zcleo          = zcln * zclcm 
384
385            zrhova         = zrhoa * sf(jp_wndm)%fnow(ji,jj)
386
387            ! sensible heat flux
388            zqsb(ji,jj) = zrhova * zcsho * 1004.0  * ( zsst - ztatm ) 
389         
390            ! latent heat flux (bounded by zero)
391            zqla(ji,jj) = MAX(  0.e0, zrhova * zcleo * 2.5e+06 * ( zqsato - zqatm )  )
392            !               
393         END DO
394      END DO
395     
396      ! ----------------------------------------------------------------------------- !
397      !     III    Total FLUXES                                                       !
398      ! ----------------------------------------------------------------------------- !
399!CDIR COLLAPSE
400      qns (:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)      ! Downward Non Solar flux
401!CDIR COLLAPSE
402      emp (:,:) = zqla(:,:) / cevap - sf(jp_prec)%fnow(:,:) / rday * tmask(:,:,1)
403!CDIR COLLAPSE
404      emps(:,:) = emp(:,:)
405      !
406      IF(ln_ctl) THEN
407         CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce_clio: zqsb   : ', tab2d_2=zqlw , clinfo2=' zqlw  : ')
408         CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce_clio: zqla   : ', tab2d_2=qsr  , clinfo2=' qsr   : ')
409         CALL prt_ctl(tab2d_1=pst  , clinfo1=' blk_oce_clio: pst    : ', tab2d_2=emp  , clinfo2=' emp   : ')
410         CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce_clio: utau   : ', mask1=umask,   &
411            &         tab2d_2=vtau , clinfo2=' vtau : ', mask2=vmask )
412      ENDIF
413
414   END SUBROUTINE blk_oce_clio
415
416
417   SUBROUTINE blk_ice_clio(  pst   , pui   , pvi   , palb_cs, palb_os ,   &
418      &                      p_taui, p_tauj, p_qns , p_qsr,   &
419      &                      p_qla , p_dqns, p_dqla,          &
420      &                      p_tpr , p_spr ,                  &
421      &                      p_fr1 , p_fr2 , cd_grid  )
422      !!---------------------------------------------------------------------------
423      !!                     ***  ROUTINE blk_ice_clio  ***
424      !!                 
425      !!  ** Purpose :   Computation of the heat fluxes at ocean and snow/ice
426      !!       surface the solar heat at ocean and snow/ice surfaces and the
427      !!       sensitivity of total heat fluxes to the SST variations
428      !!         
429      !!  ** Method  :   The flux of heat at the ice and ocean surfaces are derived
430      !!       from semi-empirical ( or bulk ) formulae which relate the flux to
431      !!       the properties of the surface and of the lower atmosphere. Here, we
432      !!       follow the work of Oberhuber, 1988   
433      !!
434      !!  ** Action  :   call albedo_oce/albedo_ice to compute ocean/ice albedo
435      !!          computation of snow precipitation
436      !!          computation of solar flux at the ocean and ice surfaces
437      !!          computation of the long-wave radiation for the ocean and sea/ice
438      !!          computation of turbulent heat fluxes over water and ice
439      !!          computation of evaporation over water
440      !!          computation of total heat fluxes sensitivity over ice (dQ/dT)
441      !!          computation of latent heat flux sensitivity over ice (dQla/dT)
442      !!
443      !!----------------------------------------------------------------------
444      REAL(wp), INTENT(in   ), DIMENSION(:,:,:)   ::   pst      ! ice surface temperature                   [Kelvin]
445      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   pui      ! ice surface velocity (i-component, I-point)  [m/s]
446      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   pvi      ! ice surface velocity (j-component, I-point)  [m/s]
447      REAL(wp), INTENT(in   ), DIMENSION(:,:,:)   ::   palb_cs  ! ice albedo (clear    sky) (alb_ice_cs)         [%]
448      REAL(wp), INTENT(in   ), DIMENSION(:,:,:)   ::   palb_os  ! ice albedo (overcast sky) (alb_ice_os)         [%]
449      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_taui   ! surface ice stress at I-point (i-component) [N/m2]
450      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_tauj   ! surface ice stress at I-point (j-component) [N/m2]
451      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_qns    ! non solar heat flux over ice (T-point)      [W/m2]
452      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_qsr    !     solar heat flux over ice (T-point)      [W/m2]
453      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_qla    ! latent    heat flux over ice (T-point)      [W/m2]
454      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_dqns   ! non solar heat sensistivity  (T-point)      [W/m2]
455      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_dqla   ! latent    heat sensistivity  (T-point)      [W/m2]
456      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_tpr    ! total precipitation          (T-point)   [Kg/m2/s]
457      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_spr    ! solid precipitation          (T-point)   [Kg/m2/s]
458      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_fr1    ! 1sr fraction of qsr penetration in ice         [%]
459      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_fr2    ! 2nd fraction of qsr penetration in ice         [%]
460      CHARACTER(len=1), INTENT(in   )             ::   cd_grid  ! type of sea-ice grid ("C" or "B" grid)
461      !!
462      INTEGER  ::   ji, jj, jl    ! dummy loop indices
463      INTEGER  ::   ijpl          ! number of ice categories (size of 3rd dim of input arrays)
464      !!
465      REAL(wp) ::   zcoef, zmt1, zmt2, zmt3, ztatm3             ! temporary scalars
466      REAL(wp) ::   ztaevbk, zind1, zind2, zind3, ztamr         !    -         -
467      REAL(wp) ::   zesi, zqsati, zdesidt                       !    -         -
468      REAL(wp) ::   zdqla, zcldeff, zev, zes, zpatm, zrhova     !    -         -
469      REAL(wp) ::   zcshi, zclei, zrhovaclei, zrhovacshi        !    -         -
470      REAL(wp) ::   ztice3, zticemb, zticemb2, zdqlw, zdqsb     !    -         -
471      !!
472      REAL(wp), DIMENSION(jpi,jpj) ::   ztatm   ! Tair in Kelvin
473      REAL(wp), DIMENSION(jpi,jpj) ::   zqatm   ! specific humidity
474      REAL(wp), DIMENSION(jpi,jpj) ::   zevsqr  ! vapour pressure square-root
475      REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa   ! air density
476      REAL(wp), DIMENSION(jpi,jpj,SIZE(pst,3)) ::   z_qlw, z_qsb
477      !!---------------------------------------------------------------------
478
479      ijpl  = SIZE( pst, 3 )                 ! number of ice categories
480      zpatm = 101000.                        ! atmospheric pressure  (assumed constant  here)
481
482      !------------------------------------!
483      !   momentum fluxes  (utau, vtau )   !
484      !------------------------------------!
485
486      SELECT CASE( cd_grid )
487      CASE( 'C' )                          ! C-grid ice dynamics
488         ! Change from wind speed to wind stress over OCEAN (cao is used)
489         zcoef = cai / cao 
490!CDIR COLLAPSE
491         DO jj = 1 , jpj
492            DO ji = 1, jpi
493               p_taui(ji,jj) = zcoef * utau(ji,jj)
494               p_tauj(ji,jj) = zcoef * vtau(ji,jj)
495            END DO
496         END DO
497      CASE( 'B' )                          ! B-grid ice dynamics
498         ! stress from ocean U- and V-points to ice U,V point
499!CDIR COLLAPSE
500         DO jj = 2, jpj
501            DO ji = fs_2, jpi   ! vector opt.
502               p_taui(ji,jj) = 0.5 * ( utau(ji-1,jj  ) + utau(ji-1,jj-1) )
503               p_tauj(ji,jj) = 0.5 * ( vtau(ji  ,jj-1) + vtau(ji-1,jj-1) )
504            END DO
505         END DO
506         CALL lbc_lnk( p_taui(:,:), 'I', -1. )   ! I-point (i.e. ice U-V point)
507         CALL lbc_lnk( p_tauj(:,:), 'I', -1. )   ! I-point (i.e. ice U-V point)
508      END SELECT
509
510
511      !  Determine cloud optical depths as a function of latitude (Chou et al., 1981).
512      !  and the correction factor for taking into account  the effect of clouds
513      !------------------------------------------------------
514!CDIR NOVERRCHK
515!CDIR COLLAPSE
516      DO jj = 1, jpj
517!CDIR NOVERRCHK
518         DO ji = 1, jpi
519            ztatm (ji,jj) = sf(jp_tair)%fnow(ji,jj) + rt0            ! air temperature in Kelvins
520     
521            zrhoa(ji,jj) = zpatm / ( 287.04 * ztatm(ji,jj) )         ! air density (equation of state for dry air)
522     
523            ztamr = ztatm(ji,jj) - rtt                               ! Saturation water vapour
524            zmt1  = SIGN( 17.269,  ztamr )
525            zmt2  = SIGN( 21.875,  ztamr )
526            zmt3  = SIGN( 28.200, -ztamr )
527            zes   = 611.0 * EXP(  ABS( ztamr ) * MIN ( zmt1, zmt2 )   &
528               &                / ( ztatm(ji,jj) - 35.86  + MAX( 0.e0, zmt3 ) )  )
529
530            zev = sf(jp_humi)%fnow(ji,jj) * zes                      ! vapour pressure 
531            zevsqr(ji,jj) = SQRT( zev * 0.01 )                       ! square-root of vapour pressure
532            zqatm(ji,jj) = 0.622 * zev / ( zpatm - 0.378 * zev )     ! specific humidity
533
534            !----------------------------------------------------
535            !   Computation of snow precipitation (Ledley, 1985) |
536            !----------------------------------------------------
537            zmt1  =   253.0 - ztatm(ji,jj)            ;   zind1 = MAX( 0.e0, SIGN( 1.e0, zmt1 ) )
538            zmt2  = ( 272.0 - ztatm(ji,jj) ) / 38.0   ;   zind2 = MAX( 0.e0, SIGN( 1.e0, zmt2 ) )
539            zmt3  = ( 281.0 - ztatm(ji,jj) ) / 18.0   ;   zind3 = MAX( 0.e0, SIGN( 1.e0, zmt3 ) )
540            p_spr(ji,jj) = sf(jp_prec)%fnow(ji,jj) / rday   &        ! rday = converte mm/day to kg/m2/s
541               &         * (          zind1      &                   ! solid  (snow) precipitation [kg/m2/s]
542               &            + ( 1.0 - zind1 ) * (          zind2   * ( 0.5 + zmt2 )   &
543               &                                 + ( 1.0 - zind2 ) *  zind3 * zmt3  )   ) 
544
545            !----------------------------------------------------!
546            !  fraction of net penetrative shortwave radiation   !
547            !----------------------------------------------------!
548            ! fraction of qsr_ice which is NOT absorbed in the thin surface layer
549            ! and thus which penetrates inside the ice cover ( Maykut and Untersteiner, 1971 ; Elbert anbd Curry, 1993 )
550            p_fr1(ji,jj) = 0.18  * ( 1.e0 - sf(jp_ccov)%fnow(ji,jj) ) + 0.35 * sf(jp_ccov)%fnow(ji,jj) 
551            p_fr2(ji,jj) = 0.82  * ( 1.e0 - sf(jp_ccov)%fnow(ji,jj) ) + 0.65 * sf(jp_ccov)%fnow(ji,jj)
552         END DO
553      END DO
554     
555      !-----------------------------------------------------------!
556      !  snow/ice Shortwave radiation   (abedo already computed)  !
557      !-----------------------------------------------------------!
558      CALL blk_clio_qsr_ice( palb_cs, palb_os, p_qsr )
559
560      !                                     ! ========================== !
561      DO jl = 1, ijpl                       !  Loop over ice categories  !
562         !                                  ! ========================== !
563!CDIR NOVERRCHK
564!CDIR COLLAPSE
565         DO jj = 1 , jpj
566!CDIR NOVERRCHK
567            DO ji = 1, jpi
568               !-------------------------------------------!
569               !  long-wave radiation over ice categories  !  ( Berliand 1952 ; all latitudes )
570               !-------------------------------------------!
571               ztatm3  = ztatm(ji,jj) * ztatm(ji,jj) * ztatm(ji,jj)
572               zcldeff = 1.0 - sbudyko(ji,jj) * sf(jp_ccov)%fnow(ji,jj) * sf(jp_ccov)%fnow(ji,jj)   
573               ztaevbk = ztatm3 * ztatm(ji,jj) * zcldeff * ( 0.39 - 0.05 * zevsqr(ji,jj) ) 
574               !
575               z_qlw(ji,jj,jl) = - emic * stefan * ( ztaevbk + 4. * ztatm3 * ( pst(ji,jj,jl) - ztatm(ji,jj) ) ) 
576
577               !----------------------------------------
578               !  Turbulent heat fluxes over snow/ice     ( Latent and sensible )
579               !----------------------------------------       
580
581               ! vapour pressure at saturation of ice (tmask to avoid overflow in the exponential)
582               zesi =  611.0 * EXP( 21.8745587 * tmask(ji,jj,1) * ( pst(ji,jj,jl) - rtt )/ ( pst(ji,jj,jl) - 7.66 ) )
583               ! humidity close to the ice surface (at saturation)
584               zqsati   = ( 0.622 * zesi ) / ( zpatm - 0.378 * zesi )
585               
586               !  computation of intermediate values
587               zticemb  = pst(ji,jj,jl) - 7.66
588               zticemb2 = zticemb * zticemb 
589               ztice3   = pst(ji,jj,jl) * pst(ji,jj,jl) * pst(ji,jj,jl)
590               zdesidt  = zesi * ( 9.5 * LOG( 10.0 ) * ( rtt - 7.66 )  / zticemb2 )
591               
592               !  Transfer cofficients assumed to be constant (Parkinson 1979 ; Maykut 1982)
593               zcshi    = 1.75e-03
594               zclei    = zcshi
595               
596               !  sensible and latent fluxes over ice
597               zrhova     = zrhoa(ji,jj) * sf(jp_wndm)%fnow(ji,jj)      ! computation of intermediate values
598               zrhovaclei = zrhova * zcshi * 2.834e+06
599               zrhovacshi = zrhova * zclei * 1004.0
600           
601               !  sensible heat flux
602               z_qsb(ji,jj,jl) = zrhovacshi * ( pst(ji,jj,jl) - ztatm(ji,jj) )
603           
604               !  latent heat flux
605               p_qla(ji,jj,jl) = MAX(  0.e0, zrhovaclei * ( zqsati - zqatm(ji,jj) )  )
606             
607               !  sensitivity of non solar fluxes (dQ/dT) (long-wave, sensible and latent fluxes)
608               zdqlw = 4.0 * emic * stefan * ztice3
609               zdqsb = zrhovacshi
610               zdqla = zrhovaclei * ( zdesidt * ( zqsati * zqsati / ( zesi * zesi ) ) * ( zpatm / 0.622 ) )   
611               !
612               p_dqla(ji,jj,jl) = zdqla                           ! latent flux sensitivity
613               p_dqns(ji,jj,jl) = -( zdqlw + zdqsb + zdqla )      !  total non solar sensitivity
614            END DO
615            !
616         END DO
617         !
618      END DO
619      !
620      ! ----------------------------------------------------------------------------- !
621      !    Total FLUXES                                                       !
622      ! ----------------------------------------------------------------------------- !
623      !
624!CDIR COLLAPSE
625      p_qns(:,:,:) = z_qlw (:,:,:) - z_qsb (:,:,:) - p_qla (:,:,:)      ! Downward Non Solar flux
626!CDIR COLLAPSE
627      p_tpr(:,:)   = sf(jp_prec)%fnow(:,:) / rday                       ! total precipitation [kg/m2/s]
628      !
629!!gm : not necessary as all input data are lbc_lnk...
630      CALL lbc_lnk( p_fr1  (:,:) , 'T', 1. )
631      CALL lbc_lnk( p_fr2  (:,:) , 'T', 1. )
632      DO jl = 1, ijpl
633         CALL lbc_lnk( p_qns (:,:,jl) , 'T', 1. )
634         CALL lbc_lnk( p_dqns(:,:,jl) , 'T', 1. )
635         CALL lbc_lnk( p_qla (:,:,jl) , 'T', 1. )
636         CALL lbc_lnk( p_dqla(:,:,jl) , 'T', 1. )
637      END DO
638
639!!gm : mask is not required on forcing
640      DO jl = 1, ijpl
641         p_qns (:,:,jl) = p_qns (:,:,jl) * tmask(:,:,1)
642         p_qla (:,:,jl) = p_qla (:,:,jl) * tmask(:,:,1)
643         p_dqns(:,:,jl) = p_dqns(:,:,jl) * tmask(:,:,1)
644         p_dqla(:,:,jl) = p_dqla(:,:,jl) * tmask(:,:,1)
645      END DO
646
647      IF(ln_ctl) THEN
648         CALL prt_ctl(tab3d_1=z_qsb  , clinfo1=' blk_ice_clio: z_qsb  : ', tab3d_2=z_qlw  , clinfo2=' z_qlw  : ', kdim=ijpl)
649         CALL prt_ctl(tab3d_1=p_qla  , clinfo1=' blk_ice_clio: z_qla  : ', tab3d_2=p_qsr  , clinfo2=' p_qsr  : ', kdim=ijpl)
650         CALL prt_ctl(tab3d_1=p_dqns , clinfo1=' blk_ice_clio: p_dqns : ', tab3d_2=p_qns  , clinfo2=' p_qns  : ', kdim=ijpl)
651         CALL prt_ctl(tab3d_1=p_dqla , clinfo1=' blk_ice_clio: p_dqla : ', tab3d_2=pst    , clinfo2=' pst    : ', kdim=ijpl)
652         CALL prt_ctl(tab2d_1=p_tpr  , clinfo1=' blk_ice_clio: p_tpr  : ', tab2d_2=p_spr  , clinfo2=' p_spr  : ')
653         CALL prt_ctl(tab2d_1=p_taui , clinfo1=' blk_ice_clio: p_taui : ', tab2d_2=p_tauj , clinfo2=' p_tauj : ')
654      ENDIF
655
656
657   END SUBROUTINE blk_ice_clio
658
659
660   SUBROUTINE blk_clio_qsr_oce( pqsr_oce )
661      !!---------------------------------------------------------------------------
662      !!                     ***  ROUTINE blk_clio_qsr_oce  ***
663      !!                 
664      !!  ** Purpose :   Computation of the shortwave radiation at the ocean and the
665      !!               snow/ice surfaces.
666      !!         
667      !!  ** Method  : - computed qsr from the cloud cover for both ice and ocean
668      !!               - also initialise sbudyko and stauc once for all
669      !!----------------------------------------------------------------------
670      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)     ::   pqsr_oce    ! shortwave radiation  over the ocean
671      !!
672      INTEGER, PARAMETER  ::   jp24 = 24   ! sampling of the daylight period (sunrise to sunset) into 24 equal parts
673      !!     
674      INTEGER  ::   ji, jj, jt    ! dummy loop indices
675      INTEGER  ::   indaet            !  = -1, 0, 1 for odd, normal and leap years resp.
676      INTEGER  ::   iday              ! integer part of day
677      INTEGER  ::   indxb, indxc      ! index for cloud depth coefficient
678
679      REAL(wp)  ::   zalat , zclat, zcmue, zcmue2    ! local scalars
680      REAL(wp)  ::   zmt1, zmt2, zmt3                !
681      REAL(wp)  ::   zdecl, zsdecl , zcdecl          !
682      REAL(wp)  ::   za_oce, ztamr                   !
683
684      REAL(wp) ::   zdl, zlha                        ! local scalars
685      REAL(wp) ::   zlmunoon, zcldcor, zdaycor       !   
686      REAL(wp) ::   zxday, zdist, zcoef, zcoef1      !
687      REAL(wp) ::   zes
688      !!
689      REAL(wp), DIMENSION(jpi,jpj) ::   zev          ! vapour pressure
690      REAL(wp), DIMENSION(jpi,jpj) ::   zdlha, zlsrise, zlsset     ! 2D workspace
691
692      REAL(wp), DIMENSION(jpi,jpj) ::   zps, zpc   ! sine (cosine) of latitude per sine (cosine) of solar declination
693      !!---------------------------------------------------------------------
694
695
696      IF( lbulk_init ) THEN             !   Initilization at first time step only
697         rdtbs2 = nn_fsbc * rdt * 0.5
698         ! cloud optical depths as a function of latitude (Chou et al., 1981).
699         ! and the correction factor for taking into account  the effect of clouds
700         DO jj = 1, jpj
701            DO ji = 1 , jpi
702               zalat          = ( 90.e0 - ABS( gphit(ji,jj) ) ) /  5.e0
703               zclat          = ( 95.e0 -      gphit(ji,jj)   ) / 10.e0
704               indxb          = 1 + INT( zalat )
705               indxc          = 1 + INT( zclat )
706               zdl            = zclat - INT( zclat )
707               !  correction factor to account for the effect of clouds
708               sbudyko(ji,jj) = budyko(indxb)
709               stauc  (ji,jj) = ( 1.e0 - zdl ) * tauco( indxc ) + zdl * tauco( indxc + 1 )
710            END DO
711         END DO
712         IF    ( nleapy == 1 ) THEN   ;   yearday = 366.e0
713         ELSEIF( nleapy == 0 ) THEN   ;   yearday = 365.e0
714         ELSEIF( nleapy == 30) THEN   ;   yearday = 360.e0
715         ENDIF
716         lbulk_init = .FALSE.
717      ENDIF
718
719
720      ! Saturated water vapour and vapour pressure
721      ! ------------------------------------------
722!CDIR NOVERRCHK
723!CDIR COLLAPSE
724      DO jj = 1, jpj
725!CDIR NOVERRCHK
726         DO ji = 1, jpi
727            ztamr = sf(jp_tair)%fnow(ji,jj) + rt0 - rtt
728            zmt1  = SIGN( 17.269,  ztamr )
729            zmt2  = SIGN( 21.875,  ztamr )
730            zmt3  = SIGN( 28.200, -ztamr )
731            zes = 611.0 * EXP(  ABS( ztamr ) * MIN ( zmt1, zmt2 )   &              ! Saturation water vapour
732               &                     / ( sf(jp_tair)%fnow(ji,jj) + rt0 - 35.86  + MAX( 0.e0, zmt3 ) )  )
733            zev(ji,jj) = sf(jp_humi)%fnow(ji,jj) * zes * 1.0e-05                   ! vapour pressure 
734         END DO
735      END DO
736
737      !-----------------------------------!
738      !  Computation of solar irradiance  !
739      !-----------------------------------!
740!!gm : hard coded  leap year ???
741      indaet   = 1                                    ! = -1, 0, 1 for odd, normal and leap years resp.
742      zxday = nday_year + rdtbs2 / rday               ! day of the year at which the fluxes are calculated
743      iday  = INT( zxday )                            ! (centred at the middle of the ice time step)
744      CALL flx_blk_declin( indaet, iday, zdecl )      ! solar declination of the current day
745      zsdecl = SIN( zdecl * rad )                     ! its sine
746      zcdecl = COS( zdecl * rad )                     ! its cosine
747
748
749      !  correction factor added for computation of shortwave flux to take into account the variation of
750      !  the distance between the sun and the earth during the year (Oberhuber 1988)
751      zdist    = zxday * 2. * rpi / yearday
752      zdaycor  = 1.0 + 0.0013 * SIN( zdist ) + 0.0342 * COS( zdist )
753
754!CDIR NOVERRCHK
755      DO jj = 1, jpj
756!CDIR NOVERRCHK
757         DO ji = 1, jpi
758            !  product of sine (cosine) of latitude and sine (cosine) of solar declination
759            zps(ji,jj) = SIN( gphit(ji,jj) * rad ) * zsdecl
760            zpc(ji,jj) = COS( gphit(ji,jj) * rad ) * zcdecl
761            !  computation of the both local time of sunrise and sunset
762            zlsrise(ji,jj) = ACOS( - SIGN( 1.e0, zps(ji,jj) )    &
763               &                   * MIN(  1.e0, SIGN( 1.e0, zps(ji,jj) ) * ( zps(ji,jj) / zpc(ji,jj) )  )   )
764            zlsset (ji,jj) = - zlsrise(ji,jj)
765            !  dividing the solar day into jp24 segments of length zdlha
766            zdlha  (ji,jj) = ( zlsrise(ji,jj) - zlsset(ji,jj) ) / REAL( jp24 )
767         END DO
768      END DO
769
770
771      !---------------------------------------------!
772      !  shortwave radiation absorbed by the ocean  !
773      !---------------------------------------------!
774      pqsr_oce(:,:)   = 0.e0      ! set ocean qsr to zero     
775
776      ! compute and sum ocean qsr over the daylight (i.e. between sunrise and sunset)
777!CDIR NOVERRCHK   
778      DO jt = 1, jp24
779         zcoef = FLOAT( jt ) - 0.5
780!CDIR NOVERRCHK     
781!CDIR COLLAPSE
782         DO jj = 1, jpj
783!CDIR NOVERRCHK
784            DO ji = 1, jpi
785               zlha = COS(  zlsrise(ji,jj) - zcoef * zdlha(ji,jj)  )                  ! local hour angle
786               zcmue              = MAX( 0.e0 ,   zps(ji,jj) + zpc(ji,jj) * zlha  )   ! cos of local solar altitude
787               zcmue2             = 1368.0 * zcmue * zcmue
788
789               ! ocean albedo depending on the cloud cover (Payne, 1972)
790               za_oce     = ( 1.0 - sf(jp_ccov)%fnow(ji,jj) ) * 0.05 / ( 1.1 * zcmue**1.4 + 0.15 )   &   ! clear sky
791                  &       +         sf(jp_ccov)%fnow(ji,jj)   * 0.06                                     ! overcast
792
793                  ! solar heat flux absorbed by the ocean (Zillman, 1972)
794               pqsr_oce(ji,jj) = pqsr_oce(ji,jj)                                         &
795                  &            + ( 1.0 - za_oce ) * zdlha(ji,jj) * zcmue2                &
796                  &            / ( ( zcmue + 2.7 ) * zev(ji,jj) + 1.085 * zcmue +  0.10 )
797            END DO
798         END DO
799      END DO
800      ! Taking into account the ellipsity of the earth orbit, the clouds AND masked if sea-ice cover > 0%
801      zcoef1 = srgamma * zdaycor / ( 2. * rpi )
802!CDIR COLLAPSE
803      DO jj = 1, jpj
804         DO ji = 1, jpi
805            zlmunoon = ASIN( zps(ji,jj) + zpc(ji,jj) ) / rad                         ! local noon solar altitude
806            zcldcor  = MIN(  1.e0, ( 1.e0 - 0.62 * sf(jp_ccov)%fnow(ji,jj)   &       ! cloud correction (Reed 1977)
807               &                          + 0.0019 * zlmunoon )                 )
808            pqsr_oce(ji,jj) = zcoef1 * zcldcor * pqsr_oce(ji,jj) * tmask(ji,jj,1)   ! and zcoef1: ellipsity
809         END DO
810      END DO
811
812   END SUBROUTINE blk_clio_qsr_oce
813
814
815   SUBROUTINE blk_clio_qsr_ice( pa_ice_cs, pa_ice_os, pqsr_ice )
816      !!---------------------------------------------------------------------------
817      !!                     ***  ROUTINE blk_clio_qsr_ice  ***
818      !!                 
819      !!  ** Purpose :   Computation of the shortwave radiation at the ocean and the
820      !!               snow/ice surfaces.
821      !!         
822      !!  ** Method  : - computed qsr from the cloud cover for both ice and ocean
823      !!               - also initialise sbudyko and stauc once for all
824      !!----------------------------------------------------------------------
825      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pa_ice_cs   ! albedo of ice under clear sky
826      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pa_ice_os   ! albedo of ice under overcast sky
827      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pqsr_ice    ! shortwave radiation over the ice/snow
828      !!
829      INTEGER, PARAMETER  ::   jp24 = 24   ! sampling of the daylight period (sunrise to sunset) into 24 equal parts
830      !!
831      INTEGER  ::   ji, jj, jl, jt    ! dummy loop indices
832      INTEGER  ::   ijpl              ! number of ice categories (3rd dim of pqsr_ice)
833      INTEGER  ::   indaet            !  = -1, 0, 1 for odd, normal and leap years resp.
834      INTEGER  ::   iday              ! integer part of day
835      !!
836      REAL(wp) ::   zcmue, zcmue2, ztamr          ! temporary scalars
837      REAL(wp) ::   zmt1, zmt2, zmt3              !    -         -
838      REAL(wp) ::   zdecl, zsdecl, zcdecl         !    -         -
839      REAL(wp) ::   zlha, zdaycor, zes            !    -         -
840      REAL(wp) ::   zxday, zdist, zcoef, zcoef1   !    -         -
841      REAL(wp) ::   zqsr_ice_cs, zqsr_ice_os      !    -         -
842      !!
843      REAL(wp), DIMENSION(jpi,jpj) ::   zev                      ! vapour pressure
844      REAL(wp), DIMENSION(jpi,jpj) ::   zdlha, zlsrise, zlsset   ! 2D workspace
845      REAL(wp), DIMENSION(jpi,jpj) ::   zps, zpc   ! sine (cosine) of latitude per sine (cosine) of solar declination
846      !!---------------------------------------------------------------------
847
848      ijpl = SIZE(pqsr_ice, 3 )      ! number of ice categories
849     
850      ! Saturated water vapour and vapour pressure
851      ! ------------------------------------------
852!CDIR NOVERRCHK
853!CDIR COLLAPSE
854      DO jj = 1, jpj
855!CDIR NOVERRCHK
856         DO ji = 1, jpi           
857            ztamr = sf(jp_tair)%fnow(ji,jj) + rt0 - rtt           
858            zmt1  = SIGN( 17.269,  ztamr )
859            zmt2  = SIGN( 21.875,  ztamr )
860            zmt3  = SIGN( 28.200, -ztamr )
861            zes = 611.0 * EXP(  ABS( ztamr ) * MIN ( zmt1, zmt2 )   &              ! Saturation water vapour
862               &                     / ( sf(jp_tair)%fnow(ji,jj) + rt0 - 35.86  + MAX( 0.e0, zmt3 ) )  )
863            zev(ji,jj) = sf(jp_humi)%fnow(ji,jj) * zes * 1.0e-05                   ! vapour pressure 
864         END DO
865      END DO
866
867      !-----------------------------------!
868      !  Computation of solar irradiance  !
869      !-----------------------------------!
870!!gm : hard coded  leap year ???
871      indaet   = 1                                    ! = -1, 0, 1 for odd, normal and leap years resp.
872      zxday = nday_year + rdtbs2 / rday               ! day of the year at which the fluxes are calculated
873      iday  = INT( zxday )                            ! (centred at the middle of the ice time step)
874      CALL flx_blk_declin( indaet, iday, zdecl )      ! solar declination of the current day
875      zsdecl = SIN( zdecl * rad )                     ! its sine
876      zcdecl = COS( zdecl * rad )                     ! its cosine
877
878     
879      !  correction factor added for computation of shortwave flux to take into account the variation of
880      !  the distance between the sun and the earth during the year (Oberhuber 1988)
881      zdist    = zxday * 2. * rpi / yearday
882      zdaycor  = 1.0 + 0.0013 * SIN( zdist ) + 0.0342 * COS( zdist )
883
884!CDIR NOVERRCHK
885      DO jj = 1, jpj
886!CDIR NOVERRCHK
887         DO ji = 1, jpi
888            !  product of sine (cosine) of latitude and sine (cosine) of solar declination
889            zps(ji,jj) = SIN( gphit(ji,jj) * rad ) * zsdecl
890            zpc(ji,jj) = COS( gphit(ji,jj) * rad ) * zcdecl
891            !  computation of the both local time of sunrise and sunset
892            zlsrise(ji,jj) = ACOS( - SIGN( 1.e0, zps(ji,jj) )    &
893               &                   * MIN(  1.e0, SIGN( 1.e0, zps(ji,jj) ) * ( zps(ji,jj) / zpc(ji,jj) )  )   ) 
894            zlsset (ji,jj) = - zlsrise(ji,jj)
895            !  dividing the solar day into jp24 segments of length zdlha
896            zdlha  (ji,jj) = ( zlsrise(ji,jj) - zlsset(ji,jj) ) / REAL( jp24 )
897         END DO
898      END DO
899
900
901      !---------------------------------------------!
902      !  shortwave radiation absorbed by the ice    !
903      !---------------------------------------------!
904      ! compute and sum ice qsr over the daylight for each ice categories
905      pqsr_ice(:,:,:) = 0.e0
906      zcoef1 = zdaycor / ( 2. * rpi )       ! Correction for the ellipsity of the earth orbit
907     
908      !                    !----------------------------!
909      DO jl = 1, ijpl      !  loop over ice categories  !
910         !                 !----------------------------!
911!CDIR NOVERRCHK   
912         DO jt = 1, jp24   
913            zcoef = FLOAT( jt ) - 0.5
914!CDIR NOVERRCHK     
915!CDIR COLLAPSE
916            DO jj = 1, jpj
917!CDIR NOVERRCHK
918               DO ji = 1, jpi
919                  zlha = COS(  zlsrise(ji,jj) - zcoef * zdlha(ji,jj)  )                  ! local hour angle
920                  zcmue              = MAX( 0.e0 ,   zps(ji,jj) + zpc(ji,jj) * zlha  )   ! cos of local solar altitude
921                  zcmue2             = 1368.0 * zcmue * zcmue
922                 
923                  !  solar heat flux absorbed by the ice/snow system (Shine and Crane 1984 adapted to high albedo)
924                  zqsr_ice_cs =  ( 1.0 - pa_ice_cs(ji,jj,jl) ) * zdlha(ji,jj) * zcmue2        &   ! clear sky
925                     &        / ( ( 1.0 + zcmue ) * zev(ji,jj) + 1.2 * zcmue + 0.0455 )
926                  zqsr_ice_os = zdlha(ji,jj) * SQRT( zcmue )                                  &   ! overcast sky
927                     &        * ( 53.5 + 1274.5 * zcmue )      * ( 1.0 - 0.996  * pa_ice_os(ji,jj,jl) )    &
928                     &        / (  1.0 + 0.139  * stauc(ji,jj) * ( 1.0 - 0.9435 * pa_ice_os(ji,jj,jl) ) )       
929             
930                  pqsr_ice(ji,jj,jl) = pqsr_ice(ji,jj,jl) + (  ( 1.0 - sf(jp_ccov)%fnow(ji,jj) ) * zqsr_ice_cs    &
931                     &                                       +         sf(jp_ccov)%fnow(ji,jj)   * zqsr_ice_os  )
932               END DO
933            END DO
934         END DO
935         !
936         ! Correction : Taking into account the ellipsity of the earth orbit
937         pqsr_ice(:,:,jl) = pqsr_ice(:,:,jl) * zcoef1 * tmask(:,:,1)
938         !
939         !                 !--------------------------------!
940      END DO               !  end loop over ice categories  !
941      !                    !--------------------------------!
942
943
944!!gm  : this should be suppress as input data have been passed through lbc_lnk
945      DO jl = 1, ijpl
946         CALL lbc_lnk( pqsr_ice(:,:,jl) , 'T', 1. )
947      END DO
948      !
949   END SUBROUTINE blk_clio_qsr_ice
950
951
952   SUBROUTINE flx_blk_declin( ky, kday, pdecl )
953      !!---------------------------------------------------------------------------
954      !!               ***  ROUTINE flx_blk_declin  ***
955      !!         
956      !! ** Purpose :   Computation of the solar declination for the day
957      !!       
958      !! ** Method  :   ???
959      !!---------------------------------------------------------------------
960      INTEGER , INTENT(in   ) ::   ky      ! = -1, 0, 1 for odd, normal and leap years resp.
961      INTEGER , INTENT(in   ) ::   kday    ! day of the year ( kday = 1 on january 1)
962      REAL(wp), INTENT(  out) ::   pdecl   ! solar declination
963      !!
964      REAL(wp) ::   a0  =  0.39507671      ! coefficients for solar declinaison computation
965      REAL(wp) ::   a1  = 22.85684301      !     "              ""                 "
966      REAL(wp) ::   a2  = -0.38637317      !     "              ""                 "
967      REAL(wp) ::   a3  =  0.15096535      !     "              ""                 "
968      REAL(wp) ::   a4  = -0.00961411      !     "              ""                 "
969      REAL(wp) ::   b1  = -4.29692073      !     "              ""                 "
970      REAL(wp) ::   b2  =  0.05702074      !     "              ""                 "
971      REAL(wp) ::   b3  = -0.09028607      !     "              ""                 "
972      REAL(wp) ::   b4  =  0.00592797
973      !!
974      REAL(wp) ::   zday   ! corresponding day of type year (cf. ky)
975      REAL(wp) ::   zp     ! temporary scalars
976      !!---------------------------------------------------------------------
977           
978      IF    ( ky == 1 )  THEN   ;   zday = REAL( kday ) - 0.5
979      ELSEIF( ky == 3 )  THEN   ;   zday = REAL( kday ) - 1.
980      ELSE                      ;   zday = REAL( kday )
981      ENDIF
982     
983      zp = rpi * ( 2.0 * zday - 367.0 ) / yearday
984     
985      pdecl  = a0                                                                      &
986         &   + a1 * COS( zp ) + a2 * COS( 2. * zp ) + a3 * COS( 3. * zp ) + a4 * COS( 4. * zp )   &
987         &   + b1 * SIN( zp ) + b2 * SIN( 2. * zp ) + b3 * SIN( 3. * zp ) + b4 * SIN( 4. * zp )
988      !
989   END SUBROUTINE flx_blk_declin
990
991   !!======================================================================
992END MODULE sbcblk_clio
Note: See TracBrowser for help on using the repository browser.