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.
Changeset 840 – NEMO

Changeset 840


Ignore:
Timestamp:
2008-03-11T14:51:35+01:00 (16 years ago)
Author:
ctlod
Message:

dev_001_SBC : add the good interface for CLIO bulk under the new surface module

Location:
branches/dev_001_SBC/NEMO/OPA_SRC/SBC
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/dev_001_SBC/NEMO/OPA_SRC/SBC/sbcblk_clio.F90

    r710 r840  
    11MODULE sbcblk_clio 
    22   !!====================================================================== 
    3    !!                       ***  MODULE  sbcblk_clio  *** 
    4    !! Ocean forcing:  momentum, heat and freshwater flux formulation 
     3   !!                   ***  MODULE  sbcblk_clio  *** 
     4   !! Ocean forcing:  bulk thermohaline forcing of the ocean (or ice) 
    55   !!===================================================================== 
    6    !! History :  8.0  !  01-04  (Louvain-La-Neuve)  Original code 
    7    !!            8.5  !  02-09  (C. Ethe , G. Madec )  F90: Free form and module 
    8    !!            9.0  !  06-06  (G. Madec)  surface module 
     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 
    910   !!---------------------------------------------------------------------- 
    10     
    11    !!---------------------------------------------------------------------- 
    12    !!   sbc_blk_clio  : bulk formulation as ocean surface boundary condition 
    13    !!                   (forced mode, CORE bulk formulea) 
    14    !!   blk_oce_clio  : ocean: computes momentum, heat and freshwater fluxes 
    15    !!   blk_ice_clio  : ice  : computes momentum, heat and freshwater fluxes 
    16    !!---------------------------------------------------------------------- 
    17    !!   flx_blk_declin : Computation of the solar declination 
     11   !!   sbc_blk_clio   : CLIO bulk formulation: read and update required input fields 
     12   !!   blk_clio_oce   : ocean CLIO bulk formulea: compute momentum, heat and freswater fluxes for the ocean 
     13   !!   blk_ice_clio   : ice   CLIO bulk formulea: compute momentum, heat and freswater fluxes for the sea-ice 
     14   !!   blk_clio_qsr_oce : shortwave radiation for ocean computed from the cloud cover 
     15   !!   blk_clio_qsr_ice : shortwave radiation for ice   computed from the cloud cover 
     16   !!   flx_blk_declin : solar declinaison 
    1817   !!---------------------------------------------------------------------- 
    1918   USE oce             ! ocean dynamics and tracers 
    2019   USE dom_oce         ! ocean space and time domain 
    21    USE cpl_oce         ! ??? 
    2220   USE phycst          ! physical constants 
    23    USE daymod 
    24  
     21   USE daymod          ! calendar 
     22   USE ocfzpt          ! ocean freezing point 
     23   USE fldread         ! read input fields 
    2524   USE sbc_oce         ! Surface boundary condition: ocean fields 
    26    USE sbc_ice         ! Surface boundary condition: ocean fields 
    27  
    28    USE fldread         ! read input fields 
    29  
    30    USE ocfzpt          ! ??? 
    31  
    32    USE iom 
    33    USE in_out_manager 
    34    USE lbclnk 
     25   USE iom             ! I/O manager library 
     26   USE in_out_manager  ! I/O manager 
     27   USE lib_mpp         ! distribued memory computing library 
     28   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     29 
     30   USE albedo 
    3531   USE prtctl          ! Print control 
    36  
     32#if defined key_lim3 
     33   USE par_ice 
     34#elif defined key_lim2 
     35   USE ice_2 
     36#endif 
    3737   IMPLICIT NONE 
    3838   PRIVATE 
    3939 
    40    PUBLIC   sbc_blk_clio   ! routine called in sbcmod 
    41    PUBLIC   blk_ice_clio   ! routine called in sbcice_lim module 
    42  
    43    INTEGER , PARAMETER ::    & 
    44       jpfld   = 7,     &  ! number of files to read 
    45       jp_utau = 1,     &  ! index of wind stress (i-component)  (m/s)    at U-point 
    46       jp_vtau = 2,     &  ! index of wind stress (j-component)  (m/s)    at V-point 
    47       jp_wspd = 3,     &  ! index of XXm wind module            (m/s)    at T-point 
    48       jp_tair = 4,     &  ! index of  2m air temperature        (Celcius) 
    49       jp_humi = 5,     &  ! index of specific humidity          ( % ) 
    50       jp_cldc = 6,     &  ! Cloud cover                       ( % ) 
    51       jp_prec = 7         ! index of total precipitation        (kg/m2/s) 
    52    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (informations on files, fields read)  
    53  
    54    !! * CLIO bulk parameters 
     40   PUBLIC sbc_blk_clio        ! routine called by flx.F90  
     41   PUBLIC blk_ice_clio        ! routine called by flx.F90  
     42 
     43   INTEGER , PARAMETER ::   jpfld   = 7           ! maximum number of files to read  
     44   INTEGER , PARAMETER ::   jp_wndi = 1           ! index of 10m wind velocity (i-component) (m/s)    at U-point 
     45   INTEGER , PARAMETER ::   jp_wndj = 2           ! index of 10m wind velocity (j-component) (m/s)    at V-point 
     46   INTEGER , PARAMETER ::   jp_wndm = 3           ! index of 10m wind module                 (m/s)    at T-point 
     47   INTEGER , PARAMETER ::   jp_humi = 4           ! index of specific humidity               ( % ) 
     48   INTEGER , PARAMETER ::   jp_ccov = 5           ! index of cloud cover                     ( % ) 
     49   INTEGER , PARAMETER ::   jp_tair = 6           ! index of 10m air temperature             (Kelvin) 
     50   INTEGER , PARAMETER ::   jp_prec = 7           ! index of total precipitation (rain+snow) (Kg/m2/s) 
     51   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read) 
     52 
     53 
     54   !!  
     55!!!!gm  to be moved 
     56   INTEGER, PARAMETER  ::   jpl = 1          ! number of layer in the ice   
     57!!!!gm  to be moved 
     58 
     59 
    5560   INTEGER, PARAMETER  ::   jpintsr = 24          ! number of time step between sunrise and sunset 
    56       !                                           ! uses for heat flux computation 
    57    LOGICAL             ::   lbulk_init = .TRUE.   ! flag, bulk initialization done or not) 
    58        
    59    REAL(wp), DIMENSION(jpi,jpj) ::   stauc        ! cloud optical depth  
    60    REAL(wp), DIMENSION(jpi,jpj) ::   sbudyko      ! ??? 
    61        
    62    !! * constants for bulk computation (flx_blk)  
    63    REAL(wp), DIMENSION(19)  ::  budyko            ! BUDYKO's coefficient 
    64    !                                              ! BUDYKO's coefficient (cloudiness effect on LW radiation): 
    65    DATA budyko / 1.00, 0.98, 0.95, 0.92, 0.89, 0.86, 0.83, 0.80, 0.78, 0.75,  & 
     61   !                                              ! uses for heat flux computation 
     62   LOGICAL ::   lbulk_init = .TRUE.   ! flag, bulk initialization done or not) 
     63 
     64   REAL(wp) ::   cai = 1.40e-3 ! best estimate of atm drag in order to get correct FS export in ORCA2-LIM 
     65   REAL(wp) ::   cao = 1.00e-3 ! chosen by default  ==> should depends on many things...  !!gmto be updated 
     66 
     67   REAL(wp) ::   yearday     !: number of days per year    
     68   REAL(wp) ::   rdtbs2      !: number of days per year    
     69    
     70   REAL(wp), DIMENSION(19)  ::  budyko            ! BUDYKO's coefficient (cloudiness effect on LW radiation) 
     71   DATA budyko / 1.00, 0.98, 0.95, 0.92, 0.89, 0.86, 0.83, 0.80, 0.78, 0.75,   & 
    6672      &          0.72, 0.69, 0.67, 0.64, 0.61, 0.58, 0.56, 0.53, 0.50 / 
    67    REAL(wp), DIMENSION(20)  ::   tauco            ! cloud optical depth coefficient 
    68    !                                              ! Cloud optical depth coefficient 
     73   REAL(wp), DIMENSION(20)  :: tauco              ! cloud optical depth coefficient 
    6974   DATA tauco / 6.6, 6.6, 7.0, 7.2, 7.1, 6.8, 6.5, 6.6, 7.1, 7.6,   & 
    7075      &         6.6, 6.1, 5.6, 5.5, 5.8, 5.8, 5.6, 5.6, 5.6, 5.6 / 
    71  
    72    REAL(wp) ::   zeps    = 1.e-20  ,  &  ! constant values 
    73       &          zeps0   = 1.e-13  ,  & 
    74       &          zeps1   = 1.e-06  ,  & 
    75       &          zzero   = 0.e0    ,  & 
    76       &          zone    = 1.0 
    77  
    78    REAL(wp) ::   yearday   !  
    79  
     76   !! 
     77   REAL(wp), DIMENSION(jpi,jpj) ::   sbudyko      ! cloudiness effect on LW radiation 
     78   REAL(wp), DIMENSION(jpi,jpj) ::   stauc        ! cloud optical depth  
     79    
     80   REAL(wp)  ::   zeps    = 1.e-20                ! constant values 
     81   REAL(wp)  ::   zeps0   = 1.e-13   
     82 
     83#  include "vectopt_loop_substitute.h90" 
    8084   !!---------------------------------------------------------------------- 
    81    !!   OPA 9.0 , LOCEAN-IPSL (2006) 
    82    !! $Header: $ 
     85   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)  
     86   !! $Id:$  
    8387   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    8488   !!---------------------------------------------------------------------- 
     
    8690CONTAINS 
    8791 
     92 
    8893   SUBROUTINE sbc_blk_clio( kt ) 
    8994      !!--------------------------------------------------------------------- 
    90       !!                    ***  ROUTINE sbc_blk_clio  *** 
    91       !! 
     95      !!                    ***  ROUTINE sbc_blk_core  *** 
     96      !!                    
    9297      !! ** Purpose :   provide at each time step the surface ocean fluxes 
    93       !!      (momentum, heat, freshwater and runoff) 
     98      !!      (momentum, heat, freshwater and runoff)  
    9499      !! 
    95100      !! ** Method  :   READ each fluxes in NetCDF files 
     
    118123      INTEGER, INTENT( in  ) ::   kt   ! ocean time step 
    119124      !! 
    120       INTEGER  ::   jf                  ! dummy indices 
     125      INTEGER  ::   jf       ! dummy indices 
    121126      INTEGER  ::   ierror   ! return error code 
    122127      !! 
    123       CHARACTER(len=100) ::  cn_dir   !   Root directory for location of clio files 
    124       TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i         ! array of namelist informations on the fields to read 
    125       TYPE(FLD_N) ::   sn_utau, sn_vtau, sn_wspd, sn_tair,  &  ! informations about the fields to be read 
    126          &             sn_humi, sn_cldc, sn_prec 
    127       NAMELIST/namsbc_clio/ cn_dir, sn_utau, sn_vtau, sn_wspd, sn_tair,   & 
    128          &                          sn_humi, sn_cldc, sn_prec 
     128      CHARACTER(len=100) ::  cn_dir                            !   Root directory for location of CLIO files 
     129      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read 
     130      TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_wndm, sn_tair      ! informations about the fields to be read 
     131      TYPE(FLD_N) ::   sn_humi, sn_ccov, sn_prec               !   "                                 " 
     132      !! 
     133      NAMELIST/namsbc_core/ cn_dir, sn_wndi, sn_wndj, sn_wndm, sn_humi,   & 
     134         &                          sn_ccov, sn_tair, sn_prec 
    129135      !!--------------------------------------------------------------------- 
    130136 
     
    134140         ! set file information (default values) 
    135141         cn_dir = './'       ! directory in which the model is executed 
     142 
    136143         ! (NB: frequency positive => hours, negative => months) 
    137          !              !  file  ! frequency !  variable  ! time intep !  clim  ! starting ! 
    138          !              !  name  !  (hours)  !   name     !   (T/F)    !  (0/1) !  record  ! 
    139          sn_utau = FLD_N( 'utau' ,    24.    , 'utau'     ,  .FALSE.   ,    0   ,    0     ) 
    140          sn_vtau = FLD_N( 'vtau' ,    24.    , 'vtau'     ,  .FALSE.   ,    0   ,    0     ) 
    141          sn_wspd = FLD_N( 'wspd' ,    24.    , 'wspd'     ,  .FALSE.   ,    0   ,    0     ) 
    142          sn_tair = FLD_N( 'tair' ,    24.    , 'Tair'     ,  .FALSE.   ,    0   ,    0     ) 
    143          sn_humi = FLD_N( 'humi' ,   -12.    , 'humi'     ,  .FALSE.   ,    0   ,    0     ) 
    144          sn_cldc = FLD_N( 'cloud',   -12.    , 'cloud'    ,  .FALSE.   ,    0   ,    0     ) 
    145          sn_prec = FLD_N( 'rain' ,   -12.    , 'precip'   ,  .FALSE.   ,    0   ,    0     ) 
    146  
    147          REWIND ( numnam )                    ! ... read in namlist namsbc_clio 
    148          READ   ( numnam, namsbc_clio ) 
     144         !            !    file     ! frequency !  variable  ! time intep !  clim  ! starting ! 
     145         !            !    name     !  (hours)  !   name     !   (T/F)    !  (0/1) !  record  ! 
     146         sn_wndi = FLD_N( 'uwnd10m' ,    24.    ,  'u_10'    ,  .true.    ,    0   ,     0    )  
     147         sn_wndj = FLD_N( 'vwnd10m' ,    24.    ,  'v_10'    ,  .true.    ,    0   ,     0    )  
     148         sn_wndm = FLD_N( 'mwnd10m' ,    24.    ,  'm_10'    ,  .true.    ,    0   ,     0    )  
     149         sn_tair = FLD_N( 'tair10m' ,    24.    ,  't_10'    ,  .FALSE.   ,    0   ,     0    )  
     150         sn_humi = FLD_N( 'humi10m' ,    24.    ,  'q_10'    ,  .FALSE.   ,    0   ,     0    )  
     151         sn_ccov = FLD_N( 'ccover'  ,   -12.    ,  'cloud'   ,  .TRUE.    ,    0   ,     0    )  
     152         sn_prec = FLD_N( 'precip'  ,   -12.    ,  'precip'  ,  .TRUE.    ,    0   ,     0    )  
     153 
     154         REWIND( numnam )                    ! ... read in namlist namsbc_core 
     155         READ  ( numnam, namsbc_core ) 
    149156 
    150157         ! store namelist information in an array 
    151          slf_i(jp_utau) = sn_utau   ;   slf_i(jp_vtau) = sn_vtau 
    152          slf_i(jp_wspd) = sn_wspd   ;   slf_i(jp_tair) = sn_tair 
    153          slf_i(jp_humi) = sn_humi   ;   slf_i(jp_cldc) = sn_cldc 
    154          slf_i(jp_prec) = sn_prec 
    155  
     158         slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj   ;   slf_i(jp_wndm) = sn_wndm 
     159         slf_i(jp_tair) = sn_tair   ;   slf_i(jp_humi) = sn_humi 
     160         slf_i(jp_ccov) = sn_ccov   ;   slf_i(jp_prec) = sn_prec 
     161          
    156162         ! set sf structure 
    157163         ALLOCATE( sf(jpfld), STAT=ierror ) 
    158164         IF( ierror > 0 ) THEN 
    159             CALL ctl_stop( 'sbc_blk_clio: unable to allocate sf_sst structure' )   ;   RETURN 
     165            CALL ctl_stop( 'sbc_blk_clio: unable to allocate sf structure' )   ;   RETURN 
    160166         ENDIF 
    161  
     167         ! 
    162168         DO jf = 1, jpfld 
    163169            WRITE(sf(jf)%clrootname,'(a,a)' )   TRIM( cn_dir ), TRIM( slf_i(jf)%clname ) 
     
    170176 
    171177         IF(lwp) THEN      ! control print 
    172             WRITE(numout,*) 
    173             WRITE(numout,*) 'sbc_blk_clio : CLIO bulk formulation for ocean surface boundary condition' 
     178            WRITE(numout,*)             
     179            WRITE(numout,*) 'sbc_blk_clio : flux formulattion for ocean surface boundary condition' 
    174180            WRITE(numout,*) '~~~~~~~~~~~~ ' 
    175             WRITE(numout,*) '          namsbc_clio Namelist' 
     181            WRITE(numout,*) '          namsbc_core Namelist' 
    176182            WRITE(numout,*) '          list of files and frequency (>0: in hours ; <0 in months)' 
    177183            DO jf = 1, jpfld 
    178                 WRITE(numout,*) '               root filename: '  , trim( sf(jf)%clrootname ),   & 
    179                    &                          ' variable name: '  , trim( sf(jf)%clvar      ) 
     184                WRITE(numout,*) '               file root name: ' , TRIM( sf(jf)%clrootname ),   & 
     185                   &                          ' variable name: '  , TRIM( sf(jf)%clvar      ) 
    180186                WRITE(numout,*) '               frequency: '      ,       sf(jf)%freqh       ,   & 
    181187                   &                          ' time interp: '    ,       sf(jf)%ln_tint     ,   & 
     
    186192         ! 
    187193      ENDIF 
    188  
    189       CALL fld_read( kt, nn_fsbc, sf )           ! Read input fields and provides the 
    190       !                                          ! input fields at the current time-step 
    191 !CDIR COLLAPSE 
    192       utau(:,:) = sf(jp_utau)%fnow(:,:)          ! set surface ocean stresses directly  
    193 !CDIR COLLAPSE 
    194       vtau(:,:) = sf(jp_vtau)%fnow(:,:)          ! from the input values 
    195  
    196       CALL blk_oce_clio( )                       ! set the ocean surface heat and freshwater fluxes 
    197       !                                          ! using CLIO bulk formulea 
    198  
    199 ! temporary staff : set fluxes to zero.... 
    200       qns (:,:)= 0.e0 
    201       qsr (:,:)= 0.e0 
    202       emp (:,:)= 0.e0 
    203       emps(:,:)= 0.e0 
    204  
    205       ! control print (if less than 100 time-step asked) 
    206 !!!   IF( nitend-nit000 <= 100 .AND. lwp ) THEN 
    207       IF( kt == nit000 .AND. lwp ) THEN 
    208          WRITE(numout,*) 
    209          WRITE(numout,*) '        CLIO bulk fields at nit000' 
    210          DO jf = 1, jpfld 
    211             WRITE(numout,*) 
    212             WRITE(numout,*) ' day: ', ndastp , TRIM(sf(jf)%clvar) 
    213             CALL prihre( sf(jf)%fnow, jpi, jpj, 1, jpi, 20, 1, jpj, 10, 1., numout ) 
    214          END DO 
    215          CALL FLUSH(numout) 
     194      !                                         ! ====================== ! 
     195      !                                         !    At each time-step   ! 
     196      !                                         ! ====================== ! 
     197      ! 
     198      CALL fld_read( kt, nn_fsbc, sf )                ! input fields provided at the current time-step 
     199      ! 
     200      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 
     201          CALL blk_oce_clio( sst_m, ssu_m, ssv_m )    ! compute the surface ocean fluxes using CLIO bulk formulea 
     202      ENDIF                                           !  
     203      ! 
     204   END SUBROUTINE sbc_blk_clio 
     205 
     206 
     207   SUBROUTINE blk_oce_clio( pst, pu, pv ) 
     208      !!--------------------------------------------------------------------------- 
     209      !!                     ***  ROUTINE blk_oce_clio  *** 
     210      !!                  
     211      !!  ** Purpose :   Compute momentum, heat and freshwater fluxes at ocean surface 
     212      !!               using CLIO bulk formulea 
     213      !!          
     214      !!  ** Method  :   The flux of heat at the ocean surfaces are derived 
     215      !!       from semi-empirical ( or bulk ) formulae which relate the flux to  
     216      !!       the properties of the surface and of the lower atmosphere. Here, we 
     217      !!       follow the work of Oberhuber, 1988    
     218      !!               - momentum flux (stresses) directly read in files at U- and V-points 
     219      !!               - compute ocean and ice albedos (call flx_blk_albedo)   
     220      !!               - compute shortwave radiation for ocean (call blk_clio_qsr_oce) 
     221      !!               - compute long-wave radiation for the ocean 
     222      !!               - compute the turbulent heat fluxes over the ocean 
     223      !!               - deduce the evaporation over the ocean 
     224      !!  ** Action  :   Fluxes over the ocean: 
     225      !!               - utau, vtau  i- and j-component of the wind stress 
     226      !!               - qns, qsr    non-slor and solar heat flux 
     227      !!               - emp, emps   evaporation minus precipitation 
     228      !!---------------------------------------------------------------------- 
     229      REAL(wp), INTENT(in), DIMENSION(jpi,jpj) ::   pst   ! surface temperature                      [Celcius] 
     230      REAL(wp), INTENT(in), DIMENSION(jpi,jpj) ::   pu    ! surface current at U-point (i-component) [m/s] 
     231      REAL(wp), INTENT(in), DIMENSION(jpi,jpj) ::   pv    ! surface current at V-point (j-component) [m/s] 
     232      !! 
     233      INTEGER  ::   ji, jj   ! dummy loop indices 
     234      !! 
     235      REAL(wp) ::   zrhova, zcsho, zcleo, zcldeff               ! temporary scalars 
     236      REAL(wp) ::   zqsato, zdteta, zdeltaq, ztvmoy, zobouks    !    -         - 
     237      REAL(wp) ::   zpsims, zpsihs, zpsils, zobouku, zxins, zpsimu   !    -         - 
     238      REAL(wp) ::   zpsihu, zpsilu, zstab,zpsim, zpsih, zpsil   !    -         - 
     239      REAL(wp) ::   zvatmg, zcmn, zchn, zcln, zcmcmn, zdenum    !    -         - 
     240      REAL(wp) ::   zdtetar, ztvmoyr, zlxins, zchcm, zclcm      !    -         - 
     241      REAL(wp) ::   zmt1, zmt2, zmt3, ztatm3, ztamr, ztaevbk    !    -         - 
     242      REAL(wp) ::   zcoef, zsst, ztatm, zcco1, zpatm            !    -         - 
     243      REAL(wp) ::   zrhoa, zev, zes, zeso, zqatm, zevsqr        !    -         - 
     244      !! 
     245      REAL(wp), DIMENSION(jpi,jpj) ::   zqlw        ! long-wave heat flux over ocean 
     246      REAL(wp), DIMENSION(jpi,jpj) ::   zqla        ! latent heat flux over ocean 
     247      REAL(wp), DIMENSION(jpi,jpj) ::   zqsb        ! sensible heat flux over ocean 
     248      !!--------------------------------------------------------------------- 
     249 
     250      zpatm = 101000.      ! atmospheric pressure  (assumed constant here) 
     251 
     252      !------------------------------------! 
     253      !   momentum fluxes  (utau, vtau )   ! 
     254      !------------------------------------! 
     255 
     256      ! Change from wind speed to wind stress over OCEAN (cao is used) 
     257      ! note: 1.3 = air density) 
     258!!gm  the module is wrong as U and V points are not the same 
     259!!gm  the surface currents are not taken into account... 
     260!CDIR COLLAPSE 
     261      DO jj = 1 , jpj 
     262         DO ji = 1, jpi 
     263            zcoef = 1.3 * cao * SQRT(  sf(jp_wndi)%fnow(ji,jj)*sf(jp_wndi)%fnow(ji,jj)   & 
     264               &                     + sf(jp_wndj)%fnow(ji,jj)*sf(jp_wndj)%fnow(ji,jj) ) 
     265            utau(ji,jj) = 1.3 * cao * zcoef * sf(jp_wndi)%fnow(ji,jj) 
     266            vtau(ji,jj) = 1.3 * cao * zcoef * sf(jp_wndj)%fnow(ji,jj) 
     267         END DO 
     268      END DO 
     269!!gm  better coding:  use the wind module and averaging to U and V points 
     270!     zcoef = 1.3 * cao * 0.5 
     271!     DO jj = 1 , jpjm1 
     272!        DO ji = 1, fs_jpim1 
     273!           utau(ji,jj) = zcoef * ( sf(jp_wndm)%fnow(ji+1,jj) + sf(jp_wndm)%fnow(ji,jj) ) * sf(jp_wndi)%fnow(ji,jj) 
     274!           vtau(ji,jj) = zcoef * ( sf(jp_wndm)%fnow(ji,jj+1) + sf(jp_wndm)%fnow(ji,jj) ) * sf(jp_wndj)%fnow(ji,jj) 
     275!          END DO 
     276!       END DO 
     277!      CALL lbc_lnk( utau(:,:), 'U', -1. ) 
     278!      CALL lbc_lnk( vtau(:,:), 'V', -1. ) 
     279!!gm end coding 
     280 
     281       
     282      !------------------------------------------------! 
     283      !   Shortwave radiation for ocean and snow/ice   ! 
     284      !------------------------------------------------! 
     285       
     286      CALL blk_clio_qsr_oce( qsr ) 
     287 
     288 
     289      !------------------------! 
     290      !   Other ocean fluxes   ! 
     291      !------------------------! 
     292!CDIR NOVERRCHK 
     293!CDIR COLLAPSE 
     294      DO jj = 1, jpj 
     295!CDIR NOVERRCHK 
     296         DO ji = 1, jpi 
     297            ! 
     298            zsst  = pst(ji,jj)              + rt0           ! converte Celcius to Kelvin the SST and Tair 
     299            ztatm = sf(jp_tair)%fnow(ji,jj) + rt0           ! and set minimum value far above 0 K (=rt0 over land) 
     300            zcco1 = 1.0 - sf(jp_ccov)%fnow(ji,jj)           ! fraction of clear sky ( 1 - cloud cover) 
     301            zrhoa = zpatm / ( 287.04 * ztatm )              ! air density (equation of state for dry air)  
     302            ztamr = ztatm - rtt                             ! Saturation water vapour 
     303            zmt1  = SIGN( 17.269,  ztamr )                  !           || 
     304            zmt2  = SIGN( 21.875,  ztamr )                  !          \  / 
     305            zmt3  = SIGN( 28.200, -ztamr )                  !           \/ 
     306            zes   = 611.0 * EXP(  ABS( ztamr ) * MIN ( zmt1, zmt2 ) / ( ztatm - 35.86  + MAX( 0.e0, zmt3 ) )  ) 
     307            zev    = sf(jp_humi)%fnow(ji,jj) * zes          ! vapour pressure   
     308            zevsqr = SQRT( zev * 0.01 )                     ! square-root of vapour pressure 
     309            zqatm = 0.622 * zev / ( zpatm - 0.378 * zev )   ! specific humidity  
     310 
     311            !--------------------------------------! 
     312            !  long-wave radiation over the ocean  !  ( Berliand 1952 ; all latitudes ) 
     313            !--------------------------------------! 
     314            ztatm3  = ztatm * ztatm * ztatm 
     315            zcldeff = 1.0 - sbudyko(ji,jj) * sf(jp_ccov)%fnow(ji,jj) * sf(jp_ccov)%fnow(ji,jj)     
     316            ztaevbk = ztatm * ztatm3 * zcldeff * ( 0.39 - 0.05 * zevsqr )  
     317            ! 
     318            zqlw(ji,jj) = - emic * stefan * ( ztaevbk + 4. * ztatm3 * ( zsst - ztatm ) )  
     319 
     320            !-------------------------------------------------- 
     321            !  Latent and sensible heat fluxes over the ocean  
     322            !-------------------------------------------------- 
     323            !                                                          ! vapour pressure at saturation of ocean 
     324            zeso =  611.0 * EXP ( 17.2693884 * ( zsst - rtt ) * tmask(ji,jj,1) / ( zsst - 35.86 ) ) 
     325 
     326            zqsato = ( 0.622 * zeso ) / ( zpatm - 0.378 * zeso )       ! humidity close to the ocean surface (at saturation) 
     327 
     328            ! Drag coefficients from Large and Pond (1981,1982) 
     329            !                                                          ! Stability parameters 
     330            zdteta  = zsst - ztatm 
     331            zdeltaq = zqatm - zqsato 
     332            ztvmoy  = ztatm * ( 1. + 2.2e-3 * ztatm * zqatm ) 
     333            zdenum  = MAX( sf(jp_wndm)%fnow(ji,jj) * sf(jp_wndm)%fnow(ji,jj) * ztvmoy, zeps ) 
     334            zdtetar = zdteta / zdenum 
     335            ztvmoyr = ztvmoy * ztvmoy * zdeltaq / zdenum 
     336            !                                                          ! case of stable atmospheric conditions 
     337            zobouks = -70.0 * 10. * ( zdtetar + 3.2e-3 * ztvmoyr ) 
     338            zobouks = MAX( 0.e0, zobouks ) 
     339            zpsims = -7.0 * zobouks 
     340            zpsihs =  zpsims 
     341            zpsils =  zpsims 
     342            !                                                          ! case of unstable atmospheric conditions 
     343            zobouku = MIN(  0.e0, -100.0 * 10.0 * ( zdtetar + 2.2e-3 * ztvmoyr )  ) 
     344            zxins   = ( 1. - 16. * zobouku )**0.25 
     345            zlxins  = LOG( ( 1. + zxins * zxins ) / 2. ) 
     346            zpsimu  = 2. * LOG( ( 1 + zxins ) * 0.5 )  + zlxins - 2. * ATAN( zxins ) + rpi * 0.5 
     347            zpsihu  = 2. * zlxins 
     348            zpsilu  = zpsihu 
     349            !                                                          ! intermediate values 
     350            zstab   = MAX( 0.e0, SIGN( 1.e0, zdteta ) ) 
     351            zpsim   = zstab * zpsimu + ( 1.0 - zstab ) * zpsims 
     352            zpsih   = zstab * zpsihu + ( 1.0 - zstab ) * zpsihs 
     353            zpsil   = zpsih 
     354             
     355            zvatmg         = MAX( 0.032 * 1.5e-3 * sf(jp_wndm)%fnow(ji,jj) * sf(jp_wndm)%fnow(ji,jj) / grav, zeps ) 
     356            zcmn           = vkarmn / LOG ( 10. / zvatmg ) 
     357            zchn           = 0.0327 * zcmn 
     358            zcln           = 0.0346 * zcmn 
     359            zcmcmn         = 1. / ( 1. - zcmn * zpsim / vkarmn ) 
     360            zchcm          = zcmcmn / ( 1. - zchn * zpsih / ( vkarmn * zcmn ) ) 
     361            zclcm          = zchcm 
     362            !                                                          ! transfert coef. (Large and Pond 1981,1982) 
     363            zcsho          = zchn * zchcm                                 
     364            zcleo          = zcln * zclcm  
     365 
     366            zrhova         = zrhoa * sf(jp_wndm)%fnow(ji,jj) 
     367 
     368            ! sensible heat flux 
     369            zqsb(ji,jj) = zrhova * zcsho * 1004.0  * ( zsst - ztatm )   
     370          
     371            ! latent heat flux (bounded by zero) 
     372            zqla(ji,jj) = MAX(  0.e0, zrhova * zcleo * 2.5e+06 * ( zqsato - zqatm )  ) 
     373            !                
     374         END DO 
     375      END DO 
     376       
     377      ! ----------------------------------------------------------------------------- ! 
     378      !     III    Total FLUXES                                                       ! 
     379      ! ----------------------------------------------------------------------------- ! 
     380!CDIR COLLAPSE 
     381      qns (:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)      ! Downward Non Solar flux 
     382!CDIR COLLAPSE 
     383      emp (:,:) = zqla(:,:) / cevap - sf(jp_prec)%fnow(:,:) / rday * tmask(:,:,1) 
     384!CDIR COLLAPSE 
     385      emps(:,:) = emp(:,:) 
     386      ! 
     387   END SUBROUTINE blk_oce_clio 
     388 
     389 
     390   SUBROUTINE blk_ice_clio(  pst   , pui   , pvi   , palb_cs, palb_os ,   & 
     391      &                      p_taui, p_tauj, p_qns , p_qsr,   & 
     392      &                      p_qla , p_dqns, p_dqla,          & 
     393      &                      p_tpr , p_spr ,                  & 
     394      &                      p_fr1 , p_fr2 , cd_grid  ) 
     395      !!--------------------------------------------------------------------------- 
     396      !!                     ***  ROUTINE blk_ice_clio  *** 
     397      !!                  
     398      !!  ** Purpose :   Computation of the heat fluxes at ocean and snow/ice 
     399      !!       surface the solar heat at ocean and snow/ice surfaces and the  
     400      !!       sensitivity of total heat fluxes to the SST variations 
     401      !!          
     402      !!  ** Method  :   The flux of heat at the ice and ocean surfaces are derived 
     403      !!       from semi-empirical ( or bulk ) formulae which relate the flux to  
     404      !!       the properties of the surface and of the lower atmosphere. Here, we 
     405      !!       follow the work of Oberhuber, 1988    
     406      !! 
     407      !!  ** Action  :   call flx_blk_albedo to compute ocean and ice albedo  
     408      !!          computation of snow precipitation 
     409      !!          computation of solar flux at the ocean and ice surfaces 
     410      !!          computation of the long-wave radiation for the ocean and sea/ice 
     411      !!          computation of turbulent heat fluxes over water and ice 
     412      !!          computation of evaporation over water 
     413      !!          computation of total heat fluxes sensitivity over ice (dQ/dT) 
     414      !!          computation of latent heat flux sensitivity over ice (dQla/dT) 
     415      !! 
     416      !!---------------------------------------------------------------------- 
     417      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpl) ::   pst      ! ice surface temperature                   [Kelvin] 
     418      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)     ::   pui      ! ice surface velocity (i-component, I-point)  [m/s] 
     419      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)     ::   pvi      ! ice surface velocity (j-component, I-point)  [m/s] 
     420      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpl) ::   palb_cs  ! ice albedo (clear sky) (alb_ice_cs)            [%] 
     421      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpl) ::   palb_os  ! ice albedo (overcast sky) (alb_ice_cs)         [%] 
     422      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)     ::   p_taui   ! surface ice stress at I-point (i-component) [N/m2] 
     423      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)     ::   p_tauj   ! surface ice stress at I-point (j-component) [N/m2] 
     424      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpl) ::   p_qns    ! non solar heat flux over ice (T-point)      [W/m2] 
     425      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpl) ::   p_qsr    !     solar heat flux over ice (T-point)      [W/m2] 
     426      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpl) ::   p_qla    ! latent    heat flux over ice (T-point)      [W/m2] 
     427      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpl) ::   p_dqns   ! non solar heat sensistivity  (T-point)      [W/m2] 
     428      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpl) ::   p_dqla   ! latent    heat sensistivity  (T-point)      [W/m2] 
     429      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)     ::   p_tpr    ! total precipitation          (T-point)       [Kg/m2/s] 
     430      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)     ::   p_spr    ! solid precipitation          (T-point)       [Kg/m2/s] 
     431      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)     ::   p_fr1    ! 1sr fraction of qsr penetration in ice             [%] 
     432      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)     ::   p_fr2    ! 2nd fraction of qsr penetration in ice             [%] 
     433      CHARACTER(len=1), INTENT(in   )                 ::   cd_grid  ! type of sea-ice grid ("C" or "B" grid) 
     434      !! 
     435      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
     436      !! 
     437      REAL(wp) ::   zcoef, zmt1, zmt2, zmt3, ztatm3             ! temporary scalars 
     438      REAL(wp) ::   ztaevbk, zind1, zind2, zind3, ztamr         !    -         - 
     439      REAL(wp) ::   zesi, zqsati, zdesidt                       !    -         - 
     440      REAL(wp) ::   zdqla, zcldeff, zev, zes, zpatm, zrhova     !    -         - 
     441      REAL(wp) ::   zcshi, zclei, zrhovaclei, zrhovacshi        !    -         - 
     442      REAL(wp) ::   ztice3, zticemb, zticemb2, zdqlw, zdqsb     !    -         - 
     443      !! 
     444      REAL(wp), DIMENSION(jpi,jpj) ::   ztatm   ! Tair in Kelvin 
     445      REAL(wp), DIMENSION(jpi,jpj) ::   zqatm   ! specific humidity 
     446      REAL(wp), DIMENSION(jpi,jpj) ::   zevsqr  ! vapour pressure square-root 
     447      REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa   ! air density 
     448      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qlw, z_qsb 
     449      !!--------------------------------------------------------------------- 
     450 
     451      zpatm = 101000.                        ! atmospheric pressure  (assumed constant  here) 
     452 
     453      !------------------------------------! 
     454      !   momentum fluxes  (utau, vtau )   ! 
     455      !------------------------------------! 
     456 
     457      SELECT CASE( cd_grid ) 
     458      CASE( 'C' )                          ! C-grid ice dynamics 
     459         ! Change from wind speed to wind stress over OCEAN (cao is used) 
     460         zcoef = cai / cao  
     461!CDIR COLLAPSE 
     462         DO jj = 1 , jpj 
     463            DO ji = 1, jpi 
     464               p_taui(ji,jj) = zcoef * utau(ji,jj) 
     465               p_tauj(ji,jj) = zcoef * vtau(ji,jj) 
     466            END DO 
     467         END DO 
     468      CASE( 'B' )                          ! B-grid ice dynamics 
     469         ! stress from ocean U- and V-points to ice U,V point 
     470         zcoef = cai / cao * 0.5 
     471!CDIR COLLAPSE 
     472         DO jj = 2, jpj 
     473            DO ji = fs_2, jpi   ! vector opt. 
     474               p_taui(ji,jj) = 0.5 * ( utau(ji-1,jj  ) + utau(ji-1,jj-1) ) 
     475               p_tauj(ji,jj) = 0.5 * ( vtau(ji  ,jj-1) + vtau(ji-1,jj-1) ) 
     476            END DO 
     477         END DO 
     478         CALL lbc_lnk( p_taui(:,:), 'I', -1. )   ! I-point (i.e. ice U-V point) 
     479         CALL lbc_lnk( p_tauj(:,:), 'I', -1. )   ! I-point (i.e. ice U-V point) 
     480      END SELECT 
     481 
     482 
     483      !  Determine cloud optical depths as a function of latitude (Chou et al., 1981). 
     484      !  and the correction factor for taking into account  the effect of clouds  
     485      !------------------------------------------------------ 
     486!CDIR NOVERRCHK 
     487!CDIR COLLAPSE 
     488      DO jj = 1, jpj 
     489!CDIR NOVERRCHK 
     490         DO ji = 1, jpi 
     491            ztatm (ji,jj) = sf(jp_tair)%fnow(ji,jj) + rt0            ! air temperature in Kelvins  
     492       
     493            zrhoa(ji,jj) = zpatm / ( 287.04 * ztatm(ji,jj) )         ! air density (equation of state for dry air)  
     494       
     495            ztamr = ztatm(ji,jj) - rtt                               ! Saturation water vapour 
     496            zmt1  = SIGN( 17.269,  ztamr ) 
     497            zmt2  = SIGN( 21.875,  ztamr ) 
     498            zmt3  = SIGN( 28.200, -ztamr ) 
     499            zes   = 611.0 * EXP(  ABS( ztamr ) * MIN ( zmt1, zmt2 )   & 
     500               &                / ( ztatm(ji,jj) - 35.86  + MAX( 0.e0, zmt3 ) )  ) 
     501 
     502            zev = sf(jp_humi)%fnow(ji,jj) * zes                      ! vapour pressure   
     503            zevsqr(ji,jj) = SQRT( zev * 0.01 )                       ! square-root of vapour pressure 
     504            zqatm(ji,jj) = 0.622 * zev / ( zpatm - 0.378 * zev )     ! specific humidity  
     505 
     506            !---------------------------------------------------- 
     507            !   Computation of snow precipitation (Ledley, 1985) | 
     508            !---------------------------------------------------- 
     509            zmt1  =   253.0 - ztatm(ji,jj)            ;   zind1 = MAX( 0.e0, SIGN( 1.e0, zmt1 ) ) 
     510            zmt2  = ( 272.0 - ztatm(ji,jj) ) / 38.0   ;   zind2 = MAX( 0.e0, SIGN( 1.e0, zmt2 ) ) 
     511            zmt3  = ( 281.0 - ztatm(ji,jj) ) / 18.0   ;   zind3 = MAX( 0.e0, SIGN( 1.e0, zmt3 ) ) 
     512            p_spr(ji,jj) = sf(jp_prec)%fnow(ji,jj) / rday   &        ! rday = converte mm/day to kg/m2/s 
     513               &         * (          zind1      &                   ! solid  (snow) precipitation [kg/m2/s] 
     514               &            + ( 1.0 - zind1 ) * (          zind2   * ( 0.5 + zmt2 )   & 
     515               &                                 + ( 1.0 - zind2 ) *  zind3 * zmt3  )   )  
     516 
     517            !----------------------------------------------------! 
     518            !  fraction of net penetrative shortwave radiation   ! 
     519            !----------------------------------------------------! 
     520            ! fraction of qsr_ice which is NOT absorbed in the thin surface layer 
     521            ! and thus which penetrates inside the ice cover ( Maykut and Untersteiner, 1971 ; Elbert anbd Curry, 1993 ) 
     522            p_fr1(ji,jj) = 0.18  * ( 1.e0 - sf(jp_ccov)%fnow(ji,jj) ) + 0.35 * sf(jp_ccov)%fnow(ji,jj)  
     523            p_fr2(ji,jj) = 0.82  * ( 1.e0 - sf(jp_ccov)%fnow(ji,jj) ) + 0.65 * sf(jp_ccov)%fnow(ji,jj) 
     524         END DO 
     525      END DO 
     526       
     527      !-----------------------------------------------------------! 
     528      !  snow/ice Shortwave radiation   (abedo already computed)  ! 
     529      !-----------------------------------------------------------! 
     530      CALL blk_clio_qsr_ice( palb_cs, palb_os, p_qsr ) 
     531 
     532      !                                     ! ========================== ! 
     533      DO jl = 1, jpl                        !  Loop over ice categories  ! 
     534         !                                  ! ========================== ! 
     535!CDIR NOVERRCHK 
     536!CDIR COLLAPSE 
     537         DO jj = 1 , jpj 
     538!CDIR NOVERRCHK 
     539            DO ji = 1, jpi 
     540               !-------------------------------------------! 
     541               !  long-wave radiation over ice categories  !  ( Berliand 1952 ; all latitudes ) 
     542               !-------------------------------------------! 
     543               ztatm3  = ztatm(ji,jj) * ztatm(ji,jj) * ztatm(ji,jj) 
     544               zcldeff = 1.0 - sbudyko(ji,jj) * sf(jp_ccov)%fnow(ji,jj) * sf(jp_ccov)%fnow(ji,jj)     
     545               ztaevbk = ztatm3 * ztatm(ji,jj) * zcldeff * ( 0.39 - 0.05 * zevsqr(ji,jj) )  
     546               ! 
     547               z_qlw(ji,jj,jl) = - emic * stefan * ( ztaevbk + 4. * ztatm3 * ( pst(ji,jj,jl) - ztatm(ji,jj) ) )  
     548 
     549               !---------------------------------------- 
     550               !  Turbulent heat fluxes over snow/ice     ( Latent and sensible )  
     551               !----------------------------------------         
     552 
     553               ! vapour pressure at saturation of ice (tmask to avoid overflow in the exponential) 
     554               zesi =  611.0 * EXP( 21.8745587 * tmask(ji,jj,1) * ( pst(ji,jj,jl) - rtt )/ ( pst(ji,jj,jl) - 7.66 ) ) 
     555               ! humidity close to the ice surface (at saturation) 
     556               zqsati   = ( 0.622 * zesi ) / ( zpatm - 0.378 * zesi ) 
     557                
     558               !  computation of intermediate values 
     559               zticemb  = pst(ji,jj,jl) - 7.66 
     560               zticemb2 = zticemb * zticemb   
     561               ztice3   = pst(ji,jj,jl) * pst(ji,jj,jl) * pst(ji,jj,jl) 
     562               zdesidt  = zesi * ( 9.5 * LOG( 10.0 ) * ( rtt - 7.66 )  / zticemb2 ) 
     563                
     564               !  Transfer cofficients assumed to be constant (Parkinson 1979 ; Maykut 1982) 
     565               zcshi    = 1.75e-03 
     566               zclei    = zcshi 
     567                
     568               !  sensible and latent fluxes over ice 
     569               zrhova     = zrhoa(ji,jj) * sf(jp_wndm)%fnow(ji,jj)      ! computation of intermediate values 
     570               zrhovaclei = zrhova * zcshi * 2.834e+06 
     571               zrhovacshi = zrhova * zclei * 1004.0 
     572             
     573               !  sensible heat flux 
     574               z_qsb(ji,jj,jl) = zrhovacshi * ( pst(ji,jj,jl) - ztatm(ji,jj) ) 
     575             
     576               !  latent heat flux  
     577               p_qla(ji,jj,jl) = MAX(  0.e0, zrhovaclei * ( zqsati - zqatm(ji,jj) )  ) 
     578               
     579               !  sensitivity of non solar fluxes (dQ/dT) (long-wave, sensible and latent fluxes) 
     580               zdqlw = 4.0 * emic * stefan * ztice3 
     581               zdqsb = zrhovacshi 
     582               zdqla = zrhovaclei * ( zdesidt * ( zqsati * zqsati / ( zesi * zesi ) ) * ( zpatm / 0.622 ) )    
     583               ! 
     584               p_dqla(ji,jj,jl) = zdqla                           ! latent flux sensitivity 
     585               p_dqns(ji,jj,jl) = -( zdqlw + zdqsb + zdqla )      !  total non solar sensitivity 
     586            END DO 
     587         END DO 
     588      END DO 
     589 
     590      ! 
     591      ! ----------------------------------------------------------------------------- ! 
     592      !     III    Total FLUXES                                                       ! 
     593      ! ----------------------------------------------------------------------------- ! 
     594      ! 
     595!CDIR COLLAPSE 
     596      p_qns(:,:,:) =     z_qlw (:,:,:) - z_qsb (:,:,:) - p_qla (:,:,:)      ! Downward Non Solar flux 
     597!CDIR COLLAPSE 
     598      p_tpr(:,:) = sf(jp_prec)%fnow(:,:) / rday     ! total precipitation [kg/m2/s] 
     599      ! 
     600!!gm : not necessary as all input data are lbc_lnk... 
     601      CALL lbc_lnk( p_fr1  (:,:) , 'T', 1. ) 
     602      CALL lbc_lnk( p_fr2  (:,:) , 'T', 1. ) 
     603      DO jl = 1, jpl 
     604         CALL lbc_lnk( p_qns (:,:,jl) , 'T', 1. ) 
     605         CALL lbc_lnk( p_dqns(:,:,jl) , 'T', 1. ) 
     606         CALL lbc_lnk( p_qla (:,:,jl) , 'T', 1. ) 
     607         CALL lbc_lnk( p_dqla(:,:,jl) , 'T', 1. ) 
     608      END DO 
     609 
     610!!gm : mask is not required on forcing 
     611      DO jl = 1, jpl 
     612         p_qns (:,:,jl) = p_qns (:,:,jl) * tmask(:,:,1) 
     613         p_qla (:,:,jl) = p_qla (:,:,jl) * tmask(:,:,1) 
     614         p_dqns(:,:,jl) = p_dqns(:,:,jl) * tmask(:,:,1) 
     615         p_dqla(:,:,jl) = p_dqla(:,:,jl) * tmask(:,:,1) 
     616      END DO 
     617   END SUBROUTINE blk_ice_clio 
     618 
     619 
     620   SUBROUTINE blk_clio_qsr_oce( pqsr_oce ) 
     621      !!--------------------------------------------------------------------------- 
     622      !!                     ***  ROUTINE blk_clio_qsr_oce  *** 
     623      !!                  
     624      !!  ** Purpose :   Computation of the shortwave radiation at the ocean and the 
     625      !!               snow/ice surfaces.  
     626      !!          
     627      !!  ** Method  : - computed qsr from the cloud cover for both ice and ocean  
     628      !!               - also initialise sbudyko and stauc once for all  
     629      !!---------------------------------------------------------------------- 
     630      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)     ::   pqsr_oce    ! shortwave radiation  over the ocean 
     631      !! 
     632      INTEGER, PARAMETER  ::   jp24 = 24   ! sampling of the daylight period (sunrise to sunset) into 24 equal parts 
     633      !!       
     634      INTEGER  ::   ji, jj, jt    ! dummy loop indices  
     635      INTEGER  ::   indaet            !  = -1, 0, 1 for odd, normal and leap years resp. 
     636      INTEGER  ::   iday              ! integer part of day 
     637      INTEGER  ::   indxb, indxc      ! index for cloud depth coefficient 
     638 
     639      REAL(wp)  ::   zalat , zclat, zcmue, zcmue2    ! local scalars  
     640      REAL(wp)  ::   zmt1, zmt2, zmt3                !  
     641      REAL(wp)  ::   zdecl, zsdecl , zcdecl          !  
     642      REAL(wp)  ::   za_oce, ztamr, zinda             ! 
     643 
     644      REAL(wp) ::   zdl, zlha     ! local scalars 
     645      REAL(wp) ::   zlmunoon, zcldcor, zdaycor            !    
     646      REAL(wp) ::   zxday, zdist, zcoef, zcoef1           ! 
     647      REAL(wp) ::   zes 
     648      !! 
     649      REAL(wp), DIMENSION(jpi,jpj) ::   zev         ! vapour pressure 
     650      REAL(wp), DIMENSION(jpi,jpj) ::   zdlha, zlsrise, zlsset     ! 2D workspace 
     651 
     652      REAL(wp), DIMENSION(jpi,jpj) ::   zps, zpc   ! sine (cosine) of latitude per sine (cosine) of solar declination  
     653      !!--------------------------------------------------------------------- 
     654 
     655 
     656      IF( lbulk_init ) THEN             !   Initilization at first time step only 
     657         rdtbs2 = nn_fsbc * rdt * 0.5 
     658         ! cloud optical depths as a function of latitude (Chou et al., 1981). 
     659         ! and the correction factor for taking into account  the effect of clouds  
     660         DO jj = 1, jpj 
     661            DO ji = 1 , jpi 
     662               zalat          = ( 90.e0 - ABS( gphit(ji,jj) ) ) /  5.e0 
     663               zclat          = ( 95.e0 -      gphit(ji,jj)   ) / 10.e0 
     664               indxb          = 1 + INT( zalat ) 
     665               indxc          = 1 + INT( zclat ) 
     666               zdl            = zclat - INT( zclat ) 
     667               !  correction factor to account for the effect of clouds 
     668               sbudyko(ji,jj) = budyko(indxb) 
     669               stauc  (ji,jj) = ( 1.e0 - zdl ) * tauco( indxc ) + zdl * tauco( indxc + 1 ) 
     670            END DO 
     671         END DO 
     672         IF    ( nleapy == 1 ) THEN   ;   yearday = 366.e0 
     673         ELSEIF( nleapy == 0 ) THEN   ;   yearday = 365.e0 
     674         ELSEIF( nleapy == 30) THEN   ;   yearday = 360.e0 
     675         ENDIF 
     676         lbulk_init = .FALSE. 
    216677      ENDIF 
    217678 
    218    END SUBROUTINE sbc_blk_clio 
    219  
    220  
    221    SUBROUTINE blk_oce_clio( ) 
    222    END SUBROUTINE blk_oce_clio 
    223  
    224  
    225    SUBROUTINE blk_ice_clio 
    226    END SUBROUTINE blk_ice_clio 
    227  
    228  
    229    SUBROUTINE blk_qsr_clio( ) 
    230    END SUBROUTINE blk_qsr_clio 
     679 
     680      ! Saturated water vapour and vapour pressure 
     681      ! ------------------------------------------ 
     682!CDIR NOVERRCHK 
     683!CDIR COLLAPSE 
     684      DO jj = 1, jpj 
     685!CDIR NOVERRCHK 
     686         DO ji = 1, jpi 
     687            ztamr = sf(jp_tair)%fnow(ji,jj) + rt0 - rtt 
     688            zmt1  = SIGN( 17.269,  ztamr ) 
     689            zmt2  = SIGN( 21.875,  ztamr ) 
     690            zmt3  = SIGN( 28.200, -ztamr ) 
     691            zes = 611.0 * EXP(  ABS( ztamr ) * MIN ( zmt1, zmt2 )   &              ! Saturation water vapour 
     692               &                     / ( sf(jp_tair)%fnow(ji,jj) + rt0 - 35.86  + MAX( 0.e0, zmt3 ) )  ) 
     693            zev(ji,jj) = sf(jp_humi)%fnow(ji,jj) * zes * 1.0e-05                   ! vapour pressure   
     694         END DO 
     695      END DO 
     696 
     697      !-----------------------------------! 
     698      !  Computation of solar irradiance  ! 
     699      !-----------------------------------! 
     700!!gm : hard coded  leap year ??? 
     701      indaet   = 1                                    ! = -1, 0, 1 for odd, normal and leap years resp. 
     702      zxday = nday_year + rdtbs2 / rday               ! day of the year at which the fluxes are calculated 
     703      iday  = INT( zxday )                            ! (centred at the middle of the ice time step) 
     704      CALL flx_blk_declin( indaet, iday, zdecl )      ! solar declination of the current day 
     705      zsdecl = SIN( zdecl * rad )                     ! its sine 
     706      zcdecl = COS( zdecl * rad )                     ! its cosine 
     707 
     708 
     709      !  correction factor added for computation of shortwave flux to take into account the variation of 
     710      !  the distance between the sun and the earth during the year (Oberhuber 1988) 
     711      zdist    = zxday * 2. * rpi / yearday 
     712      zdaycor  = 1.0 + 0.0013 * SIN( zdist ) + 0.0342 * COS( zdist ) 
     713 
     714!CDIR NOVERRCHK 
     715      DO jj = 1, jpj 
     716!CDIR NOVERRCHK 
     717         DO ji = 1, jpi 
     718            !  product of sine (cosine) of latitude and sine (cosine) of solar declination 
     719            zps(ji,jj) = SIN( gphit(ji,jj) * rad ) * zsdecl 
     720            zpc(ji,jj) = COS( gphit(ji,jj) * rad ) * zcdecl 
     721            !  computation of the both local time of sunrise and sunset 
     722            zlsrise(ji,jj) = ACOS( - SIGN( 1.e0, zps(ji,jj) )    & 
     723               &                   * MIN(  1.e0, SIGN( 1.e0, zps(ji,jj) ) * ( zps(ji,jj) / zpc(ji,jj) )  )   ) 
     724            zlsset (ji,jj) = - zlsrise(ji,jj) 
     725            !  dividing the solar day into jp24 segments of length zdlha 
     726            zdlha  (ji,jj) = ( zlsrise(ji,jj) - zlsset(ji,jj) ) / REAL( jp24 ) 
     727         END DO 
     728      END DO 
     729 
     730 
     731      !---------------------------------------------! 
     732      !  shortwave radiation absorbed by the ocean  ! 
     733      !---------------------------------------------! 
     734      pqsr_oce(:,:)   = 0.e0      ! set ocean qsr to zero       
     735 
     736      ! compute and sum ocean qsr over the daylight (i.e. between sunrise and sunset) 
     737!CDIR NOVERRCHK    
     738      DO jt = 1, jp24 
     739         zcoef = FLOAT( jt ) - 0.5 
     740!CDIR NOVERRCHK      
     741!CDIR COLLAPSE 
     742         DO jj = 1, jpj 
     743!CDIR NOVERRCHK 
     744            DO ji = 1, jpi 
     745               zlha = COS(  zlsrise(ji,jj) - zcoef * zdlha(ji,jj)  )                  ! local hour angle 
     746               zcmue              = MAX( 0.e0 ,   zps(ji,jj) + zpc(ji,jj) * zlha  )   ! cos of local solar altitude 
     747               zcmue2             = 1368.0 * zcmue * zcmue 
     748 
     749               ! ocean albedo depending on the cloud cover (Payne, 1972) 
     750               za_oce     = ( 1.0 - sf(jp_ccov)%fnow(ji,jj) ) * 0.05 / ( 1.1 * zcmue**1.4 + 0.15 )   &   ! clear sky 
     751                  &       +         sf(jp_ccov)%fnow(ji,jj)   * 0.06                                     ! overcast 
     752 
     753                  ! solar heat flux absorbed by the ocean (Zillman, 1972) 
     754               pqsr_oce(ji,jj) = pqsr_oce(ji,jj)                                         & 
     755                  &            + ( 1.0 - za_oce ) * zdlha(ji,jj) * zcmue2                & 
     756                  &            / ( ( zcmue + 2.7 ) * zev(ji,jj) + 1.085 * zcmue +  0.10 ) 
     757            END DO 
     758         END DO 
     759      END DO 
     760      ! Taking into account the ellipsity of the earth orbit, the clouds AND masked if sea-ice cover > 0% 
     761!!gm : bug zinda is always 0 si ice.... 
     762      zcoef1 = srgamma * zdaycor / ( 2. * rpi ) 
     763!CDIR COLLAPSE 
     764      DO jj = 1, jpj 
     765         DO ji = 1, jpi 
     766            zlmunoon = ASIN( zps(ji,jj) + zpc(ji,jj) ) / rad                         ! local noon solar altitude 
     767            zcldcor  = MIN(  1.e0, ( 1.e0 - 0.62 * sf(jp_ccov)%fnow(ji,jj)   &       ! cloud correction (Reed 1977) 
     768               &                          + 0.0019 * zlmunoon )                 ) 
     769            zinda    = MAX(  0.e0, SIGN(  1.e0, -( -1.5 + freeze(ji,jj) )  )   )            ! 0 if more than 0% of ice 
     770            pqsr_oce(ji,jj) = zcoef1 * zcldcor * zinda * pqsr_oce(ji,jj) * tmask(ji,jj,1)   ! and zcoef1: ellipsity 
     771         END DO 
     772      END DO 
     773 
     774   END SUBROUTINE blk_clio_qsr_oce 
     775 
     776 
     777   SUBROUTINE blk_clio_qsr_ice( pa_ice_cs, pa_ice_os, pqsr_ice ) 
     778      !!--------------------------------------------------------------------------- 
     779      !!                     ***  ROUTINE blk_clio_qsr_ice  *** 
     780      !!                  
     781      !!  ** Purpose :   Computation of the shortwave radiation at the ocean and the 
     782      !!               snow/ice surfaces.  
     783      !!          
     784      !!  ** Method  : - computed qsr from the cloud cover for both ice and ocean  
     785      !!               - also initialise sbudyko and stauc once for all  
     786      !!---------------------------------------------------------------------- 
     787      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpl) ::   pa_ice_cs   ! albedo of ice under clear sky 
     788      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpl) ::   pa_ice_os   ! albedo of ice under overcast sky 
     789      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpl) ::   pqsr_ice    ! shortwave radiation over the ice/snow 
     790      !! 
     791      INTEGER, PARAMETER  ::   jp24 = 24   ! sampling of the daylight period (sunrise to sunset) into 24 equal parts 
     792      !! 
     793      INTEGER  ::   ji, jj, jl, jt    ! dummy loop indices 
     794      INTEGER  ::   indaet            !  = -1, 0, 1 for odd, normal and leap years resp. 
     795      INTEGER  ::   iday              ! integer part of day 
     796 
     797      REAL(wp)  ::   zcmue, zcmue2    ! local scalars  
     798      REAL(wp)  ::   zmt1, zmt2, zmt3                !  
     799      REAL(wp)  ::   zdecl, zsdecl , zcdecl          !  
     800      REAL(wp)  ::   ztamr             ! 
     801 
     802      REAL(wp) ::   zlha     ! local scalars 
     803      REAL(wp) ::   zdaycor            ! 
     804      REAL(wp) ::   zxday, zdist, zcoef, zcoef1           ! 
     805      REAL(wp) ::   zes 
     806      REAL(wp) ::   zqsr_ice_cs, zqsr_ice_os 
     807      !! 
     808      REAL(wp), DIMENSION(jpi,jpj) ::   zev         ! vapour pressure 
     809      REAL(wp), DIMENSION(jpi,jpj) ::   zdlha, zlsrise, zlsset     ! 2D workspace 
     810        
     811      REAL(wp), DIMENSION(jpi,jpj) ::   zps, zpc   ! sine (cosine) of latitude per sine (cosine) of solar declination  
     812      !!--------------------------------------------------------------------- 
     813       
     814      ! Saturated water vapour and vapour pressure 
     815      ! ------------------------------------------ 
     816!CDIR NOVERRCHK 
     817!CDIR COLLAPSE 
     818      DO jj = 1, jpj 
     819!CDIR NOVERRCHK 
     820         DO ji = 1, jpi            
     821            ztamr = sf(jp_tair)%fnow(ji,jj) + rt0 - rtt            
     822            zmt1  = SIGN( 17.269,  ztamr ) 
     823            zmt2  = SIGN( 21.875,  ztamr ) 
     824            zmt3  = SIGN( 28.200, -ztamr ) 
     825            zes = 611.0 * EXP(  ABS( ztamr ) * MIN ( zmt1, zmt2 )   &              ! Saturation water vapour 
     826               &                     / ( sf(jp_tair)%fnow(ji,jj) + rt0 - 35.86  + MAX( 0.e0, zmt3 ) )  ) 
     827            zev(ji,jj) = sf(jp_humi)%fnow(ji,jj) * zes * 1.0e-05                   ! vapour pressure   
     828         END DO 
     829      END DO 
     830 
     831      !-----------------------------------! 
     832      !  Computation of solar irradiance  ! 
     833      !-----------------------------------! 
     834!!gm : hard coded  leap year ??? 
     835      indaet   = 1                                    ! = -1, 0, 1 for odd, normal and leap years resp. 
     836      zxday = nday_year + rdtbs2 / rday               ! day of the year at which the fluxes are calculated 
     837      iday  = INT( zxday )                            ! (centred at the middle of the ice time step) 
     838      CALL flx_blk_declin( indaet, iday, zdecl )      ! solar declination of the current day 
     839      zsdecl = SIN( zdecl * rad )                     ! its sine 
     840      zcdecl = COS( zdecl * rad )                     ! its cosine 
     841 
     842       
     843      !  correction factor added for computation of shortwave flux to take into account the variation of 
     844      !  the distance between the sun and the earth during the year (Oberhuber 1988) 
     845      zdist    = zxday * 2. * rpi / yearday 
     846      zdaycor  = 1.0 + 0.0013 * SIN( zdist ) + 0.0342 * COS( zdist ) 
     847 
     848!CDIR NOVERRCHK 
     849      DO jj = 1, jpj 
     850!CDIR NOVERRCHK 
     851         DO ji = 1, jpi 
     852            !  product of sine (cosine) of latitude and sine (cosine) of solar declination 
     853            zps(ji,jj) = SIN( gphit(ji,jj) * rad ) * zsdecl 
     854            zpc(ji,jj) = COS( gphit(ji,jj) * rad ) * zcdecl 
     855            !  computation of the both local time of sunrise and sunset 
     856            zlsrise(ji,jj) = ACOS( - SIGN( 1.e0, zps(ji,jj) )    & 
     857               &                   * MIN(  1.e0, SIGN( 1.e0, zps(ji,jj) ) * ( zps(ji,jj) / zpc(ji,jj) )  )   )  
     858            zlsset (ji,jj) = - zlsrise(ji,jj) 
     859            !  dividing the solar day into jp24 segments of length zdlha 
     860            zdlha  (ji,jj) = ( zlsrise(ji,jj) - zlsset(ji,jj) ) / REAL( jp24 ) 
     861         END DO 
     862      END DO 
     863 
     864 
     865      !---------------------------------------------! 
     866      !  shortwave radiation absorbed by the ice    ! 
     867      !---------------------------------------------! 
     868      ! compute and sum ice qsr over the daylight for each ice categories 
     869      pqsr_ice(:,:,:) = 0.e0 
     870      zcoef1 = zdaycor / ( 2. * rpi ) 
     871       
     872      !                    !----------------------------!  
     873      DO jl = 1, jpl       !  loop over ice categories  ! 
     874         !                 !----------------------------!  
     875!CDIR NOVERRCHK    
     876         DO jt = 1, jp24    
     877            zcoef = FLOAT( jt ) - 0.5 
     878!CDIR NOVERRCHK      
     879!CDIR COLLAPSE 
     880            DO jj = 1, jpj 
     881!CDIR NOVERRCHK 
     882               DO ji = 1, jpi 
     883                  zlha = COS(  zlsrise(ji,jj) - zcoef * zdlha(ji,jj)  )                  ! local hour angle 
     884                  zcmue              = MAX( 0.e0 ,   zps(ji,jj) + zpc(ji,jj) * zlha  )   ! cos of local solar altitude 
     885                  zcmue2             = 1368.0 * zcmue * zcmue 
     886                   
     887                  !  solar heat flux absorbed by the ice/snow system (Shine and Crane 1984 adapted to high albedo)  
     888                  zqsr_ice_cs =  ( 1.0 - pa_ice_cs(ji,jj,jl) ) * zdlha(ji,jj) * zcmue2        &   ! clear sky 
     889                     &        / ( ( 1.0 + zcmue ) * zev(ji,jj) + 1.2 * zcmue + 0.0455 ) 
     890                  zqsr_ice_os = zdlha(ji,jj) * SQRT( zcmue )                                  &   ! overcast sky 
     891                     &        * ( 53.5 + 1274.5 * zcmue )      * ( 1.0 - 0.996  * pa_ice_os(ji,jj,jl) )    & 
     892                     &        / (  1.0 + 0.139  * stauc(ji,jj) * ( 1.0 - 0.9435 * pa_ice_os(ji,jj,jl) ) )        
     893              
     894                  pqsr_ice(ji,jj,jl) = pqsr_ice(ji,jj,jl) + (  ( 1.0 - sf(jp_ccov)%fnow(ji,jj) ) * zqsr_ice_cs    & 
     895                     &                                       +         sf(jp_ccov)%fnow(ji,jj)   * zqsr_ice_os  ) 
     896               END DO 
     897            END DO 
     898         END DO 
     899         ! 
     900         ! Correction : Taking into account the ellipsity of the earth orbit 
     901         pqsr_ice(:,:,jl) = pqsr_ice(:,:,jl) * zcoef1 * tmask(:,:,1) 
     902         ! 
     903         !                 !--------------------------------!  
     904      END DO               !  end loop over ice categories  ! 
     905         !                 !--------------------------------!  
     906 
     907 
     908!!gm  : this should be suppress as input data have been passed through lbc_lnk 
     909      DO jl = 1, jpl 
     910         CALL lbc_lnk( pqsr_ice(:,:,jl) , 'T', 1. ) 
     911      END DO 
     912      ! 
     913   END SUBROUTINE blk_clio_qsr_ice 
    231914 
    232915 
     
    236919      !!           
    237920      !! ** Purpose :   Computation of the solar declination for the day 
    238       !!         kday ( in decimal degrees ). 
    239921      !!        
    240       !! ** Method  : 
     922      !! ** Method  :   ??? 
    241923      !!--------------------------------------------------------------------- 
    242924      INTEGER , INTENT(in   ) ::   ky      ! = -1, 0, 1 for odd, normal and leap years resp. 
    243925      INTEGER , INTENT(in   ) ::   kday    ! day of the year ( kday = 1 on january 1) 
    244926      REAL(wp), INTENT(  out) ::   pdecl   ! solar declination 
    245  
    246       REAL(wp) ::   zday              ,  &  ! corresponding day of type year (cf. ky) 
    247          &          zp1, zp2, zp3, zp4      ! temporary scalars 
    248       REAL(wp) ::   a0  =  0.39507671 ,  &  ! constants used in  solar  
    249          &          a1  = 22.85684301 ,  &  ! declinaison computation 
    250          &          a2  = -0.38637317 ,  & 
    251          &          a3  =  0.15096535 ,  & 
    252          &          a4  = -0.00961411 ,  & 
    253          &          b1  = -4.29692073 ,  & 
    254          &          b2  =  0.05702074 ,  & 
    255          &          b3  = -0.09028607 ,  & 
    256          &          b4  =  0.00592797 
     927      !! 
     928      REAL(wp) ::   a0  =  0.39507671      ! coefficients for solar declinaison computation 
     929      REAL(wp) ::   a1  = 22.85684301      !     "              ""                 " 
     930      REAL(wp) ::   a2  = -0.38637317      !     "              ""                 " 
     931      REAL(wp) ::   a3  =  0.15096535      !     "              ""                 " 
     932      REAL(wp) ::   a4  = -0.00961411      !     "              ""                 " 
     933      REAL(wp) ::   b1  = -4.29692073      !     "              ""                 " 
     934      REAL(wp) ::   b2  =  0.05702074      !     "              ""                 " 
     935      REAL(wp) ::   b3  = -0.09028607      !     "              ""                 " 
     936      REAL(wp) ::   b4  =  0.00592797 
     937      !! 
     938      REAL(wp) ::   zday   ! corresponding day of type year (cf. ky) 
     939      REAL(wp) ::   zp     ! temporary scalars 
    257940      !!--------------------------------------------------------------------- 
    258       ! 
    259       SELECT CASE ( ky ) 
    260       CASE ( 1 ) 
    261          zday = REAL( kday, wp ) - 0.5 
    262       CASE ( 3 ) 
    263          zday = REAL( kday, wp ) - 1.0 
    264       CASE DEFAULT 
    265          zday = REAL( kday, wp )  
    266       END SELECT 
    267941             
    268       zp1 = rpi * ( 2.0 * zday - 367.0 ) / yearday 
    269       zp2 = 2. * zp1 
    270       zp3 = 3. * zp1 
    271       zp4 = 4. * zp1 
     942      IF    ( ky == 1 )  THEN   ;   zday = REAL( kday ) - 0.5 
     943      ELSEIF( ky == 3 )  THEN   ;   zday = REAL( kday ) - 1. 
     944      ELSE                      ;   zday = REAL( kday ) 
     945      ENDIF 
     946       
     947      zp = rpi * ( 2.0 * zday - 367.0 ) / yearday 
    272948       
    273949      pdecl  = a0                                                                      & 
    274          &   + a1 * COS( zp1 ) + a2 * COS( zp2 ) + a3 * COS( zp3 ) + a4 * COS( zp4 )   & 
    275          &   + b1 * SIN( zp1 ) + b2 * SIN( zp2 ) + b3 * SIN( zp3 ) + b4 * SIN( zp4 ) 
     950         &   + a1 * COS( zp ) + a2 * COS( 2. * zp ) + a3 * COS( 3. * zp ) + a4 * COS( 4. * zp )   & 
     951         &   + b1 * SIN( zp ) + b2 * SIN( 2. * zp ) + b3 * SIN( 3. * zp ) + b4 * SIN( 4. * zp ) 
    276952      ! 
    277953   END SUBROUTINE flx_blk_declin 
  • branches/dev_001_SBC/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r751 r840  
    125125         SELECT CASE( kblk ) 
    126126         CASE( 3 )           ! CLIO bulk formulation 
    127             CALL blk_ice_clio() 
     127            CALL blk_ice_clio( sist     , ui_ice   , vi_ice   , alb_ice_cs, alb_ice_os,            & 
     128               &                                     utaui_ice, vtaui_ice , qns_ice   , qsr_ice,   & 
     129               &                                     qla_ice  , dqns_ice  , dqla_ice  ,            & 
     130               &                                     tprecip  , sprecip   ,                        & 
     131               &                                     fr1_i0   , fr2_i0    , 'B'  ) 
    128132         CASE( 4 )           ! CORE bulk formulation 
    129133            CALL blk_ice_core( sist     , ui_ice   , vi_ice   , alb_ice_cs,                      & 
Note: See TracChangeset for help on using the changeset viewer.