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_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90 @ 10775

Last change on this file since 10775 was 10759, checked in by andmirek, 5 years ago

GMED 450 write output.namelist.dyn only for nprint > 2

File size: 57.5 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 .AND. nprint > 2) 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      zsnw(:,:) = 0._wp
666      CALL lim_thd_snwblow( pfrld, zsnw )          ! snow redistribution by wind
667      emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * ( 1._wp - zsnw )
668      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
669      emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
670
671      ! --- heat flux associated with emp --- !
672      qemp_oce(:,:) = - pfrld(:,:) * zevap(:,:) * sst_m(:,:) * rcp                               & ! evap
673         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip
674         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip
675         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
676      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only)
677         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
678
679      ! --- total solar and non solar fluxes --- !
680      qns_tot(:,:) = pfrld(:,:) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) + qemp_ice(:,:) + qemp_oce(:,:)
681      qsr_tot(:,:) = pfrld(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
682
683      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
684      qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
685
686      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- !
687      DO jl = 1, jpl
688         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) - lfus )
689                                   ! but then qemp_ice should also include sublimation
690      END DO
691
692      CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 
693#endif
694
695!!gm : not necessary as all input data are lbc_lnk...
696      CALL lbc_lnk( fr1_i0  (:,:) , 'T', 1. )
697      CALL lbc_lnk( fr2_i0  (:,:) , 'T', 1. )
698      DO jl = 1, jpl
699         CALL lbc_lnk( qns_ice (:,:,jl) , 'T', 1. )
700         CALL lbc_lnk( dqns_ice(:,:,jl) , 'T', 1. )
701         CALL lbc_lnk( qla_ice (:,:,jl) , 'T', 1. )
702         CALL lbc_lnk( dqla_ice(:,:,jl) , 'T', 1. )
703      END DO
704
705!!gm : mask is not required on forcing
706      DO jl = 1, jpl
707         qns_ice (:,:,jl) = qns_ice (:,:,jl) * tmask(:,:,1)
708         qla_ice (:,:,jl) = qla_ice (:,:,jl) * tmask(:,:,1)
709         dqns_ice(:,:,jl) = dqns_ice(:,:,jl) * tmask(:,:,1)
710         dqla_ice(:,:,jl) = dqla_ice(:,:,jl) * tmask(:,:,1)
711      END DO
712
713      CALL wrk_dealloc( jpi,jpj, ztatm, zqatm, zevsqr, zrhoa )
714      CALL wrk_dealloc( jpi,jpj, jpl  , z_qlw, z_qsb )
715
716      IF(ln_ctl) THEN
717         CALL prt_ctl(tab3d_1=z_qsb  , clinfo1=' blk_ice_clio: z_qsb  : ', tab3d_2=z_qlw  , clinfo2=' z_qlw  : ', kdim=jpl)
718         CALL prt_ctl(tab3d_1=qla_ice  , clinfo1=' blk_ice_clio: z_qla  : ', tab3d_2=qsr_ice  , clinfo2=' qsr_ice  : ', kdim=jpl)
719         CALL prt_ctl(tab3d_1=dqns_ice , clinfo1=' blk_ice_clio: dqns_ice : ', tab3d_2=qns_ice  , clinfo2=' qns_ice  : ', kdim=jpl)
720         CALL prt_ctl(tab3d_1=dqla_ice , clinfo1=' blk_ice_clio: dqla_ice : ', tab3d_2=ptsu    , clinfo2=' ptsu    : ', kdim=jpl)
721         CALL prt_ctl(tab2d_1=tprecip  , clinfo1=' blk_ice_clio: tprecip  : ', tab2d_2=sprecip  , clinfo2=' sprecip  : ')
722      ENDIF
723
724      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_clio_flx')
725      !
726   END SUBROUTINE blk_ice_clio_flx
727
728#endif
729
730   SUBROUTINE blk_clio_qsr_oce( pqsr_oce )
731      !!---------------------------------------------------------------------------
732      !!                     ***  ROUTINE blk_clio_qsr_oce  ***
733      !!                 
734      !!  ** Purpose :   Computation of the shortwave radiation at the ocean and the
735      !!               snow/ice surfaces.
736      !!         
737      !!  ** Method  : - computed qsr from the cloud cover for both ice and ocean
738      !!               - also initialise sbudyko and stauc once for all
739      !!----------------------------------------------------------------------
740      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)     ::   pqsr_oce    ! shortwave radiation  over the ocean
741      !!
742      INTEGER, PARAMETER  ::   jp24 = 24   ! sampling of the daylight period (sunrise to sunset) into 24 equal parts
743      !!     
744      INTEGER  ::   ji, jj, jt    ! dummy loop indices
745      INTEGER  ::   indaet            !  = -1, 0, 1 for odd, normal and leap years resp.
746      INTEGER  ::   iday              ! integer part of day
747      INTEGER  ::   indxb, indxc      ! index for cloud depth coefficient
748
749      REAL(wp)  ::   zalat , zclat, zcmue, zcmue2    ! local scalars
750      REAL(wp)  ::   zmt1, zmt2, zmt3                !
751      REAL(wp)  ::   zdecl, zsdecl , zcdecl          !
752      REAL(wp)  ::   za_oce, ztamr                   !
753
754      REAL(wp) ::   zdl, zlha                        ! local scalars
755      REAL(wp) ::   zlmunoon, zcldcor, zdaycor       !   
756      REAL(wp) ::   zxday, zdist, zcoef, zcoef1      !
757      REAL(wp) ::   zes
758     
759      REAL(wp), DIMENSION(:,:), POINTER ::   zev          ! vapour pressure
760      REAL(wp), DIMENSION(:,:), POINTER ::   zdlha, zlsrise, zlsset     ! 2D workspace
761      REAL(wp), DIMENSION(:,:), POINTER ::   zps, zpc   ! sine (cosine) of latitude per sine (cosine) of solar declination
762      !!---------------------------------------------------------------------
763      !
764      IF( nn_timing == 1 )  CALL timing_start('blk_clio_qsr_oce')
765      !
766      CALL wrk_alloc( jpi,jpj, zev, zdlha, zlsrise, zlsset, zps, zpc )
767
768      IF( lbulk_init ) THEN             !   Initilization at first time step only
769         rdtbs2 = nn_fsbc * rdt * 0.5
770         ! cloud optical depths as a function of latitude (Chou et al., 1981).
771         ! and the correction factor for taking into account  the effect of clouds
772         DO jj = 1, jpj
773            DO ji = 1 , jpi
774               zalat          = ( 90.e0 - ABS( gphit(ji,jj) ) ) /  5.e0
775               zclat          = ( 95.e0 -      gphit(ji,jj)   ) / 10.e0
776               indxb          = 1 + INT( zalat )
777               indxc          = 1 + INT( zclat )
778               zdl            = zclat - INT( zclat )
779               !  correction factor to account for the effect of clouds
780               sbudyko(ji,jj) = budyko(indxb)
781               stauc  (ji,jj) = ( 1.e0 - zdl ) * tauco( indxc ) + zdl * tauco( indxc + 1 )
782            END DO
783         END DO
784         lbulk_init = .FALSE.
785      ENDIF
786
787
788      ! Saturated water vapour and vapour pressure
789      ! ------------------------------------------
790!CDIR NOVERRCHK
791!CDIR COLLAPSE
792      DO jj = 1, jpj
793!CDIR NOVERRCHK
794         DO ji = 1, jpi
795            ztamr = sf(jp_tair)%fnow(ji,jj,1) - rtt
796            zmt1  = SIGN( 17.269,  ztamr )
797            zmt2  = SIGN( 21.875,  ztamr )
798            zmt3  = SIGN( 28.200, -ztamr )
799            zes = 611.0 * EXP(  ABS( ztamr ) * MIN ( zmt1, zmt2 )   &              ! Saturation water vapour
800               &                     / ( sf(jp_tair)%fnow(ji,jj,1) - 35.86  + MAX( 0.e0, zmt3 ) )  )
801            zev(ji,jj) = sf(jp_humi)%fnow(ji,jj,1) * zes * 1.0e-05                 ! vapour pressure 
802         END DO
803      END DO
804
805      !-----------------------------------!
806      !  Computation of solar irradiance  !
807      !-----------------------------------!
808!!gm : hard coded  leap year ???
809      indaet   = 1                                    ! = -1, 0, 1 for odd, normal and leap years resp.
810      zxday = nday_year + rdtbs2 / rday               ! day of the year at which the fluxes are calculated
811      iday  = INT( zxday )                            ! (centred at the middle of the ice time step)
812      CALL flx_blk_declin( indaet, iday, zdecl )      ! solar declination of the current day
813      zsdecl = SIN( zdecl * rad )                     ! its sine
814      zcdecl = COS( zdecl * rad )                     ! its cosine
815
816
817      !  correction factor added for computation of shortwave flux to take into account the variation of
818      !  the distance between the sun and the earth during the year (Oberhuber 1988)
819      zdist    = zxday * 2. * rpi / REAL(nyear_len(1), wp)
820      zdaycor  = 1.0 + 0.0013 * SIN( zdist ) + 0.0342 * COS( zdist )
821
822!CDIR NOVERRCHK
823      DO jj = 1, jpj
824!CDIR NOVERRCHK
825         DO ji = 1, jpi
826            !  product of sine (cosine) of latitude and sine (cosine) of solar declination
827            zps(ji,jj) = SIN( gphit(ji,jj) * rad ) * zsdecl
828            zpc(ji,jj) = COS( gphit(ji,jj) * rad ) * zcdecl
829            !  computation of the both local time of sunrise and sunset
830            zlsrise(ji,jj) = ACOS( - SIGN( 1.e0, zps(ji,jj) )    &
831               &                   * MIN(  1.e0, SIGN( 1.e0, zps(ji,jj) ) * ( zps(ji,jj) / zpc(ji,jj) )  )   )
832            zlsset (ji,jj) = - zlsrise(ji,jj)
833            !  dividing the solar day into jp24 segments of length zdlha
834            zdlha  (ji,jj) = ( zlsrise(ji,jj) - zlsset(ji,jj) ) / REAL( jp24, wp )
835         END DO
836      END DO
837
838
839      !---------------------------------------------!
840      !  shortwave radiation absorbed by the ocean  !
841      !---------------------------------------------!
842      pqsr_oce(:,:)   = 0.e0      ! set ocean qsr to zero     
843
844      ! compute and sum ocean qsr over the daylight (i.e. between sunrise and sunset)
845!CDIR NOVERRCHK   
846      DO jt = 1, jp24
847         zcoef = FLOAT( jt ) - 0.5
848!CDIR NOVERRCHK     
849!CDIR COLLAPSE
850         DO jj = 1, jpj
851!CDIR NOVERRCHK
852            DO ji = 1, jpi
853               zlha = COS(  zlsrise(ji,jj) - zcoef * zdlha(ji,jj)  )                  ! local hour angle
854               zcmue              = MAX( 0.e0 ,   zps(ji,jj) + zpc(ji,jj) * zlha  )   ! cos of local solar altitude
855               zcmue2             = 1368.0 * zcmue * zcmue
856
857               ! ocean albedo depending on the cloud cover (Payne, 1972)
858               za_oce     = ( 1.0 - sf(jp_ccov)%fnow(ji,jj,1) ) * 0.05 / ( 1.1 * zcmue**1.4 + 0.15 )   &   ! clear sky
859                  &       +         sf(jp_ccov)%fnow(ji,jj,1)   * 0.06                                     ! overcast
860
861                  ! solar heat flux absorbed by the ocean (Zillman, 1972)
862               pqsr_oce(ji,jj) = pqsr_oce(ji,jj)                                         &
863                  &            + ( 1.0 - za_oce ) * zdlha(ji,jj) * zcmue2                &
864                  &            / ( ( zcmue + 2.7 ) * zev(ji,jj) + 1.085 * zcmue +  0.10 )
865            END DO
866         END DO
867      END DO
868      ! Taking into account the ellipsity of the earth orbit, the clouds AND masked if sea-ice cover > 0%
869      zcoef1 = srgamma * zdaycor / ( 2. * rpi )
870!CDIR COLLAPSE
871      DO jj = 1, jpj
872         DO ji = 1, jpi
873            zlmunoon = ASIN( zps(ji,jj) + zpc(ji,jj) ) / rad                         ! local noon solar altitude
874            zcldcor  = MIN(  1.e0, ( 1.e0 - 0.62 * sf(jp_ccov)%fnow(ji,jj,1)   &     ! cloud correction (Reed 1977)
875               &                          + 0.0019 * zlmunoon )                 )
876            pqsr_oce(ji,jj) = zcoef1 * zcldcor * pqsr_oce(ji,jj) * tmask(ji,jj,1)    ! and zcoef1: ellipsity
877         END DO
878      END DO
879
880      CALL wrk_dealloc( jpi,jpj, zev, zdlha, zlsrise, zlsset, zps, zpc )
881      !
882      IF( nn_timing == 1 )  CALL timing_stop('blk_clio_qsr_oce')
883      !
884   END SUBROUTINE blk_clio_qsr_oce
885
886
887   SUBROUTINE blk_clio_qsr_ice( pa_ice_cs, pa_ice_os, pqsr_ice )
888      !!---------------------------------------------------------------------------
889      !!                     ***  ROUTINE blk_clio_qsr_ice  ***
890      !!                 
891      !!  ** Purpose :   Computation of the shortwave radiation at the ocean and the
892      !!               snow/ice surfaces.
893      !!         
894      !!  ** Method  : - computed qsr from the cloud cover for both ice and ocean
895      !!               - also initialise sbudyko and stauc once for all
896      !!----------------------------------------------------------------------
897      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pa_ice_cs   ! albedo of ice under clear sky
898      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pa_ice_os   ! albedo of ice under overcast sky
899      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pqsr_ice    ! shortwave radiation over the ice/snow
900      !!
901      INTEGER, PARAMETER  ::   jp24 = 24   ! sampling of the daylight period (sunrise to sunset) into 24 equal parts
902      !!
903      INTEGER  ::   ji, jj, jl, jt    ! dummy loop indices
904      INTEGER  ::   ijpl              ! number of ice categories (3rd dim of pqsr_ice)
905      INTEGER  ::   indaet            !  = -1, 0, 1 for odd, normal and leap years resp.
906      INTEGER  ::   iday              ! integer part of day
907      !!
908      REAL(wp) ::   zcmue, zcmue2, ztamr          ! temporary scalars
909      REAL(wp) ::   zmt1, zmt2, zmt3              !    -         -
910      REAL(wp) ::   zdecl, zsdecl, zcdecl         !    -         -
911      REAL(wp) ::   zlha, zdaycor, zes            !    -         -
912      REAL(wp) ::   zxday, zdist, zcoef, zcoef1   !    -         -
913      REAL(wp) ::   zqsr_ice_cs, zqsr_ice_os      !    -         -
914
915      REAL(wp), DIMENSION(:,:), POINTER ::   zev                      ! vapour pressure
916      REAL(wp), DIMENSION(:,:), POINTER ::   zdlha, zlsrise, zlsset   ! 2D workspace
917      REAL(wp), DIMENSION(:,:), POINTER ::   zps, zpc   ! sine (cosine) of latitude per sine (cosine) of solar declination
918      !!---------------------------------------------------------------------
919      !
920      IF( nn_timing == 1 )  CALL timing_start('blk_clio_qsr_ice')
921      !
922      CALL wrk_alloc( jpi,jpj, zev, zdlha, zlsrise, zlsset, zps, zpc )
923
924      ijpl = SIZE(pqsr_ice, 3 )      ! number of ice categories
925     
926      ! Saturated water vapour and vapour pressure
927      ! ------------------------------------------
928!CDIR NOVERRCHK
929!CDIR COLLAPSE
930      DO jj = 1, jpj
931!CDIR NOVERRCHK
932         DO ji = 1, jpi           
933            ztamr = sf(jp_tair)%fnow(ji,jj,1) - rtt           
934            zmt1  = SIGN( 17.269,  ztamr )
935            zmt2  = SIGN( 21.875,  ztamr )
936            zmt3  = SIGN( 28.200, -ztamr )
937            zes = 611.0 * EXP(  ABS( ztamr ) * MIN ( zmt1, zmt2 )   &              ! Saturation water vapour
938               &                     / ( sf(jp_tair)%fnow(ji,jj,1) - 35.86  + MAX( 0.e0, zmt3 ) )  )
939            zev(ji,jj) = sf(jp_humi)%fnow(ji,jj,1) * zes * 1.0e-05                 ! vapour pressure 
940         END DO
941      END DO
942
943      !-----------------------------------!
944      !  Computation of solar irradiance  !
945      !-----------------------------------!
946!!gm : hard coded  leap year ???
947      indaet   = 1                                    ! = -1, 0, 1 for odd, normal and leap years resp.
948      zxday = nday_year + rdtbs2 / rday               ! day of the year at which the fluxes are calculated
949      iday  = INT( zxday )                            ! (centred at the middle of the ice time step)
950      CALL flx_blk_declin( indaet, iday, zdecl )      ! solar declination of the current day
951      zsdecl = SIN( zdecl * rad )                     ! its sine
952      zcdecl = COS( zdecl * rad )                     ! its cosine
953
954     
955      !  correction factor added for computation of shortwave flux to take into account the variation of
956      !  the distance between the sun and the earth during the year (Oberhuber 1988)
957      zdist    = zxday * 2. * rpi / REAL(nyear_len(1), wp)
958      zdaycor  = 1.0 + 0.0013 * SIN( zdist ) + 0.0342 * COS( zdist )
959
960!CDIR NOVERRCHK
961      DO jj = 1, jpj
962!CDIR NOVERRCHK
963         DO ji = 1, jpi
964            !  product of sine (cosine) of latitude and sine (cosine) of solar declination
965            zps(ji,jj) = SIN( gphit(ji,jj) * rad ) * zsdecl
966            zpc(ji,jj) = COS( gphit(ji,jj) * rad ) * zcdecl
967            !  computation of the both local time of sunrise and sunset
968            zlsrise(ji,jj) = ACOS( - SIGN( 1.e0, zps(ji,jj) )    &
969               &                   * MIN(  1.e0, SIGN( 1.e0, zps(ji,jj) ) * ( zps(ji,jj) / zpc(ji,jj) )  )   ) 
970            zlsset (ji,jj) = - zlsrise(ji,jj)
971            !  dividing the solar day into jp24 segments of length zdlha
972            zdlha  (ji,jj) = ( zlsrise(ji,jj) - zlsset(ji,jj) ) / REAL( jp24, wp )
973         END DO
974      END DO
975
976
977      !---------------------------------------------!
978      !  shortwave radiation absorbed by the ice    !
979      !---------------------------------------------!
980      ! compute and sum ice qsr over the daylight for each ice categories
981      pqsr_ice(:,:,:) = 0.e0
982      zcoef1 = zdaycor / ( 2. * rpi )       ! Correction for the ellipsity of the earth orbit
983     
984      !                    !----------------------------!
985      DO jl = 1, ijpl      !  loop over ice categories  !
986         !                 !----------------------------!
987!CDIR NOVERRCHK   
988         DO jt = 1, jp24   
989            zcoef = FLOAT( jt ) - 0.5
990!CDIR NOVERRCHK     
991!CDIR COLLAPSE
992            DO jj = 1, jpj
993!CDIR NOVERRCHK
994               DO ji = 1, jpi
995                  zlha = COS(  zlsrise(ji,jj) - zcoef * zdlha(ji,jj)  )                  ! local hour angle
996                  zcmue              = MAX( 0.e0 ,   zps(ji,jj) + zpc(ji,jj) * zlha  )   ! cos of local solar altitude
997                  zcmue2             = 1368.0 * zcmue * zcmue
998                 
999                  !  solar heat flux absorbed by the ice/snow system (Shine and Crane 1984 adapted to high albedo)
1000                  zqsr_ice_cs =  ( 1.0 - pa_ice_cs(ji,jj,jl) ) * zdlha(ji,jj) * zcmue2        &   ! clear sky
1001                     &        / ( ( 1.0 + zcmue ) * zev(ji,jj) + 1.2 * zcmue + 0.0455 )
1002                  zqsr_ice_os = zdlha(ji,jj) * SQRT( zcmue )                                  &   ! overcast sky
1003                     &        * ( 53.5 + 1274.5 * zcmue )      * ( 1.0 - 0.996  * pa_ice_os(ji,jj,jl) )    &
1004                     &        / (  1.0 + 0.139  * stauc(ji,jj) * ( 1.0 - 0.9435 * pa_ice_os(ji,jj,jl) ) )       
1005             
1006                  pqsr_ice(ji,jj,jl) = pqsr_ice(ji,jj,jl) + (  ( 1.0 - sf(jp_ccov)%fnow(ji,jj,1) ) * zqsr_ice_cs    &
1007                     &                                       +         sf(jp_ccov)%fnow(ji,jj,1)   * zqsr_ice_os  )
1008               END DO
1009            END DO
1010         END DO
1011         !
1012         ! Correction : Taking into account the ellipsity of the earth orbit
1013         pqsr_ice(:,:,jl) = pqsr_ice(:,:,jl) * zcoef1 * tmask(:,:,1)
1014         !
1015         !                 !--------------------------------!
1016      END DO               !  end loop over ice categories  !
1017      !                    !--------------------------------!
1018
1019
1020!!gm  : this should be suppress as input data have been passed through lbc_lnk
1021      DO jl = 1, ijpl
1022         CALL lbc_lnk( pqsr_ice(:,:,jl) , 'T', 1. )
1023      END DO
1024      !
1025      CALL wrk_dealloc( jpi,jpj, zev, zdlha, zlsrise, zlsset, zps, zpc )
1026      !
1027      IF( nn_timing == 1 )  CALL timing_stop('blk_clio_qsr_ice')
1028      !
1029   END SUBROUTINE blk_clio_qsr_ice
1030
1031
1032   SUBROUTINE flx_blk_declin( ky, kday, pdecl )
1033      !!---------------------------------------------------------------------------
1034      !!               ***  ROUTINE flx_blk_declin  ***
1035      !!         
1036      !! ** Purpose :   Computation of the solar declination for the day
1037      !!       
1038      !! ** Method  :   ???
1039      !!---------------------------------------------------------------------
1040      INTEGER , INTENT(in   ) ::   ky      ! = -1, 0, 1 for odd, normal and leap years resp.
1041      INTEGER , INTENT(in   ) ::   kday    ! day of the year ( kday = 1 on january 1)
1042      REAL(wp), INTENT(  out) ::   pdecl   ! solar declination
1043      !!
1044      REAL(wp) ::   a0  =  0.39507671      ! coefficients for solar declinaison computation
1045      REAL(wp) ::   a1  = 22.85684301      !     "              ""                 "
1046      REAL(wp) ::   a2  = -0.38637317      !     "              ""                 "
1047      REAL(wp) ::   a3  =  0.15096535      !     "              ""                 "
1048      REAL(wp) ::   a4  = -0.00961411      !     "              ""                 "
1049      REAL(wp) ::   b1  = -4.29692073      !     "              ""                 "
1050      REAL(wp) ::   b2  =  0.05702074      !     "              ""                 "
1051      REAL(wp) ::   b3  = -0.09028607      !     "              ""                 "
1052      REAL(wp) ::   b4  =  0.00592797
1053      !!
1054      REAL(wp) ::   zday   ! corresponding day of type year (cf. ky)
1055      REAL(wp) ::   zp     ! temporary scalars
1056      !!---------------------------------------------------------------------
1057           
1058      IF    ( ky == 1 )  THEN   ;   zday = REAL( kday, wp ) - 0.5
1059      ELSEIF( ky == 3 )  THEN   ;   zday = REAL( kday, wp ) - 1.
1060      ELSE                      ;   zday = REAL( kday, wp )
1061      ENDIF
1062     
1063      zp = rpi * ( 2.0 * zday - 367.0 ) / REAL(nyear_len(1), wp)
1064     
1065      pdecl  = a0                                                                      &
1066         &   + a1 * COS( zp ) + a2 * COS( 2. * zp ) + a3 * COS( 3. * zp ) + a4 * COS( 4. * zp )   &
1067         &   + b1 * SIN( zp ) + b2 * SIN( 2. * zp ) + b3 * SIN( 3. * zp ) + b4 * SIN( 4. * zp )
1068      !
1069   END SUBROUTINE flx_blk_declin
1070
1071   !!======================================================================
1072END MODULE sbcblk_clio
Note: See TracBrowser for help on using the repository browser.