source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90 @ 8529

Last change on this file since 8529 was 8529, checked in by cbricaud, 3 years ago

fix ticket #1667 in nemo_v3_6_STABLE

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