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 branches/devmercator2010/NEMO/OPA_SRC/SBC – NEMO

source: branches/devmercator2010/NEMO/OPA_SRC/SBC/sbcblk_clio.F90 @ 2077

Last change on this file since 2077 was 2077, checked in by cbricaud, 13 years ago

commit change from DEV_r1784_3DF

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