Changeset 1858


Ignore:
Timestamp:
2010-05-04T10:39:48+02:00 (10 years ago)
Author:
gm
Message:

ticket:#665 : step 1 - heat content of freezing-melting ice

Location:
branches/DEV_r1837_mass_heat_salt_fluxes/NEMO
Files:
14 edited

Legend:

Unmodified
Added
Removed
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/ice_2.F90

    r1857 r1858  
    44   !! Sea Ice physics:  diagnostics variables of ice defined in memory 
    55   !!===================================================================== 
    6    !! History :  2.0  !  03-08  (C. Ethe)  F90: Free form and module 
     6   !! History :  LIM  ! 2003-08  (C. Ethe)  F90: Free form and module 
     7   !!            2.0  ! 2010-05  (Y. Aksenov, M. Vancoppenolle, G. Madec) add heat content exchanges 
    78   !!---------------------------------------------------------------------- 
    89#if defined key_lim2 
     
    6869   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qstoif        !: Energy stored in the brine pockets 
    6970   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fbif          !: Heat flux at the ice base 
    70    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdmsnif       !: Variation of snow mass 
    71    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdmicif       !: Variation of ice mass 
     71   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdm_snw       !: Variation of snow mass over 1 time step           [Kg/m2] 
     72   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdq_snw       !: heat content associated to rdm_snw                [J/m2] 
     73   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdq_ice       !: Variation of ice  mass over 1 time step           [Kg/m2] 
     74   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdq_ice       !: heat content associated to rdm_snw                [J/m2] 
    7275   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qldif         !: heat balance of the lead (or of the open ocean) 
    7376   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qcmif         !: Energy needed to bring the ocean surface layer until its freezing  
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limsbc_2.F90

    r1857 r1858  
    44   !!           computation of the flux at the sea ice/ocean interface 
    55   !!====================================================================== 
    6    !! History : 00-01 (H. Goosse) Original code 
    7    !!           02-07 (C. Ethe, G. Madec) re-writing F90 
    8    !!           06-07 (G. Madec) surface module 
     6   !! History : LIM  ! 2000-01 (H. Goosse) Original code 
     7   !!           2.0  ! 2002-07 (C. Ethe, G. Madec) re-writing F90 
     8   !!            -   ! 2006-07 (G. Madec) surface module 
     9   !!           2.1  ! 2010-05  (Y. Aksenov, M. Vancoppenolle, G. Madec) add heat content exchanges 
    910   !!---------------------------------------------------------------------- 
    1011#if defined key_lim2 
    1112   !!---------------------------------------------------------------------- 
    1213   !!   'key_lim2'                                    LIM 2.0 sea-ice model 
    13    !!---------------------------------------------------------------------- 
    1414   !!---------------------------------------------------------------------- 
    1515   !!   lim_sbc_2  : flux at the ice / ocean interface 
     
    211211#if defined key_coupled 
    212212          zemp = emp_tot(ji,jj) - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  
    213              &   + rdmsnif(ji,jj) * zrdtir                                      !  freshwaterflux due to snow melting  
     213             &   + rdm_snw(ji,jj) * zrdtir                                      !  freshwaterflux due to snow melting  
    214214#else 
    215215!!$            !  computing freshwater exchanges at the ice/ocean interface 
     
    217217!!$               &   + tprecip(ji,jj)                           &   !  total precipitation 
    218218!!$               &   - sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )   &   !  remov. snow precip over ice 
    219 !!$               &   - rdmsnif(ji,jj) / rdt_ice                     !  freshwaterflux due to snow melting  
     219!!$               &   - rdm_snw(ji,jj) / rdt_ice                     !  freshwaterflux due to snow melting  
    220220            !  computing freshwater exchanges at the ice/ocean interface 
    221221            zemp = + emp(ji,jj)     *         frld(ji,jj)      &   !  e-p budget over open ocean fraction  
    222222               &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &   !  liquid precipitation reaches directly the ocean 
    223223               &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  taking into account change in ice cover within the time step 
    224                &   + rdmsnif(ji,jj) * zrdtir                       !  freshwaterflux due to snow melting  
     224               &   + rdm_snw(ji,jj) * zrdtir                       !  freshwaterflux due to snow melting  
    225225               !                                                   !  ice-covered fraction: 
    226226#endif             
    227227 
    228228            !  computing salt exchanges at the ice/ocean interface 
    229             zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdmicif(ji,jj) * zrdtir )  
     229            zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdm_ice(ji,jj) * zrdtir )  
    230230             
    231231            !  converting the salt flux from ice to a freshwater flux from ocean 
     
    239239 
    240240      IF( lk_diaar5 ) THEN 
    241          CALL iom_put( 'isnwmlt_cea'  ,                 rdmsnif(:,:) * zrdtir ) 
    242          CALL iom_put( 'fsal_virt_cea',   soce_r(:,:) * rdmicif(:,:) * zrdtir ) 
    243          CALL iom_put( 'fsal_real_cea', - sice_r(:,:) * rdmicif(:,:) * zrdtir ) 
     241         CALL iom_put( 'isnwmlt_cea'  ,                 rdm_snw(:,:) * zrdtir ) 
     242         CALL iom_put( 'fsal_virt_cea',   soce_r(:,:) * rdm_ice(:,:) * zrdtir ) 
     243         CALL iom_put( 'fsal_real_cea', - sice_r(:,:) * rdm_ice(:,:) * zrdtir ) 
    244244      ENDIF 
    245245 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_2.F90

    r1857 r1858  
    44   !!              LIM thermo ice model : ice thermodynamic 
    55   !!====================================================================== 
    6    !! History :  1.0  ! 2000-01 (LIM) 
    7    !!            2.0  ! 2002-07 (C. Ethe, G. Madec) F90 
    8    !!            2.0  ! 2003-08 (C. Ethe)  add lim_thd_init 
    9    !!             -   ! 2008-2008  (A. Caubel, G. Madec, E. Maisonnave, S. Masson ) generic coupled interface 
     6   !! History :  1.0  !  2000-01  (LIM) Original code 
     7   !!            2.0  !  2002-07  (C. Ethe, G. Madec) F90 
     8   !!            2.0  !  2003-08  (C. Ethe)  add lim_thd_init 
     9   !!             -   !  2008-08  (A. Caubel, G. Madec, E. Maisonnave, S. Masson ) generic coupled interface 
     10   !!            2.1  !  2010-05  (Y. Aksenov G. Madec) add heat associated with mass exchanges 
    1011   !!--------------------------------------------------------------------- 
    1112#if defined key_lim2 
     
    8889      REAL(wp) ::   za , zh, zthsnice    ! 
    8990      REAL(wp) ::   zfric_u              ! friction velocity  
    90       REAL(wp) ::   zfnsol               ! total non solar heat 
    91       REAL(wp) ::   zfontn               ! heat flux from snow thickness 
    9291      REAL(wp) ::   zfntlat, zpareff     ! test. the val. of lead heat budget 
    9392      REAL(wp), DIMENSION(jpi,jpj)     ::   ztmp      ! 2D workspace 
     
    132131      ffltbif(:,:) = 0.e0   ! linked with fstric 
    133132      qfvbq  (:,:) = 0.e0   ! linked with fstric 
    134       rdmsnif(:,:) = 0.e0   ! variation of snow mass per unit area 
    135       rdmicif(:,:) = 0.e0   ! variation of ice mass per unit area 
     133      rdm_snw(:,:) = 0.e0   ! variation of snow mass over 1 time step 
     134      rdq_snw(:,:) = 0.e0   ! heat content associated with rdm_snw 
     135      rdm_ice(:,:) = 0.e0   ! variation of ice mass over 1 time step 
     136      rdq_ice(:,:) = 0.e0   ! heat content associated with rdm_ice 
    136137      zmsk (:,:,:) = 0.e0 
    137138 
     
    193194      !      Partial computation of forcing for the thermodynamic sea ice model. 
    194195      !-------------------------------------------------------------------------- 
    195  
    196       sst_m(:,:) = sst_m(:,:) + rt0 
    197  
    198 !CDIR NOVERRCHK 
    199       DO jj = 1, jpj 
    200 !CDIR NOVERRCHK 
     196      !CDIR NOVERRCHK 
     197      DO jj = 1, jpj 
     198         !CDIR NOVERRCHK 
    201199         DO ji = 1, jpi 
    202200            zthsnice       = hsnif(ji,jj) + hicif(ji,jj) 
     
    211209            !  net downward heat flux from the ice to the ocean, expressed as a function of ocean  
    212210            !  temperature and turbulent mixing (McPhee, 1992) 
    213             zfric_u        = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin )  ! friction velocity 
    214             fdtcn(ji,jj)  = zindb * rau0 * rcp * 0.006  * zfric_u * ( sst_m(ji,jj) - tfu(ji,jj) )  
     211            zfric_u       = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin )  ! friction velocity 
     212            fdtcn(ji,jj)  = zindb * rau0 * rcp * 0.006  * zfric_u * ( sst_m(ji,jj) + rt0 - tfu(ji,jj) )  
    215213            qdtcn(ji,jj)  = zindb * fdtcn(ji,jj) * frld(ji,jj) * rdt_ice 
    216214                         
    217215            !  partial computation of the lead energy budget (qldif) 
    218216#if defined key_coupled  
    219             qldif(ji,jj)   = tms(ji,jj) * rdt_ice                                             & 
     217            qldif(ji,jj)   = tms(ji,jj) * rdt_ice                                                  & 
    220218               &    * (   ( qsr_tot(ji,jj) - qsr_ice(ji,jj,1) * zfricp ) * ( 1.0 - thcm(ji,jj) )   & 
    221219               &        + ( qns_tot(ji,jj) - qns_ice(ji,jj,1) * zfricp )                           & 
    222220               &        + frld(ji,jj) * ( fdtcn(ji,jj) + ( 1.0 - zindb ) * fsbbq(ji,jj) )   ) 
    223221#else 
    224             zfontn         = ( sprecip(ji,jj) / rhosn ) * xlsn  !   energy for melting solid precipitation 
    225             zfnsol         = qns(ji,jj)                         !  total non solar flux over the ocean 
    226             qldif(ji,jj)   = tms(ji,jj) * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) )   & 
    227                &                               + zfnsol + fdtcn(ji,jj) - zfontn     & 
    228                &                               + ( 1.0 - zindb ) * fsbbq(ji,jj) )   & 
    229                &                        * frld(ji,jj) * rdt_ice     
    230 !!$            qldif(ji,jj)   = tms(ji,jj) * rdt_ice * frld(ji,jj)  
    231 !!$               &           * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) )      & 
    232 !!$               &             + qns(ji,jj)  + fdtcn(ji,jj) - zfontn     & 
    233 !!$               &             + ( 1.0 - zindb ) * fsbbq(ji,jj)      )   & 
     222            qldif(ji,jj)   = tms(ji,jj) * rdt_ice * frld(ji,jj)                    & 
     223               &                        * (  qsr(ji,jj) * ( 1.0 - thcm(ji,jj) )    & 
     224               &                           + qns(ji,jj)  +  fdtcn(ji,jj)           & 
     225               &                           + ( 1.0 - zindb ) * fsbbq(ji,jj)      ) 
    234226#endif 
    235227            !  parlat : percentage of energy used for lateral ablation (0.0)  
     
    241233             
    242234            !  energy needed to bring ocean surface layer until its freezing 
    243             qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj,1)   & 
    244                 &          * ( tfu(ji,jj) - sst_m(ji,jj) ) * ( 1 - zinda ) 
     235            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj,1) * ( tfu(ji,jj) - sst_m(ji,jj) - rt0 ) * ( 1 - zinda ) 
    245236             
    246237            !  calculate oceanic heat flux. 
     
    251242         END DO 
    252243      END DO 
    253        
    254       sst_m(:,:) = sst_m(:,:) - rt0 
    255                 
     244                      
    256245      !         Select icy points and fulfill arrays for the vectorial grid. 
    257246      !---------------------------------------------------------------------- 
     
    307296         CALL tab_2d_1d_2( nbpb, qldif_1d   (1:nbpb)     , qldif      , jpi, jpj, npb(1:nbpb) ) 
    308297         CALL tab_2d_1d_2( nbpb, qstbif_1d  (1:nbpb)     , qstoif     , jpi, jpj, npb(1:nbpb) ) 
    309          CALL tab_2d_1d_2( nbpb, rdmicif_1d (1:nbpb)     , rdmicif    , jpi, jpj, npb(1:nbpb) ) 
     298!!gm bug ??  rdm_snw, rdq_snw should be transformed into 1D fields ?  like rdm_ice and rdq_ice ??? 
     299!!gm         I think that  in fact the transformation is useless for both snow and ice as they are already set to zero 
     300!!gm         at the begining of the routine and only the value over icy point is updated in 1D and converted into 2D 
     301!!gm         the two lines below are thus useless. 
     302         CALL tab_2d_1d_2( nbpb, rdm_ice_1d (1:nbpb)     , rdm_ice    , jpi, jpj, npb(1:nbpb) ) 
     303         CALL tab_2d_1d_2( nbpb, rdq_ice_1d (1:nbpb)     , rdq_ice    , jpi, jpj, npb(1:nbpb) ) 
     304!!gm bug end. 
    310305         CALL tab_2d_1d_2( nbpb, dmgwi_1d   (1:nbpb)     , dmgwi      , jpi, jpj, npb(1:nbpb) ) 
    311306         CALL tab_2d_1d_2( nbpb, qlbbq_1d   (1:nbpb)     , zqlbsbq    , jpi, jpj, npb(1:nbpb) ) 
     
    327322         CALL tab_1d_2d_2( nbpb, qfvbq      , npb, qfvbq_1d  (1:nbpb)     , jpi, jpj ) 
    328323         CALL tab_1d_2d_2( nbpb, qstoif     , npb, qstbif_1d (1:nbpb)     , jpi, jpj ) 
    329          CALL tab_1d_2d_2( nbpb, rdmicif    , npb, rdmicif_1d(1:nbpb)     , jpi, jpj ) 
     324         CALL tab_1d_2d_2( nbpb, rdm_ice    , npb, rdm_ice_1d(1:nbpb)     , jpi, jpj ) 
     325         CALL tab_1d_2d_2( nbpb, rdq_ice    , npb, rdq_ice_1d(1:nbpb)     , jpi, jpj ) 
     326         CALL tab_1d_2d_2( nbpb, rdm_snw    , npb, rdm_snw_1d(1:nbpb)     , jpi, jpj ) 
     327         CALL tab_1d_2d_2( nbpb, rdq_snw    , npb, rdq_snw_1d(1:nbpb)     , jpi, jpj ) 
    330328         CALL tab_1d_2d_2( nbpb, dmgwi      , npb, dmgwi_1d  (1:nbpb)     , jpi, jpj ) 
    331          CALL tab_1d_2d_2( nbpb, rdmsnif    , npb, rdmsnif_1d(1:nbpb)     , jpi, jpj ) 
    332329         CALL tab_1d_2d_2( nbpb, zdvosif    , npb, dvsbq_1d  (1:nbpb)     , jpi, jpj ) 
    333330         CALL tab_1d_2d_2( nbpb, zdvobif    , npb, dvbbq_1d  (1:nbpb)     , jpi, jpj ) 
     
    387384      IF( nbpac > 0 ) THEN 
    388385         ! 
    389          zlicegr(:,:) = rdmicif(:,:)      ! to output the lateral sea-ice growth  
     386         zlicegr(:,:) = rdm_ice(:,:)      ! to output the lateral sea-ice growth  
    390387         !...Put the variable in a 1-D array for lateral accretion 
    391388         CALL tab_2d_1d_2( nbpac, frld_1d   (1:nbpac)     , frld       , jpi, jpj, npac(1:nbpac) ) 
     
    398395         CALL tab_2d_1d_2( nbpac, qcmif_1d  (1:nbpac)     , qcmif      , jpi, jpj, npac(1:nbpac) ) 
    399396         CALL tab_2d_1d_2( nbpac, qstbif_1d (1:nbpac)     , qstoif     , jpi, jpj, npac(1:nbpac) ) 
    400          CALL tab_2d_1d_2( nbpac, rdmicif_1d(1:nbpac)     , rdmicif    , jpi, jpj, npac(1:nbpac) ) 
     397         CALL tab_2d_1d_2( nbpac, rdm_ice_1d(1:nbpac)     , rdm_ice    , jpi, jpj, npac(1:nbpac) ) 
     398         CALL tab_2d_1d_2( nbpac, rdq_ice_1d(1:nbpac)     , rdq_ice    , jpi, jpj, npac(1:nbpac) ) 
    401399         CALL tab_2d_1d_2( nbpac, dvlbq_1d  (1:nbpac)     , zdvolif    , jpi, jpj, npac(1:nbpac) ) 
    402400         CALL tab_2d_1d_2( nbpac, tfu_1d    (1:nbpac)     , tfu        , jpi, jpj, npac(1:nbpac) ) 
     
    412410         CALL tab_1d_2d_2( nbpac, tbif(:,:,3), npac(1:nbpac), tbif_1d   (1:nbpac , 3 ), jpi, jpj ) 
    413411         CALL tab_1d_2d_2( nbpac, qstoif     , npac(1:nbpac), qstbif_1d (1:nbpac)     , jpi, jpj ) 
    414          CALL tab_1d_2d_2( nbpac, rdmicif    , npac(1:nbpac), rdmicif_1d(1:nbpac)     , jpi, jpj ) 
     412         CALL tab_1d_2d_2( nbpac, rdm_ice    , npac(1:nbpac), rdm_ice_1d(1:nbpac)     , jpi, jpj ) 
     413         CALL tab_1d_2d_2( nbpac, rdq_ice    , npac(1:nbpac), rdq_ice_1d(1:nbpac)     , jpi, jpj ) 
    415414         CALL tab_1d_2d_2( nbpac, zdvolif    , npac(1:nbpac), dvlbq_1d  (1:nbpac)     , jpi, jpj ) 
    416415         ! 
     
    443442      CALL iom_put( 'iceprod_cea' , hicifp (:,:) * zztmp     )   ! Ice produced               [m/s] 
    444443      IF( lk_diaar5 ) THEN 
    445          CALL iom_put( 'snowmel_cea' , rdmsnif(:,:) * zztmp     )   ! Snow melt                  [kg/m2/s] 
     444         CALL iom_put( 'snowmel_cea' , rdm_snw(:,:) * zztmp     )   ! Snow melt                  [kg/m2/s] 
    446445         zztmp = rhoic / rdt_ice 
    447446         CALL iom_put( 'sntoice_cea' , zdvonif(:,:) * zztmp     )   ! Snow to Ice transformation [kg/m2/s] 
    448447         CALL iom_put( 'ticemel_cea' , zdvosif(:,:) * zztmp     )   ! Melt at Sea Ice top        [kg/m2/s] 
    449448         CALL iom_put( 'bicemel_cea' , zdvomif(:,:) * zztmp     )   ! Melt at Sea Ice bottom     [kg/m2/s] 
    450          zlicegr(:,:) = MAX( 0.e0, rdmicif(:,:)-zlicegr(:,:) ) 
     449         zlicegr(:,:) = MAX( 0.e0, rdm_ice(:,:)-zlicegr(:,:) ) 
    451450         CALL iom_put( 'licepro_cea' , zlicegr(:,:) * zztmp     )   ! Latereal sea ice growth    [kg/m2/s] 
    452451      ENDIF 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_lac_2.F90

    r1857 r1858  
    11MODULE limthd_lac_2 
    2 #if defined key_lim2 
    32   !!====================================================================== 
    43   !!                       ***  MODULE limthd_lac_2   *** 
    54   !!                lateral thermodynamic growth of the ice  
    65   !!====================================================================== 
    7  
     6   !! History : LIM  ! 2000-04  (UCL) Original code 
     7   !!           2.0  ! 2002-07  (C. Ethe, G. Madec) F90 & mpp 
     8   !!            -   ! 2006-07  (G. Madec) surface module 
     9   !!           2.1  ! 2010-05  (Y. Aksenov, M. Vancoppenolle, G. Madec) add heat content exchanges 
     10   !!---------------------------------------------------------------------- 
     11#if defined key_lim2 
    812   !!---------------------------------------------------------------------- 
    913   !!   lim_lat_acr_2    : lateral accretion of ice 
    10    !! * Modules used 
     14   !!---------------------------------------------------------------------- 
    1115   USE par_oce          ! ocean parameters 
    1216   USE phycst 
     
    6064      !!               update h_snow_1d, h_ice_1d and tbif_1d(:,:)       
    6165      !!  
    62       !! ** References : 
    63       !!      M. Maqueda, 1995, PhD Thesis, Univesidad Complutense Madrid 
    64       !!      Fichefet T. and M. Maqueda 1997, J. Geo. Res., 102(C6),  
    65       !!                                                12609 -12646    
    66       !! History : 
    67       !!   1.0  !  01-04 (LIM)  original code 
    68       !!   2.0  !  02-08 (C. Ethe, G. Madec)  F90, mpp 
     66      !! References : M. Maqueda, 1995, PhD Thesis, Univesidad Complutense Madrid 
     67      !!              Fichefet T. and M. Maqueda 1997, JGR, 102(C6), 12609 -12646    
    6968      !!------------------------------------------------------------------- 
    7069      !! * Arguments 
     
    144143         frld_1d   (ji) = MAX( zfrlnew , zfrlmin(ji) ) 
    145144         !--computation of the remaining part of ice thickness which has been already used 
    146          zdhicbot(ji) =  ( frld_1d(ji) - zfrlnew ) * zhice0(ji) / ( 1.0 - zfrlmin(ji) ) &  
    147                       -  (  ( 1.0 - zfrrate ) / ( 1.0 - frld_1d(ji) ) )  * ( zqbgow(ji) / xlic )  
     145         zdhicbot(ji) =  ( frld_1d(ji) - zfrlnew ) * zhice0(ji) / ( 1.0 - zfrlmin(ji) )   &  
     146            &         -  (  ( 1.0 - zfrrate ) / ( 1.0 - frld_1d(ji) ) )  * ( zqbgow(ji) / xlic )  
    148147      END DO 
    149148  
     
    195194            &          ) / zah 
    196195          
    197          tbif_1d(ji,3) =     (  iiceform * ( zhnews2 - zdh3 )                                          * zta1  & 
     196         tbif_1d(ji,3) =     ( iiceform * ( zhnews2 - zdh3 )                                           * zta1  & 
    198197            &              + ( iiceform * zdh3 + ( 1 - iiceform ) * zdh1 )                             * zta2  & 
    199198            &              + ( iiceform * ( zhnews2 - zdh5 ) + ( 1 - iiceform ) * ( zhnews2 - zdh1 ) ) * zta3  &  
     
    213212      !           dV = Vnew - Vold 
    214213      !------------------------------------------------------------- 
    215        
    216214      DO ji = kideb , kiut 
    217215         dvlbq_1d  (ji) = ( 1. - frld_1d(ji) ) * h_ice_1d(ji) - ( 1. - zfrl_old(ji) ) * zhice_old(ji) 
    218          rdmicif_1d(ji) = rdmicif_1d(ji) + rhoic * dvlbq_1d(ji) 
     216         rdm_ice_1d(ji) = rdm_ice_1d(ji) + rhoic * dvlbq_1d(ji) 
     217         rdq_ice_1d(ji) = rdq_ice_1d(ji) + rcpic * dvlbq_1d(ji) * ( tfu_1d(ji) - rt0 )      ! heat content relative to rt0 
    219218      END DO 
    220219       
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_zdf_2.F90

    r1857 r1858  
    44   !!                thermodynamic growth and decay of the ice  
    55   !!====================================================================== 
    6    !! History :  1.0  !  01-04 (LIM) Original code 
    7    !!            2.0  !  02-08 (C. Ethe, G. Madec) F90 
     6   !! History : LIM  ! 2001-04 (UCL) Original code 
     7   !!           2.0  ! 2002-08 (C. Ethe, G. Madec) re-writing F90 
     8   !!           2.1  ! 2010-05  (Y. Aksenov, M. Vancoppenolle, G. Madec) add heat content exchanges 
    89   !!---------------------------------------------------------------------- 
    910#if defined key_lim2 
     
    158159          , zdrfrl1, zdrfrl2 &  ! tempory scalars 
    159160          , zihsn, zidhb, zihic, zihe, zihq, ziexp, ziqf, zihnf, zibmlt, ziqr, zihgnew, zind 
     161       REAL(wp) ::   ztmp   ! temporary scalar 
    160162       !!---------------------------------------------------------------------- 
    161163 
     
    168170        
    169171       DO ji = kideb , kiut 
     172          ! do nothing if the snow (ice) thickness falls below its minimum thickness 
    170173          zihsn = MAX( zzero , SIGN( zone , hsndif - h_snow_1d(ji) ) ) 
    171           zihic = MAX( zzero , SIGN( zone , hicdif - h_ice_1d(ji) ) ) 
    172           !--computation of energy due to surface melting 
    173           zqcmlt(ji,1) = ( MAX ( zzero ,  & 
    174              &                   rcpsn * h_snow_1d(ji) * ( tbif_1d(ji,1) - rt0_snow ) ) ) * ( 1.0 - zihsn ) 
    175           !--computation of energy due to bottom melting 
    176           zqcmlt(ji,2) = ( MAX( zzero , & 
    177              &                  rcpic * ( tbif_1d(ji,2) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) & 
    178              &           + MAX( zzero , & 
    179              &                  rcpic * ( tbif_1d(ji,3) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) & 
     174          zihic = MAX( zzero , SIGN( zone , hicdif - h_ice_1d (ji) ) ) 
     175          !--energy required to bring snow to its melting point (rt0_snow) 
     176          zqcmlt(ji,1) = ( MAX ( zzero , rcpsn * h_snow_1d(ji) * ( tbif_1d(ji,1) - rt0_snow ) ) ) * ( 1.0 - zihsn ) 
     177          !--energy required to bring ice to its melting point (rt0_ice) 
     178          zqcmlt(ji,2) = ( MAX( zzero , rcpic * ( tbif_1d(ji,2) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) )   & 
     179             &           + MAX( zzero , rcpic * ( tbif_1d(ji,3) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) )   & 
    180180             &           ) * ( 1.0 - zihic  ) 
    181           !--limitation of  snow/ice system internal temperature 
     181          !--limitation of snow/ice system internal temperature 
    182182          tbif_1d(ji,1)   = MIN( rt0_snow, tbif_1d(ji,1) ) 
    183183          tbif_1d(ji,2)   = MIN( rt0_ice , tbif_1d(ji,2) ) 
     
    478478          dvsbq_1d(ji) =  ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsnw_old(ji) - zsprecip(ji) ) 
    479479          dvsbq_1d(ji) =  MIN( zzero , dvsbq_1d(ji) ) 
    480           rdmsnif_1d(ji) =  rhosn * dvsbq_1d(ji) 
     480          ztmp = rhosn * dvsbq_1d(ji) 
     481          rdm_snw_1d(ji) =  ztmp 
     482          !--heat content of the water provided to the ocean (referenced to rt0) 
     483          rdq_snw_1d(ji) =  cpic * ztmp * ( rt0_snow - rt0 ) 
    481484          !-- If the snow is completely melted the remaining heat is used to melt ice 
    482485          zqsn_mlt_rem  = MAX( zzero , -zhsn ) * xlsn 
     
    621624          !---updating new ice thickness and computing the newly formed ice mass 
    622625          zhicnew   =  zihgnew * zhicnew 
    623           rdmicif_1d(ji) =  rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d(ji) ) * rhoic 
     626          ztmp    =  ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d(ji) ) * rhoic 
     627          rdm_ice_1d(ji) =  rdm_ice_1d(ji) + ztmp 
     628          !---heat content of the water provided to the ocean (referenced to rt0) 
     629          ! use of rt0_ice is OK for melting ice, in case of freezing tfu_1d should be used. This is done in 9.5 section (see below) 
     630          rdq_ice_1d(ji) =  cpic * ztmp * ( rt0_ice - rt0 ) 
    624631          !---updating new snow thickness and computing the newly formed snow mass 
    625632          zhsnfi   = zhsn + zdhsnm 
    626633          h_snow_1d(ji) = MAX( zzero , zhsnfi ) 
    627           rdmsnif_1d(ji) =  rdmsnif_1d(ji) + ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsn ) * rhosn 
     634          ztmp = ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsn ) * rhosn 
     635          rdm_snw_1d(ji) = rdm_snw_1d(ji) + ztmp 
     636          !---updating the heat content of the water provided to the ocean (referenced to rt0) 
     637          rdq_snw_1d(ji) = rdq_snw_1d(ji) + cpic * ztmp * ( rt0_snow - rt0 ) 
    628638          !--remaining energy in case of total ablation 
    629639          zqocea(ji) = - ( zihsn * xlic * zdhicm + xlsn * ( zhsnfi - h_snow_1d(ji) ) ) * ( 1.0 - frld_1d(ji) ) 
     
    657667          tbif_1d(ji,3) =  zihgnew * ztb3 + ( 1.0 - zihgnew ) * tfu_1d(ji) 
    658668          h_ice_1d(ji)  =  zhicnew 
    659        END DO 
     669          ! update the ice heat content given to the ocean in freezing case (part from rt0_ice to tfu_1d) 
     670         ztmp = ( 1. - zidhb ) * rhoic * dvbbq_1d(ji) 
     671         rdqicif_1d(ji) = rdqicif_1d(ji) + cpic * ztmp * ( tfu_1d(ji) - rt0_ice ) 
     672      END DO 
    660673 
    661674 
     
    698711          dmgwi_1d(ji) = dmgwi_1d(ji) + ( 1.0 -frld_1d(ji) ) * ( h_snow_1d(ji) - zhsnnew ) * rhosn 
    699712          !---  volume change of ice and snow (used for ocean-ice freshwater flux computation) 
    700           rdmicif_1d(ji) = rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhicnew - h_ice_1d (ji) ) * rhoic 
    701           rdmsnif_1d(ji) = rdmsnif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhsnnew - h_snow_1d(ji) ) * rhosn 
     713          ztmp = ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d (ji) ) * rhoic 
     714          rdm_ice_1d(ji) = rdm_ice_1d(ji) + ztmp 
     715          rdq_ice_1d(ji) = rdq_ice_1d(ji) + cpic * ztmp * ( tfu_1d(ji) - rt0 ) 
     716!!gm BUG ??   snow ==>  only needed for nn_ice_embd == 0  (standard levitating sea-ice) 
     717          ztmp = ( 1.0 - frld_1d(ji) )   * ( zhsnnew - h_snow_1d(ji) ) * rhosn 
     718         rdm_snw_1d(ji) = rdm_snw_1d(ji) + ztmp 
     719         rdq_snw_1d(ji) = rdq_snw_1d(ji) + cpic * ztmp * ( rt0_snow - rt0 ) 
    702720 
    703721          !---  Actualize new snow and ice thickness. 
     
    746764          !--variation of ice volume and ice mass  
    747765          dvlbq_1d(ji)   = zihic * ( zfrl_old(ji) - frld_1d(ji) ) * h_ice_1d(ji) 
    748           rdmicif_1d(ji) = rdmicif_1d(ji) + dvlbq_1d(ji) * rhoic 
     766          ztmp = dvlbq_1d(ji) * rhoic 
     767          rdm_ice_1d(ji) = rdm_ice_1d(ji) + ztmp 
     768!!gm 
     769!!gm   This should be split in two parts: 
     770!!gm         1-  heat required to bring sea-ice at tfu  : this part should be added to the heat flux taken from the ocean 
     771!!gm                 cpic * ztmp * 0.5 * ( tbif_1d(ji,2) + tbif_1d(ji,3) - 2.* rt0_ice ) 
     772!!gm         2-  heat content of lateral ablation referenced to rt0 : this part only put in rdq_ice_1d 
     773!!gm                 cpic * ztmp * ( rt0_ice - rt0 ) 
     774!!gm   Currently we put all the heat in rdq_ice_1d 
     775          rdq_ice_1d(ji) = rdq_ice_1d(ji) + cpic * ztmp * 0.5 * ( tbif_1d(ji,2) + tbif_1d(ji,3) - 2.* rt0 ) 
     776          ! 
    749777          !--variation of snow volume and snow mass  
    750           zdvsnvol    = zihsn * ( zfrl_old(ji) - frld_1d(ji) ) * h_snow_1d(ji) 
    751           rdmsnif_1d(ji) = rdmsnif_1d(ji) + zdvsnvol * rhosn 
     778          zdvsnvol = zihsn * ( zfrl_old(ji) - frld_1d(ji) ) * h_snow_1d(ji) 
     779          ztmp     = zdvsnvol * rhosn 
     780          rdm_snw_1d(ji) = rdm_snw_1d(ji) + ztmp 
     781!!gm 
     782!!gm   This should be split in two parts: 
     783!!gm         1-  heat required to bring snow at tfu  : this part should be added to the heat flux taken from the ocean 
     784!!gm                 cpic * ztmp * ( tbif_1d(ji,1) - rt0_snow ) 
     785!!gm         2-  heat content of lateral ablation referenced to rt0 : this part only put in rdqicif_1d 
     786!!gm                 cpic * ztmp * ( rt0_snow - rt0 ) 
     787!!gm   Currently we put all the heat in rdqicif_1d 
     788          rdq_snw_1d(ji) = rdq_snw_1d(ji) + cpic * ztmp * ( tbif_1d(ji,1) - rt0 ) 
     789 
    752790          h_snow_1d(ji)  = ziqf * h_snow_1d(ji) 
    753791 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/thd_ice_2.F90

    r1857 r1858  
    55   !! LIM 2.0 sea-ice :   Ice thermodynamics in 1D 
    66   !!===================================================================== 
    7    !! History : 
    8    !!   2.0  !  02-11  (C. Ethe)  F90: Free form and module 
     7   !! History : 2.0  ! 2002-11  (C. Ethe)  F90: Free form and module 
     8   !!           2.1  ! 2010-05  (Y. Aksenov, M. Vancoppenolle, G. Madec) add heat content exchanges 
    99   !!---------------------------------------------------------------------- 
    1010   !!   LIM 2.0, UCL-LOCEAN-IPSL (2005) 
     
    6767      qstbif_1d   ,     &  !:    "                  "      qstoif 
    6868      fbif_1d     ,     &  !:    "                  "      fbif 
    69       rdmicif_1d  ,     &  !:    "                  "      rdmicif 
    70       rdmsnif_1d  ,     &  !:    "                  "      rdmsnif 
     69      rdm_ice_1d  ,     &  !:    "                  "      rdm_ice 
     70      rdq_ice_1d  ,     &  !:    "                  "      rdq_ice 
     71      rdm_snw_1d  ,     &  !:    "                  "      rdm_snw 
     72      rdq_snw_1d  ,     &  !:    "                  "      rdq_snw 
    7173      qlbbq_1d    ,     &  !:    "                  "      qlbsbq 
    7274      dmgwi_1d    ,     &  !:    "                  "      dmgwi 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/DOM/phycst.F90

    r1739 r1858  
    44   !!     Definition of of both ocean and ice parameters used in the code 
    55   !!===================================================================== 
    6    !! History :        !  90-10  (C. Levy - G. Madec)  Original code 
    7    !!                  !  91-11  (G. Madec) 
    8    !!                  !  91-12  (M. Imbard) 
    9    !!             8.5  !  02-08  (G. Madec, C. Ethe)  F90, add ice constants 
    10    !!             9.0  !  06-08  (G. Madec)  style  
     6   !! History :  OPA  !  1990-10  (C. Levy - G. Madec)  Original code 
     7   !!                 !  1991-11  (G. Madec, M. Imbard) 
     8   !!  NEMO      1.0  !  2002-08  (G. Madec, C. Ethe)  F90, add ice constants 
     9   !!             -   !  2006-08  (G. Madec)  style  
     10   !!            3.0  !  2008-08  (G. Madec)  add LIM-3  
    1111   !!---------------------------------------------------------------------- 
    1212 
     
    2626   REAL(wp), PUBLIC ::   rsmall = 0.5 * EPSILON( 1. )           !: smallest real computer value 
    2727    
    28    REAL(wp), PUBLIC ::          & !: 
    29       rday = 24.*60.*60.  ,     & !: day (s) 
    30       rsiyea              ,     & !: sideral year (s) 
    31       rsiday              ,     & !: sideral day (s) 
    32       raamo =  12._wp     ,     & !: number of months in one year 
    33       rjjhh =  24._wp     ,     & !: number of hours in one day 
    34       rhhmm =  60._wp     ,     & !: number of minutes in one hour 
    35       rmmss =  60._wp     ,     & !: number of seconds in one minute 
    36 !!!   omega = 7.292115083046061e-5_wp ,  &  !: change the last digit! 
    37       omega               ,    &  !: earth rotation parameter 
    38       ra    = 6371229._wp ,    &  !: earth radius (meter) 
    39       grav  = 9.80665_wp          !: gravity (m/s2) 
     28   REAL(wp), PUBLIC ::   rday = 24.*60.*60.     !: day                                [s] 
     29   REAL(wp), PUBLIC ::   rsiyea                 !: sideral year                       [s] 
     30   REAL(wp), PUBLIC ::   rsiday                 !: sideral day                        [s] 
     31   REAL(wp), PUBLIC ::   raamo =  12._wp        !: number of months in one year 
     32   REAL(wp), PUBLIC ::   rjjhh =  24._wp        !: number of hours in one day 
     33   REAL(wp), PUBLIC ::   rhhmm =  60._wp        !: number of minutes in one hour 
     34   REAL(wp), PUBLIC ::   rmmss =  60._wp        !: number of seconds in one minute 
     35   REAL(wp), PUBLIC ::   omega                  !: earth rotation parameter           [s-1] 
     36   REAL(wp), PUBLIC ::   ra    = 6371229._wp    !: earth radius                       [m] 
     37   REAL(wp), PUBLIC ::   grav  = 9.80665_wp     !: gravity                            [m/s2] 
    4038    
    41    REAL(wp), PUBLIC ::         &  !: 
    42       rtt      = 273.16_wp  ,  &  !: triple point of temperature (Kelvin) 
    43       rt0      = 273.15_wp  ,  &  !: freezing point of water (Kelvin) 
    44 #if defined key_lim3 
    45       rt0_snow = 273.16_wp  ,  &  !: melting point of snow  (Kelvin) 
    46       rt0_ice  = 273.16_wp  ,  &  !: melting point of ice   (Kelvin) 
    47 #else 
    48       rt0_snow = 273.15_wp  ,  &  !: melting point of snow  (Kelvin) 
    49       rt0_ice  = 273.05_wp  ,  &  !: melting point of ice   (Kelvin) 
    50 #endif 
    51       rau0     = 1035._wp   ,  &  !: volumic mass of reference (kg/m3) 
    52       rauw     = 1000._wp   ,  &  !: volumic mass of pure water (kg/m3) 
    53       rcp      =    4.e+3_wp,  &  !: ocean specific heat 
    54       ro0cpr                      !: = 1. / ( rau0 * rcp ) 
    55  
    56    REAL(wp), PUBLIC ::            &  !: 
    57 #if defined key_lim3 
    58       rcdsn   =   0.31_wp     ,   &  !: thermal conductivity of snow 
    59       rcdic   =   2.034396_wp ,   &  !: thermal conductivity of fresh ice 
    60       cpic    = 2067.0        ,   & 
    61       ! add the following lines 
    62       lsub    = 2.834e+6      ,   &  !: pure ice latent heat of sublimation (J.kg-1) 
    63       lfus    = 0.334e+6      ,   &  !: latent heat of fusion of fresh ice   (J.kg-1) 
    64       rhoic   = 917._wp       ,   &  !: volumic mass of sea ice (kg/m3) 
    65       tmut    =   0.054       ,   &  !: decrease of seawater meltpoint with salinity 
    66 #else 
    67       rcdsn   =   0.22_wp     ,   &  !: conductivity of the snow 
    68       rcdic   =   2.034396_wp ,   &  !: conductivity of the ice 
    69       rcpsn   =   6.9069e+5_wp,   &  !: density times specific heat for snow 
    70       rcpic   =   1.8837e+6_wp,   &  !: volumetric latent heat fusion of sea ice 
    71       xlsn    = 110.121e+6_wp ,   &  !: volumetric latent heat fusion of snow 
    72       xlic    = 300.33e+6_wp  ,   &  !: volumetric latent heat fusion of ice 
    73       xsn     =   2.8e+6      ,   &  !: latent heat of sublimation of snow 
    74       rhoic   = 900._wp       ,   &  !: volumic mass of sea ice (kg/m3) 
    75 #endif 
    76       rhosn   = 330._wp       ,   &  !: volumic mass of snow (kg/m3) 
    77       emic    =   0.97_wp     ,   &  !: emissivity of snow or ice 
    78       sice    =   6.0_wp      ,   &  !: salinity of ice (psu) 
    79       soce    =  34.7_wp      ,   &  !: salinity of sea (psu) 
    80       cevap   =   2.5e+6_wp   ,   &  !: latent heat of evaporation (water) 
    81       srgamma =   0.9_wp      ,   &  !: correction factor for solar radiation (Oberhuber, 1974) 
    82       vkarmn  =   0.4_wp      ,   &  !: von Karman constant 
    83       stefan  =   5.67e-8_wp         !: Stefan-Boltzmann constant  
    84       !!---------------------------------------------------------------------- 
    85       !!  OPA 9.0 , LOCEAN-IPSL (2005)  
    86       !! $Id$  
    87       !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
    88       !!---------------------------------------------------------------------- 
    89     
     39   REAL(wp), PUBLIC ::   rtt      = 273.16_wp        !: triple point of temperature   [Kelvin] 
     40   REAL(wp), PUBLIC ::   rt0      = 273.15_wp        !: freezing point of fresh water [Kelvin] 
     41#if defined key_lim3 
     42   REAL(wp), PUBLIC ::   rt0_snow = 273.16_wp        !: melting point of snow         [Kelvin] 
     43   REAL(wp), PUBLIC ::   rt0_ice  = 273.16_wp        !: melting point of ice          [Kelvin] 
     44#else 
     45   REAL(wp), PUBLIC ::   rt0_snow = 273.15_wp        !: melting point of snow         [Kelvin] 
     46   REAL(wp), PUBLIC ::   rt0_ice  = 273.05_wp        !: melting point of ice          [Kelvin] 
     47#endif 
     48   REAL(wp), PUBLIC ::   rau0     = 1035._wp         !: volumic mass of reference     [kg/m3] 
     49   REAL(wp), PUBLIC ::   r1_rau0                     !: = 1. / rau0                   [m3/kg] 
     50   REAL(wp), PUBLIC ::   rauw     = 1000._wp         !: volumic mass of pure water    [m3/kg] 
     51   REAL(wp), PUBLIC ::   rcp      =    4.e3_wp       !: ocean specific heat           [J/Kelvin] 
     52   REAL(wp), PUBLIC ::   r1_rcp                      !: = 1. / rcp                    [Kelvin/J] 
     53   REAL(wp), PUBLIC ::   r1_rau0_rcp                 !: = 1. / ( rau0 * rcp ) 
     54 
     55   REAL(wp), PUBLIC ::   rhosn    =  330._wp         !: volumic mass of snow          [kg/m3] 
     56   REAL(wp), PUBLIC ::   emic     =    0.97_wp       !: emissivity of snow or ice 
     57   REAL(wp), PUBLIC ::   sice     =    6.0_wp        !: salinity of ice (psu) 
     58   REAL(wp), PUBLIC ::   soce     =   34.7_wp        !: salinity of sea (psu) 
     59 
     60#if defined key_lim3 
     61   REAL(wp), PUBLIC ::   rhoic    =  917._wp         !: volumic mass of sea ice                               [kg/m3] 
     62   REAL(wp), PUBLIC ::   rcdic    =    2.034396_wp   !: thermal conductivity of fresh ice 
     63   REAL(wp), PUBLIC ::   rcdsn    =    0.31_wp       !: thermal conductivity of snow 
     64   REAL(wp), PUBLIC ::   cpic     = 2067.0_wp        !: specific heat for ice  
     65   REAL(wp), PUBLIC ::   lsub     =    2.834e+6_wp   !: pure ice latent heat of sublimation                   [J/kg] 
     66   REAL(wp), PUBLIC ::   lfus     =    0.334e+6_wp   !: latent heat of fusion of fresh ice                    [J/kg] 
     67   REAL(wp), PUBLIC ::   tmut     =    0.054_wp      !: decrease of seawater meltpoint with salinity 
     68   REAL(wp), PUBLIC ::   xlsn                        !: = lfus*rhosn (volumetric latent heat fusion of snow)  [J/m3] 
     69#else 
     70   REAL(wp), PUBLIC ::   rhoic    =  900._wp         !: volumic mass of sea ice                               [kg/m3] 
     71   REAL(wp), PUBLIC ::   rcdic    =    2.034396_wp   !: conductivity of the ice                               [W/m/K] 
     72   REAL(wp), PUBLIC ::   rcpic    =    1.8837e+6_wp  !: volumetric specific heat for ice                      [J/m3/K] 
     73   REAL(wp), PUBLIC ::   cpic                        !: = rcpic / rhoic  (specific heat for ice)              [J/Kg/K] 
     74   REAL(wp), PUBLIC ::   rcdsn    =    0.22_wp       !: conductivity of the snow                              [W/m/K] 
     75   REAL(wp), PUBLIC ::   rcpsn    =    6.9069e+5_wp  !: volumetric specific heat for snow                     [J/m3/K] 
     76   REAL(wp), PUBLIC ::   xlsn     =  110.121e+6_wp   !: volumetric latent heat fusion of snow                 [J/m3] 
     77   REAL(wp), PUBLIC ::   lfus                        !: = xlsn / rhosn   (latent heat of fusion of fresh ice) [J/Kg] 
     78   REAL(wp), PUBLIC ::   xlic     =  300.33e+6_wp    !: volumetric latent heat fusion of ice                  [J/m3] 
     79   REAL(wp), PUBLIC ::   xsn      =    2.8e+6_wp     !: volumetric latent heat of sublimation of snow         [J/m3] 
     80#endif 
     81   REAL(wp), PUBLIC ::   cevap    =    2.5e+6_wp     !: latent heat of evaporation (water) 
     82   REAL(wp), PUBLIC ::   srgamma  =    0.9_wp        !: correction factor for solar radiation (Oberhuber, 1974) 
     83   REAL(wp), PUBLIC ::   vkarmn   =    0.4_wp        !: von Karman constant 
     84   REAL(wp), PUBLIC ::   stefan   =    5.67e-8_wp    !: Stefan-Boltzmann constant  
     85   !!---------------------------------------------------------------------- 
     86   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)  
     87   !! $Id$  
     88   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     89   !!---------------------------------------------------------------------- 
    9090CONTAINS 
    9191    
     
    153153      IF(lwp) WRITE(numout,*) '          melting point of ice             rt0_ice  = ', rt0_ice , ' K' 
    154154 
    155       ro0cpr = 1. / ( rau0 * rcp ) 
    156       IF(lwp) WRITE(numout,*) 
    157       IF(lwp) WRITE(numout,*) '          volumic mass of pure water         rauw   = ', rauw, ' kg/m^3' 
    158       IF(lwp) WRITE(numout,*) '          volumic mass of reference          rau0   = ', rau0, ' kg/m^3' 
    159       IF(lwp) WRITE(numout,*) '          ocean specific heat                rcp    = ', rcp 
    160       IF(lwp) WRITE(numout,*) '                       1. / ( rau0 * rcp ) = ro0cpr = ', ro0cpr 
     155      r1_rau0     = 1.e0 / rau0 
     156      r1_rcp      = 1.e0 / rcp 
     157      r1_rau0_rcp = 1.e0 / ( rau0 * rcp ) 
     158      IF(lwp) WRITE(numout,*) 
     159      IF(lwp) WRITE(numout,*) '          volumic mass of pure water          rauw  = ', rauw   , ' kg/m3' 
     160      IF(lwp) WRITE(numout,*) '          volumic mass of reference           rau0  = ', rau0   , ' kg/m3' 
     161      IF(lwp) WRITE(numout,*) '          1. / rau0                        r1_rau0  = ', r1_rau0, ' m3/kg' 
     162      IF(lwp) WRITE(numout,*) '          ocean specific heat                 rcp   = ', rcp    , ' J/Kelvin' 
     163      IF(lwp) WRITE(numout,*) '          1. / ( rau0 * rcp )           r1_rau0_rcp = ', r1_rau0_rcp 
     164 
     165 
     166#if defined key_lim3 
     167      xlsn = lfus * rhosn        ! volumetric latent heat fusion of snow [J/m3] 
     168#else 
     169      cpic = rcpic / rhoic       ! specific heat for ice   [J/Kg/K] 
     170      lfus = xlsn / rhosn        ! latent heat of fusion of fresh ice 
     171#endif 
    161172 
    162173      IF(lwp) THEN 
     
    164175         WRITE(numout,*) '          thermal conductivity of the snow          = ', rcdsn   , ' J/s/m/K' 
    165176         WRITE(numout,*) '          thermal conductivity of the ice           = ', rcdic   , ' J/s/m/K' 
    166 #if defined key_lim3 
    167177         WRITE(numout,*) '          fresh ice specific heat                   = ', cpic    , ' J/kg/K' 
    168178         WRITE(numout,*) '          latent heat of fusion of fresh ice / snow = ', lfus    , ' J/kg' 
     179#if defined key_lim3 
    169180         WRITE(numout,*) '          latent heat of subl.  of fresh ice / snow = ', lsub    , ' J/kg' 
    170181#else 
    171          WRITE(numout,*) '          density times specific heat for snow      = ', rcpsn   , ' J/m^3/K'  
    172          WRITE(numout,*) '          density times specific heat for ice       = ', rcpic   , ' J/m^3/K' 
     182         WRITE(numout,*) '          density times specific heat for snow      = ', rcpsn   , ' J/m3/K'  
     183         WRITE(numout,*) '          density times specific heat for ice       = ', rcpic   , ' J/m3/K' 
    173184         WRITE(numout,*) '          volumetric latent heat fusion of sea ice  = ', xlic    , ' J/m'  
    174          WRITE(numout,*) '          volumetric latent heat fusion of snow     = ', xlsn    , ' J/m'  
    175185         WRITE(numout,*) '          latent heat of sublimation of snow        = ', xsn     , ' J/kg'  
    176186#endif 
    177          WRITE(numout,*) '          density of sea ice                        = ', rhoic   , ' kg/m^3' 
    178          WRITE(numout,*) '          density of snow                           = ', rhosn   , ' kg/m^3' 
     187         WRITE(numout,*) '          volumetric latent heat fusion of snow     = ', xlsn    , ' J/m^3'  
     188         WRITE(numout,*) '          density of sea ice                        = ', rhoic   , ' kg/m3' 
     189         WRITE(numout,*) '          density of snow                           = ', rhosn   , ' kg/m3' 
    179190         WRITE(numout,*) '          emissivity of snow or ice                 = ', emic   
    180191         WRITE(numout,*) '          salinity of ice                           = ', sice    , ' psu' 
    181192         WRITE(numout,*) '          salinity of sea                           = ', soce    , ' psu' 
    182          WRITE(numout,*) '          latent heat of evaporation (water)        = ', cevap   , ' J/m^3'  
     193         WRITE(numout,*) '          latent heat of evaporation (water)        = ', cevap   , ' J/m3'  
    183194         WRITE(numout,*) '          correction factor for solar radiation     = ', srgamma  
    184195         WRITE(numout,*) '          von Karman constant                       = ', vkarmn  
     
    191202         WRITE(numout,*) '          smallest real computer value       rsmall = ', rsmall 
    192203      ENDIF 
    193  
     204      ! 
    194205   END SUBROUTINE phy_cst 
    195206 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/DYN/dynzdf_exp.F90

    r1537 r1858  
    5454      !! ** Action : - Update (ua,va) with the vertical diffusive trend 
    5555      !!--------------------------------------------------------------------- 
    56       !! * Arguments 
    57       INTEGER , INTENT( in ) ::   kt                           ! ocean time-step index 
    58       REAL(wp), INTENT( in ) ::   p2dt                         ! time-step  
    59  
    60       !! * Local declarations 
    61       INTEGER ::   ji, jj, jk, jl                              ! dummy loop indices 
    62       REAL(wp) ::   zrau0r, zlavmr, zua, zva                   ! temporary scalars 
    63       REAL(wp), DIMENSION(jpi,jpk) ::   zwx, zwy, zwz, zww     ! temporary workspace arrays 
     56      INTEGER , INTENT( in ) ::   kt     ! ocean time-step index 
     57      REAL(wp), INTENT( in ) ::   p2dt   ! time-step  
     58      !!  
     59      INTEGER  ::   ji, jj, jk, jl     ! dummy loop indices 
     60      REAL(wp) ::   zlavmr, zua, zva   ! temporary scalars 
     61      REAL(wp), DIMENSION(jpi,jpk) ::   zwx, zwy, zwz, zww   ! temporary workspace arrays 
    6462      !!---------------------------------------------------------------------- 
    6563 
     
    7270      ! Local constant initialization 
    7371      ! ----------------------------- 
    74       zrau0r = 1. / rau0                                   ! inverse of the reference density 
    7572      zlavmr = 1. / float( nn_zdfexp )                      ! inverse of the number of sub time step 
    7673 
     
    8178         ! Surface boundary condition 
    8279         DO ji = 2, jpim1 
    83             zwy(ji,1) = utau(ji,jj) * zrau0r 
    84             zww(ji,1) = vtau(ji,jj) * zrau0r 
     80            zwy(ji,1) = utau(ji,jj) * r1_rau0 
     81            zww(ji,1) = vtau(ji,jj) * r1_rau0 
    8582         END DO   
    8683 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r1662 r1858  
    5959      !!         ua = ua + dz( avmu dz(u) ) 
    6060      !! 
    61       !! ** Action : - Update (ua,va) arrays with the after vertical diffusive 
    62       !!               mixing trend. 
     61      !! ** Action : -  (ua,va) updated with the after vertical diffusive mixing trend 
    6362      !!--------------------------------------------------------------------- 
    64       !! * Modules used 
    65       USE oce, ONLY :  zwd   => ta,   &                ! use ta as workspace 
    66                        zws   => sa                     ! use sa as workspace 
    67  
    68       !! * Arguments 
    69       INTEGER , INTENT( in ) ::   kt                   ! ocean time-step index 
    70       REAL(wp), INTENT( in ) ::  p2dt                  ! vertical profile of tracer time-step 
    71  
    72       !! * Local declarations 
    73       INTEGER ::   ji, jj, jk                          ! dummy loop indices 
    74       REAL(wp) ::   zrau0r, z2dtf, zcoef, zzws, zrhs   ! temporary scalars 
    75       REAL(wp) ::   zzwi                               ! temporary scalars 
    76       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi        ! temporary workspace arrays 
     63      USE oce, ONLY :  zwd   => ta   ! use ta as workspace 
     64      USE oce, ONLY :  zws   => sa   !  -  sa         - 
     65      !! 
     66      INTEGER , INTENT( in ) ::   kt    ! ocean time-step index 
     67      REAL(wp), INTENT( in ) ::  p2dt   ! vertical profile of tracer time-step 
     68      !! 
     69      INTEGER ::   ji, jj, jk                     ! dummy loop indices 
     70      REAL(wp) ::   z2dtf, zcoef, zzws, zrhs      ! temporary scalars 
     71      REAL(wp) ::   zzwi                          !    -         - 
     72      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi   ! 3D workspace  
    7773      !!---------------------------------------------------------------------- 
    7874 
     
    8278         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
    8379      ENDIF 
    84  
    85       ! 0. Local constant initialization 
    86       ! -------------------------------- 
    87       zrau0r = 1. / rau0      ! inverse of the reference density 
    8880 
    8981      ! 1. Vertical diffusion on u 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/TRA/eosbn2.F90

    r1613 r1858  
    114114      REAL(wp) ::   zd , zc , zaw, za    !    -         - 
    115115      REAL(wp) ::   zb1, za1, zkw, zk0   !    -         - 
    116       REAL(wp) ::   zrau0r               !    -         - 
    117116      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zws   ! temporary workspace 
    118117      !!---------------------------------------------------------------------- 
     
    121120      ! 
    122121      CASE( 0 )                !==  Jackett and McDougall (1994) formulation  ==! 
    123          zrau0r = 1.e0 / rau0 
    124122!CDIR NOVERRCHK 
    125123         zws(:,:,:) = SQRT( ABS( psal(:,:,:) ) ) 
     
    162160                  ! masked in situ density anomaly 
    163161                  prd(ji,jj,jk) = (  zrhop / (  1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )  )    & 
    164                      &             - rau0  ) * zrau0r * tmask(ji,jj,jk) 
     162                     &             - rau0  ) * r1_rau0 * tmask(ji,jj,jk) 
    165163               END DO 
    166164            END DO 
     
    237235      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    238236      REAL(wp) ::   zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw   ! temporary scalars 
    239       REAL(wp) ::   zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zrau0r       !    -         - 
     237      REAL(wp) ::   zb, zd, zc, zaw, za, zb1, za1, zkw, zk0               !    -         - 
    240238      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zws   ! 3D workspace 
    241239      !!---------------------------------------------------------------------- 
     
    244242      ! 
    245243      CASE( 0 )                !==  Jackett and McDougall (1994) formulation  ==! 
    246          zrau0r = 1.e0 / rau0 
    247244!CDIR NOVERRCHK 
    248245         zws(:,:,:) = SQRT( ABS( psal(:,:,:) ) ) 
     
    288285                  ! masked in situ density anomaly 
    289286                  prd(ji,jj,jk) = (  zrhop / (  1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )  )    & 
    290                      &             - rau0  ) * zrau0r * tmask(ji,jj,jk) 
     287                     &             - rau0  ) * r1_rau0 * tmask(ji,jj,jk) 
    291288               END DO 
    292289            END DO 
     
    296293         DO jk = 1, jpkm1 
    297294            prd  (:,:,jk) = ( 0.0285 - rn_alpha * ptem(:,:,jk) )        * tmask(:,:,jk) 
    298             prhop(:,:,jk) = ( 1.e0   +          prd (:,:,jk) ) * rau0 * tmask(:,:,jk) 
     295            prhop(:,:,jk) = ( 1.e0   +            prd (:,:,jk) ) * rau0 * tmask(:,:,jk) 
    299296         END DO 
    300297         ! 
     
    302299         DO jk = 1, jpkm1 
    303300            prd  (:,:,jk) = ( rn_beta  * psal(:,:,jk) - rn_alpha * ptem(:,:,jk) )        * tmask(:,:,jk) 
    304             prhop(:,:,jk) = ( 1.e0   + prd (:,:,jk) )                         * rau0 * tmask(:,:,jk) 
     301            prhop(:,:,jk) = ( 1.e0   + prd (:,:,jk) )                             * rau0 * tmask(:,:,jk) 
    305302         END DO 
    306303         ! 
     
    409406               ! 
    410407               ! masked in situ density anomaly 
    411                prd(ji,jj) = ( zrhop / (  1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )  ) - rau0 ) / rau0 * zmask 
     408               prd(ji,jj) = ( zrhop / (  1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )  ) - rau0 ) * r1_rau0 * zmask 
    412409            END DO 
    413410         END DO 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/TRA/trabbc.F90

    r1601 r1858  
    7171      !!---------------------------------------------------------------------- 
    7272      USE oce, ONLY :   ztrdt => ua   ! use ua as 3D workspace    
    73       USE oce, ONLY :   ztrds => va   ! use va as 3D workspace    
     73      USE oce, ONLY :   ztrds => va   !  -  va     -      -   
    7474      !! 
    7575      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     
    9898            DO ji = 2, jpim1 
    9999#endif 
    100                zqgh_trd = ro0cpr * qgh_trd0(ji,jj) / fse3t(ji,jj,nbotlevt(ji,jj)) 
     100               zqgh_trd = r1_rau0_rcp * qgh_trd0(ji,jj) / fse3t(ji,jj,nbotlevt(ji,jj)) 
    101101               ta(ji,jj,nbotlevt(ji,jj)) = ta(ji,jj,nbotlevt(ji,jj)) + zqgh_trd 
    102102            END DO 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/TRA/traqsr.F90

    r1756 r1858  
    120120            DO jj = 2, jpjm1 
    121121               DO ji = fs_2, fs_jpim1   ! vector opt. 
    122                   ta(ji,jj,jk) = ta(ji,jj,jk) + ro0cpr * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / fse3t(ji,jj,jk)  
     122                  ta(ji,jj,jk) = ta(ji,jj,jk) + r1_rau0_rcp * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / fse3t(ji,jj,jk) 
    123123               END DO 
    124124            END DO 
     
    177177               ! 
    178178               DO jk = 1, nksr                                        ! compute and add qsr trend to ta 
    179                   ta(:,:,jk) = ta(:,:,jk) + ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk) 
     179                  ta(:,:,jk) = ta(:,:,jk) + r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk) 
    180180               END DO 
    181181               zea(:,:,nksr+1:jpk) = 0.e0     ! below 400m set to zero 
     
    377377               ! 
    378378               DO jk = 1, nksr 
    379                   etot3(:,:,jk) = ro0cpr * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk) 
     379                  etot3(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) ) / fse3t(:,:,jk) 
    380380               END DO 
    381381               etot3(:,:,nksr+1:jpk) = 0.e0                ! below 400m set to zero 
     
    399399                     zc0 = rn_abs * EXP( -fsdepw(ji,jj,jk  )*zsi0r ) + (1.-rn_abs) * EXP( -fsdepw(ji,jj,jk  )*zsi1r ) 
    400400                     zc1 = rn_abs * EXP( -fsdepw(ji,jj,jk+1)*zsi0r ) + (1.-rn_abs) * EXP( -fsdepw(ji,jj,jk+1)*zsi1r ) 
    401                      etot3(ji,jj,jk) = ro0cpr * (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1) ) / fse3t(ji,jj,jk) 
     401                     etot3(ji,jj,jk) = r1_rau0_rcp * ( zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1) ) / fse3t(ji,jj,jk) 
    402402                  END DO 
    403403               END DO 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/TRA/trasbc.F90

    r1739 r1858  
    9999      !!---------------------------------------------------------------------- 
    100100      USE oce, ONLY :   ztrdt => ua   ! use ua as 3D workspace    
    101       USE oce, ONLY :   ztrds => va   ! use va as 3D workspace    
     101      USE oce, ONLY :   ztrds => va   !  -  va     -      -    
    102102      !! 
    103103      INTEGER, INTENT(in) ::   kt     ! ocean time-step index 
     
    132132#endif 
    133133            IF( lk_vvl) THEN 
    134                zta = ro0cpr * qns(ji,jj) * zse3t &                   ! temperature : heat flux 
     134               zta = r1_rau0_rcp * qns(ji,jj) * zse3t &              ! temperature : heat flux 
    135135                &    - emp(ji,jj) * zsrau * tn(ji,jj,1)  * zse3t     ! & cooling/heating effet of EMP flux 
    136136               zsa = 0.e0                                            ! No salinity concent./dilut. effect 
    137137            ELSE 
    138                zta = ro0cpr * qns(ji,jj) * zse3t     ! temperature : heat flux 
     138               zta = r1_rau0_rcp * qns(ji,jj) * zse3t                ! temperature : heat flux 
    139139               zsa = emps(ji,jj) * zsrau * sn(ji,jj,1)   * zse3t     ! salinity :  concent./dilut. effect 
    140140            ENDIF 
  • branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/OPA_SRC/ZDF/zdfkpp.F90

    r1695 r1858  
    455455            zBo   (ji,jj) = grav * zthermal * qns(ji,jj) -  grav * zhalin * emps(ji,jj) 
    456456            ! Surface Temperature flux for non-local term 
    457             wt0(ji,jj) = - ( qsr(ji,jj) + qns(ji,jj) )* ro0cpr * tmask(ji,jj,1) 
     457            wt0(ji,jj) = - ( qsr(ji,jj) + qns(ji,jj) )* r1_rau0_rcp * tmask(ji,jj,1) 
    458458            ! Surface salinity flux for non-local term 
    459459            ws0(ji,jj) = - ( emps(ji,jj) * sn(ji,jj,1) * rcs ) * tmask(ji,jj,1) 
Note: See TracChangeset for help on using the changeset viewer.