New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbcblk_clio.F90 in branches/UKMO/dev_r5107_hadgem3_cplfld/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/dev_r5107_hadgem3_cplfld/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90 @ 5473

Last change on this file since 5473 was 5473, checked in by cguiavarch, 9 years ago

Clear svn keywords from UKMO/dev_r5107_hadgem3_cplfld

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