source: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90 @ 5126

Last change on this file since 5126 was 5126, checked in by vancop, 6 years ago

LIM3 bugfix: did not match with sbcblk_clio.F90

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