source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90 @ 5443

Last change on this file since 5443 was 5443, checked in by davestorkey, 5 years ago

Update 2015/dev_r5021_UKMO1_CICE_coupling branch to revision 5442 of the trunk.

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