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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90 @ 2674

Last change on this file since 2674 was 2674, checked in by rblod, 12 years ago

First set of changes in OPA_SRC to ensure AGRIF compatibility

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