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/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90 @ 5951

Last change on this file since 5951 was 5951, checked in by timgraham, 8 years ago

Merged trunk r5936 into branch

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