source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90 @ 10395

Last change on this file since 10395 was 8058, checked in by jgraham, 4 years ago

Clear keywords

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