source: branches/2015/dev_r5218_CNRS17_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90 @ 5362

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

update CLIO forcing to match with new standards

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