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

Last change on this file since 5357 was 5357, checked in by clem, 6 years ago

LIM3: change the interface between the ice and atm for both coupled and forced modes. Some work still needs to be done to deal with sublimation in coupled mode.

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