source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90 @ 3229

Last change on this file since 3229 was 3229, checked in by charris, 10 years ago

Added timing calls to most significant routines in LDF, SBC and ZDF.

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