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

Last change on this file since 3186 was 3186, checked in by smasson, 10 years ago

dev_NEMO_MERGE_2011: replace the old wrk_nemo with the new wrk_nemo

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