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

source: branches/devmercator2010_1/NEMO/OPA_SRC/SBC/sbcblk_clio.F90 @ 2134

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

remove debug print and comment

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