Ignore:
Timestamp:
2014-02-21T15:20:39+01:00 (7 years ago)
Author:
vancop
Message:

[Heat conservation in LIM3, part HC1 (LIM_SRC_3_HC17)]

Location:
branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r4332 r4506  
    266266   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fbif        !: Heat flux at the ice base 
    267267   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdm_snw     !: Variation of snow mass over 1 time step     [Kg/m2] 
    268    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdq_snw     !: Heat content associated with rdm_snw        [J/m2] 
     268   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdq_snw     !: Heat per grid cell area released by snow melt    [J/m2], >0 if to the ocean 
    269269   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdm_ice     !: Variation of ice mass over 1 time step      [Kg/m2] 
    270    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdq_ice     !: Heat content associated with rdm_ice        [J/m2] 
     270   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdq_ice     !: Heat per grid cell area taken or released by ice growth and melt [J/m2], >0 if to the ocean 
    271271   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qldif       !: heat balance of the lead (or of the open ocean) 
    272272   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qcmif       !: Energy needed to bring the ocean to freezing  
     
    284284   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_res     !: residual salt flux due to correction of ice thickness [PSU/m2/s] 
    285285   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fhbri       !: heat flux due to brine rejection 
    286    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fheat_mec   !: heat flux associated with porous ridged ice formation [???] 
    287286   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fheat_res   !: residual heat flux due to correction of ice thickness 
    288287   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fmmec       !: mass flux due to snow loss during compression         [Kg/m2/s] 
     
    455454         &      fdtcn    (jpi,jpj) , qdtcn  (jpi,jpj) , fstric (jpi,jpj) , fscmbq   (jpi,jpj) ,     & 
    456455         &      ffltbif  (jpi,jpj) , fsbbq  (jpi,jpj) , qfvbq  (jpi,jpj) , dmgwi    (jpi,jpj) ,     & 
    457          &      sfx_res  (jpi,jpj) , sfx_bri(jpi,jpj) , sfx_mec(jpi,jpj) , fheat_mec(jpi,jpj) ,     & 
    458          &      fhbri    (jpi,jpj) , fmmec  (jpi,jpj) , sfx_thd(jpi,jpj) , fhmec    (jpi,jpj) ,     & 
     456         &      sfx_res  (jpi,jpj) , sfx_bri(jpi,jpj) , sfx_mec(jpi,jpj) , fhbri    (jpi,jpj) ,     & 
     457         &      fmmec  (jpi,jpj) , sfx_thd(jpi,jpj) , fhmec    (jpi,jpj) ,                          & 
    459458         &      fheat_res(jpi,jpj)                                                            , STAT=ierr(ii) ) 
    460459 
  • branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r4345 r4506  
    363363            !-----------------------------------------------------------------------------! 
    364364            fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice     ! fresh water source for ocean 
     365            ! there is a unit problem for e_snow_mlt that we should fix 
    365366            fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice     ! heat sink for ocean 
    366367 
     
    908909      INTEGER ::   ij                ! horizontal index, combines i and j loops 
    909910      INTEGER ::   icells            ! number of cells with aicen > puny 
    910       REAL(wp) ::   zindb, zsrdg2   ! local scalar 
     911      REAL(wp) ::   zindb, zsrdg2    ! local scalar 
    911912      REAL(wp) ::   hL, hR, farea, zdummy, zdummy0, ztmelts    ! left and right limits of integration 
     913      REAL(wp) ::   zsstK            ! SST in Kelvin 
    912914 
    913915      INTEGER , POINTER, DIMENSION(:) ::   indxi, indxj   ! compressed indices 
     
    10981100            esrdg(ji,jj) = esnon_init(ji,jj,jl1) * afrac(ji,jj) 
    10991101            srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por ) 
    1100             srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 
     1102            srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) !! MV HC 2014 this line seems useless 
    11011103 
    11021104            ! rafting volumes, heat contents ... 
     
    11201122            ! Salinity 
    11211123            !------------- 
    1122             smsw(ji,jj)  = sss_m(ji,jj) * vsw(ji,jj) * rhoic / rau0       ! salt content of seawater frozen in voids 
    1123  
     1124            smsw(ji,jj)  = vsw(ji,jj) * sss_m(ji,jj)                      ! salt content of seawater frozen in voids !! MV HC2014 
    11241125            zsrdg2       = srdg1(ji,jj) + smsw(ji,jj)                     ! salt content of new ridge 
    11251126 
     
    11281129            !                                                             ! excess of salt is flushed into the ocean 
    11291130            !sfx_mec(ji,jj) = sfx_mec(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice 
     1131            ! MV HC 2014 this previous line seems ok too... Maybe units are not that good 
    11301132 
    11311133            !rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic    ! gurvan: increase in ice volume du to seawater frozen in voids              
     1134            ! MV HC 2014 this previous line seems ok, i'm not sure at this moment of the sign convention 
    11321135 
    11331136            !------------------------------------             
     
    11581161               &                                + rhosn*vsrft(ji,jj)*(1.0-fsnowrft) 
    11591162 
     1163            ! MV HC 2014 - I think the units are wrong here... shit there is a 10-9 in the snow... 
     1164            ! we should therefore do something 
    11601165            esnow_mlt(ji,jj) = esnow_mlt(ji,jj) + esrdg(ji,jj)*(1.0-fsnowrdg)         &   !rafting included 
    11611166               &                                + esrft(ji,jj)*(1.0-fsnowrft)           
     
    11871192               eirft(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj) 
    11881193               e_i  (ji,jj,jk,jl1)  = e_i(ji,jj,jk,jl1) - erdg1(ji,jj,jk) - eirft(ji,jj,jk) 
    1189                ! sea water heat content 
    1190                ztmelts          = - tmut * sss_m(ji,jj) + rtt 
    1191                ! heat content per unit volume 
    1192                zdummy0          = - rcp * ( sst_m(ji,jj) + rt0 - rtt ) * vsw(ji,jj) 
    1193  
    1194                ! corrected sea water salinity 
    1195                zindb  = MAX( 0._wp , SIGN( 1._wp , vsw(ji,jj) - epsi20 ) ) 
    1196                zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / MAX( ridge_por * vsw(ji,jj), epsi20 ) 
    1197  
    1198                ztmelts          = - tmut * zdummy + rtt 
    1199                ersw(ji,jj,jk)   = - rcp * ( ztmelts - rtt ) * vsw(ji,jj) 
    1200  
    1201                ! heat flux 
    1202                fheat_mec(ji,jj) = fheat_mec(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) * r1_rdtice 
     1194                
     1195                
     1196               ! enthalpy of the trapped seawater (J/m2, >0) 
     1197               zsstK  = sst_m(ji,jj) + rt0 
     1198               ersw(ji,jj,jk)   = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 ) 
     1199 
     1200               ! heat flux to the ocean 
     1201               rdq_ice(ji,jj) = rdq_ice(ji,jj) + ersw(ji,jj,jk) !! MV HC 2014 
    12031202 
    12041203               ! Correct dimensions to avoid big values 
     
    12061205 
    12071206               ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 
     1207               ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean  
     1208               !! MV HC 2014 
    12081209               ersw (ji,jj,jk)  = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / REAL( nlay_i ) 
    12091210 
    12101211               erdg2(ji,jj,jk)  = erdg1(ji,jj,jk) + ersw(ji,jj,jk) 
     1212 
    12111213            END DO ! ij 
    12121214         END DO !jk 
  • branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r4345 r4506  
    106106      INTEGER  ::   iflt, ial , iadv , ifral, ifrdv   !   -      - 
    107107      REAL(wp) ::   zinda, zemp, zemp_snow, zfmm      ! local scalars 
    108       REAL(wp) ::   zemp_snw                          !   -      - 
     108      REAL(wp) ::   zemp_snw, zqmass,   zcd           !   -      -   
     109      REAL(wp) ::   zswitch                           !   -      -   
    109110      REAL(wp) ::   zfcm1 , zfcm2                     !   -      - 
    110111      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp     ! 2D/3D workspace 
     
    113114       
    114115      IF( lk_cpl )   CALL wrk_alloc( jpi, jpj, jpl, zalb, zalbp ) 
     116 
     117      SELECT CASE( nn_ice_embd )                 ! levitating or embedded sea-ice option                               
     118        CASE( 0    )   ;   zswitch = 1           ! (0) standard levitating sea-ice : salt exchange only               
     119        CASE( 1, 2 )   ;   zswitch = 0           ! (1) levitating sea-ice: salt and volume exchange but no pressure effect 
     120                                                 ! (2) embedded sea-ice : salt and volume fluxes and pressure             
     121      END SELECT                                 !                                                                       
    115122 
    116123      !------------------------------------------! 
     
    172179               &    + ifrdv   * (       qfvbq(ji,jj) +             qdtcn(ji,jj) ) * r1_rdtice   & 
    173180               &    + fhmec(ji,jj)                                                              & ! snow melt when ridging 
    174                &    + fheat_mec(ji,jj)                                                          & ! ridge formation 
    175181               &    + fheat_res(ji,jj)                                                            ! residual heat flux 
    176182            ! qcmif   Energy needed to bring the ocean surface layer until its freezing (ok) 
     
    184190            fsbbq(ji,jj) = ( 1._wp - ( ifvt + iflt ) ) * fscmbq(ji,jj)      
    185191 
    186             ! used to compute the oceanic heat flux at the next time step 
     192            ! heat flux associated with ice-ocean mass exchange 
     193            zqmass = ( rdq_snw(ji,jj)                                         &  
     194                 & + rdq_ice(ji,jj) * ( 1.- zswitch) ) * r1_rdtice   !  heat flux due to snow & ice heat content  
    187195            qsr(ji,jj) = zfcm1                                       ! solar heat flux  
    188             qns(ji,jj) = zfcm2 - fdtcn(ji,jj)                        ! non solar heat flux 
    189             !                           ! fdtcn : turbulent oceanic heat flux 
     196            qns(ji,jj) = zfcm2 - fdtcn(ji,jj) + zqmass               ! non solar heat flux    
     197                               ! fdtcn : turbulent oceanic heat flux 
    190198         END DO 
    191199      END DO 
  • branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r4332 r4506  
    184184               &                            + fdtcn(ji,jj)                        &   ! turbulent ice-ocean heat 
    185185               &                            + fsbbq(ji,jj) * ( 1.0 - zinda )  )   &   ! residual heat from previous step 
     186               ! MV HC 2014 check that this is good. 
     187               ! We should remove the heat content of precip that has fallen on sea ice 
    186188               & - pfrld(ji,jj)**betas * sprecip(ji,jj) * lfus                    )   ! latent heat of sprecip melting 
     189               ! MV HC 2014 partie heat content manque 
    187190            ! 
    188191            ! Positive heat budget is used for bottom ablation 
     
    282285            CALL tab_2d_1d( nbpb, qldif_1d   (1:nbpb), qldif           , jpi, jpj, npb(1:nbpb) ) 
    283286            CALL tab_2d_1d( nbpb, rdm_ice_1d (1:nbpb), rdm_ice         , jpi, jpj, npb(1:nbpb) ) 
     287            CALL tab_2d_1d( nbpb, rdq_ice_1d (1:nbpb), rdq_ice         , jpi, jpj, npb(1:nbpb) ) 
    284288            CALL tab_2d_1d( nbpb, rdm_snw_1d (1:nbpb), rdm_snw         , jpi, jpj, npb(1:nbpb) ) 
     289            CALL tab_2d_1d( nbpb, rdq_snw_1d (1:nbpb), rdq_snw         , jpi, jpj, npb(1:nbpb) ) 
    285290            CALL tab_2d_1d( nbpb, dmgwi_1d   (1:nbpb), dmgwi           , jpi, jpj, npb(1:nbpb) ) 
    286291            CALL tab_2d_1d( nbpb, qlbbq_1d   (1:nbpb), zqlbsbq         , jpi, jpj, npb(1:nbpb) ) 
     
    349354               CALL tab_1d_2d( nbpb, qfvbq         , npb, qfvbq_1d  (1:nbpb)   , jpi, jpj ) 
    350355               CALL tab_1d_2d( nbpb, rdm_ice       , npb, rdm_ice_1d(1:nbpb)   , jpi, jpj ) 
     356               CALL tab_1d_2d( nbpb, rdq_ice       , npb, rdq_ice_1d(1:nbpb)   , jpi, jpj ) 
    351357               CALL tab_1d_2d( nbpb, rdm_snw       , npb, rdm_snw_1d(1:nbpb)   , jpi, jpj ) 
     358               CALL tab_1d_2d( nbpb, rdq_snw       , npb, rdq_snw_1d(1:nbpb)   , jpi, jpj ) 
    352359               CALL tab_1d_2d( nbpb, dmgwi         , npb, dmgwi_1d  (1:nbpb)   , jpi, jpj ) 
    353360               CALL tab_1d_2d( nbpb, rdvosif       , npb, dvsbq_1d  (1:nbpb)   , jpi, jpj ) 
  • branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r4332 r4506  
    7474      INTEGER  ::   ji , jk        ! dummy loop indices 
    7575      INTEGER  ::   ii, ij         ! 2D corresponding indices to ji 
     76      INTEGER  ::   i_use          ! debugging flag 
    7677      INTEGER  ::   isnow          ! switch for presence (1) or absence (0) of snow 
    7778      INTEGER  ::   isnowic        ! snow ice formation not 
     
    8586      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment 
    8687      REAL(wp) ::   zcoeff       ! dummy argument for snowfall partitioning over ice and leads 
    87       REAL(wp) ::   zsm_snowice  ! snow-ice salinity 
     88      REAL(wp) ::   zs_snic  ! snow-ice salinity 
    8889      REAL(wp) ::   zswi1        ! switch for computation of bottom salinity 
    8990      REAL(wp) ::   zswi12       ! switch for computation of bottom salinity 
    9091      REAL(wp) ::   zswi2        ! switch for computation of bottom salinity 
    9192      REAL(wp) ::   zgrr         ! bottom growth rate 
    92       REAL(wp) ::   ztform       ! bottom formation temperature 
    93       ! 
     93      REAL(wp) ::   zt_i_new     ! bottom formation temperature 
     94 
     95      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean 
     96      REAL(wp) ::   zEi          ! specific enthalpy of sea ice (J/kg) 
     97      REAL(wp) ::   zEw          ! specific enthalpy of exchanged water (J/kg) 
     98      REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg) 
     99      REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean 
     100      REAL(wp) ::   zsstK        ! SST in Kelvin 
     101 
    94102      REAL(wp), POINTER, DIMENSION(:) ::   zh_i        ! ice layer thickness 
    95103      REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness 
     
    119127 
    120128      ! mass and salt flux (clem) 
    121       REAL(wp) :: zdvres, zdvsur, zdvbot 
     129      REAL(wp) :: zdvres, zdvsur, zdvbot, zswitch_sal 
    122130      REAL(wp), POINTER, DIMENSION(:) ::   zviold, zvsold   ! old ice volume... 
    123131 
     
    129137      REAL(wp), POINTER, DIMENSION(:) ::   zfbase, zdq_i 
    130138      !!------------------------------------------------------------------ 
     139 
     140      ! Discriminate between varying salinity (num_sal=2) and prescribed cases (other values) 
     141      SELECT CASE( num_sal )                       ! varying salinity or not 
     142         CASE( 1, 3, 4 ) ;   zswitch_sal = 0       ! prescribed salinity profile 
     143         CASE( 2 )       ;   zswitch_sal = 1       ! varying salinity profile 
     144      END SELECT 
    131145 
    132146      CALL wrk_alloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i ) 
     
    222236      DO ji = kideb, kiut 
    223237         ! tatm_ice is now in K 
    224          zqprec   (ji)   =  rhosn * ( cpic * ( rtt - tatm_ice_1d(ji) ) + lfus )   
     238         zqprec   (ji)   =  rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )   
    225239         zqfont_su(ji)   =  z_f_surf(ji) * rdt_ice 
    226240         zdeltah  (ji,1) =  MIN( 0.e0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) ) 
     
    272286      DO jk = 1, nlay_i 
    273287         DO ji = kideb, kiut  
    274             !                                                    ! melt of layer jk 
    275             zdeltah  (ji,jk) = - zqfont_su(ji) / q_i_b(ji,jk) 
    276             !                                                    ! recompute heat available 
    277             zqfont_su(ji   ) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk)  
    278             !                                                    ! melt of layer jk cannot be higher than its thickness 
    279             zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) ) 
    280             !                                                    ! update surface melt 
    281             dh_i_surf(ji   ) = dh_i_surf(ji) + zdeltah(ji,jk)  
    282             !                                                    ! for energy conservation 
    283             zdq_i    (ji   ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 
    284             ! 
    285             ! clem 
    286             sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji)    & 
     288            zEi            = - q_i_b(ji,jk) / rhoic                ! Specific enthalpy of layer k [J/kg, <0] 
     289 
     290            ztmelts        = - tmut * s_i_b(ji,jk) + rtt           ! Melting point of layer k [K] 
     291 
     292            zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0] 
     293 
     294            zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0 
     295 
     296            zfmdt          = - zqfont_su(ji) / zdE                 ! Mass flux to the ocean [kg/m2, >0] 
     297 
     298            zdeltah(ji,jk) = - zfmdt / rhoic                       ! Melt of layer jk [m, <0] 
     299 
     300            zqfont_su(ji)  = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * rhoic * ( - zdE ) 
     301                                                                   ! Energy remaining in case of melting of the full layer [J/m2, >0] 
     302            zdeltah(ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0] 
     303 
     304            dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt 
     305 
     306            zfmdt          = - rhoic*zdeltah(ji,jk)                ! Recompute mass flux [kg/m2, >0] 
     307 
     308            zQm            = MAX ( zfmdt, 0._wp ) * zEw            ! Energy of the melt water sent to the ocean [J/m2, <0] 
     309 
     310            ! Contribution to salt flux  
     311            sfx_thd_1d(ji)   = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji)    & 
    287312               &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice 
     313            ! Contribution to heat flux 
     314            rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm * a_i_b(ji) 
     315             
     316            ! Conservation test 
     317            zdq_i(ji)      = zdq_i(ji) + zdeltah(ji,jk) * rhoic * zdE * r1_rdtice  ! ? still don't know 
    288318         END DO 
    289319      END DO 
     
    364394      !------------------------------------------------------------------------------! 
    365395      ! 
    366       ! Ice basal growth / melt is given by the ratio of heat budget over basal 
    367       ! ice heat content.  Basal heat budget is given by the difference between 
    368       ! the inner conductive flux  (fc_bo_i), from the open water heat flux  
     396      !------------------ 
     397      ! 4.1 Basal growth  
     398      !------------------ 
     399      ! Basal growth is driven by heat imbalance at the ice-ocean interface, 
     400      ! between the inner conductive flux  (fc_bo_i), from the open water heat flux  
    369401      ! (qlbbqb) and the turbulent ocean flux (fbif).  
    370402      ! fc_bo_i is positive downwards. fbif and qlbbq are positive to the ice  
    371403 
    372       !----------------------------------------------------- 
    373       ! 4.1 Basal growth - (a) salinity not varying in time  
    374       !----------------------------------------------------- 
    375       IF(  num_sal /= 2  ) THEN   ! ice salinity constant in time 
     404      ! If salinity varies in time, an iterative procedure is required, because 
     405      ! the involved quantities are inter-dependent. 
     406      ! Basal growth (dh_i_bott) depends upon new ice specific enthalpy (zEi), 
     407      ! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bott) 
     408      ! -> need for an iterative procedure, which converges quickly 
     409 
     410      IF ( num_sal == 2 ) THEN 
     411         num_iter_max = 5 
     412      ELSE 
     413         num_iter_max = 1 
     414      ENDIF 
     415 
     416      ! Initialize dh_i_bott 
     417      dh_i_bott(:) = 0.0e0 
     418 
     419      ! Iterative procedure 
     420      DO iter = 1, num_iter_max 
    376421         DO ji = kideb, kiut 
    377             IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) < 0._wp  ) THEN 
    378                s_i_new(ji)         =  sm_i_b(ji) 
    379                ! Melting point in K 
    380                ztmelts             =  - tmut * s_i_new(ji) + rtt  
    381                ! New ice heat content (Bitz and Lipscomb, 1999) 
    382                ztform              =  t_i_b(ji,nlay_i)  ! t_bo_b crashes in the 
    383                ! Baltic 
    384                q_i_b(ji,nlay_i+1)  = rhoic * (  cpic * ( ztmelts - ztform )                                & 
    385                   &                           + lfus * (  1.0 - ( ztmelts - rtt ) / ( ztform - rtt )  )    & 
    386                   &                           - rcp  * ( ztmelts - rtt )                                 ) 
    387                ! Basal growth rate = - F*dt / q 
    388                dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)  
    389                sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice 
    390             ENDIF 
    391          END DO 
    392       ENDIF 
    393  
    394       !------------------------------------------------- 
    395       ! 4.1 Basal growth - (b) salinity varying in time  
    396       !------------------------------------------------- 
    397       IF(  num_sal == 2  ) THEN 
    398          ! the growth rate (dh_i_bott) is function of the new ice heat content (q_i_b(nlay_i+1)).  
    399          ! q_i_b depends on the new ice salinity (snewice).  
    400          ! snewice depends on dh_i_bott ; it converges quickly, so, no problem 
    401          ! See Vancoppenolle et al., OM08 for more info on this 
    402  
    403          ! Initial value (tested 1D, can be anything between 1 and 20) 
    404          num_iter_max = 4 
    405          s_i_new(:)   = 4.0 
    406  
    407          ! Iterative procedure 
    408          DO iter = 1, num_iter_max 
    409             DO ji = kideb, kiut 
    410                IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0.e0  ) THEN 
    411                   ii = MOD( npb(ji) - 1, jpi ) + 1 
    412                   ij = ( npb(ji) - 1 ) / jpi + 1 
    413                   ! Melting point in K 
    414                   ztmelts             =   - tmut * s_i_new(ji) + rtt  
    415                   ! New ice heat content (Bitz and Lipscomb, 1999) 
    416                   q_i_b(ji,nlay_i+1)  =  rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                             & 
    417                      &                            + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )   & 
    418                      &                            - rcp * ( ztmelts-rtt )                                     ) 
    419                   ! Bottom growth rate = - F*dt / q 
    420                   dh_i_bott(ji) =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
    421                   ! New ice salinity ( Cox and Weeks, JGR, 1988 ) 
    422                   ! zswi2  (1) if dh_i_bott/rdt .GT. 3.6e-7 
    423                   ! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8 
    424                   ! zswi1  (1) if dh_i_bott/rdt .LT. 2.0e-8 
    425                   zgrr   = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi13 ) ) 
    426                   zswi2  = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) )  
    427                   zswi12 = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) 
    428                   zswi1  = 1. - zswi2 * zswi12  
    429                   zfracs = zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   & 
    430                      &                   + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  
    431                   zfracs = MIN( 0.5 , zfracs ) 
    432                   s_i_new(ji) = zfracs * sss_m(ii,ij) 
    433                ENDIF ! fc_bo_i 
    434             END DO ! ji 
    435          END DO ! iter 
    436  
    437          ! Final values 
    438          DO ji = kideb, kiut 
    439             IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN 
    440                ! New ice salinity must not exceed 20 psu 
    441                s_i_new(ji) = MIN( s_i_new(ji), s_i_max ) 
    442                ! Metling point in K 
    443                ztmelts     =   - tmut * s_i_new(ji) + rtt  
    444                ! New ice heat content (Bitz and Lipscomb, 1999) 
    445                q_i_b(ji,nlay_i+1)  =  rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                              & 
    446                   &                            + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )    & 
    447                   &                            - rcp * ( ztmelts - rtt )                                    ) 
    448                ! Basal growth rate = - F*dt / q 
    449                dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
    450                ! Salinity update 
    451                ! entrapment during bottom growth 
    452                sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice 
    453             ENDIF ! heat budget 
    454          END DO 
    455       ENDIF 
     422 
     423            IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0.e0  ) THEN 
     424 
     425               ! New bottom ice salinity (Cox & Weeks, JGR88 ) 
     426               !--- zswi1  if dh/dt < 2.0e-8 
     427               !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7  
     428               !--- zswi2  if dh/dt > 3.6e-7 
     429               zgrr               = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi13 ) ) 
     430               zswi2              = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) ) 
     431               zswi12             = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) 
     432               zswi1              = 1. - zswi2 * zswi12 
     433               zfracs             = MIN ( zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   & 
     434                  &               + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 ) 
     435 
     436               ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
     437 
     438               s_i_new(ji)        = zswitch_sal * zfracs * sss_m(ii,ij)  &  ! New ice salinity 
     439                                  + ( 1. - zswitch_sal ) * sm_i_b(ji)  
     440               ! New ice growth 
     441               ztmelts            = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K) 
     442 
     443               zt_i_new           = zswitch_sal * t_bo_b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i) 
     444                
     445               zEi                = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)       
     446                  &               - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) )   & 
     447                  &               + rcp  * ( ztmelts-rtt )           
     448 
     449               zEw                = rcp  * ( t_bo_b(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0) 
     450 
     451               zdE                = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0) 
     452 
     453               dh_i_bott(ji)      = rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / ( zdE * rhoic ) 
     454 
     455               !!! not sure we still need the next line... useful to keep this in memory ? check limthd_ent... 
     456               q_i_b(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0) 
     457 
     458            ENDIF ! fc_bo_i 
     459         END DO ! ji 
     460      END DO ! iter 
     461 
     462      ! Contribution to Energy and Salt Fluxes 
     463      DO ji = kideb, kiut 
     464 
     465         zEw            = rcp  * ( t_bo_b(ji) - rt0 )                   ! Specific enthalpy of seawater (J/kg, < 0) 
     466 
     467         zfmdt          = - rhoic * dh_i_bott(ji)                       ! Mass flux x time step (kg/m2, < 0) 
     468 
     469         ! Contribution to energy flux to the ocean (J/m2) 
     470         rdq_ice_1d(ji) = rdq_ice_1d(ji) + zEw         * a_i_b(ji) * zfmdt 
     471 
     472         ! Contribution to salt flux  () 
     473         sfx_thd_1d(ji) = sfx_thd_1d(ji) + s_i_new(ji) * a_i_b(ji) * zfmdt * r1_rdtice 
     474 
     475      END DO 
    456476 
    457477      !---------------- 
     
    465485         ! heat convergence at the surface > 0 
    466486         IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0._wp  ) THEN 
    467             s_i_new(ji)   =  s_i_b(ji,nlay_i) 
     487            s_i_new(ji) = s_i_b(ji,nlay_i) ! is this line still useful ? 
    468488            zqfont_bo(ji) =  rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) 
    469489            zfbase(ji)    =  zqfont_bo(ji) * r1_rdtice     ! heat conservation test 
     
    476496         DO ji = kideb, kiut 
    477497            IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji)  >=  0._wp  ) THEN 
    478                ztmelts = - tmut * s_i_b(ji,jk) + rtt  
    479                IF( t_i_b(ji,jk) >= ztmelts ) THEN   !!gm : a comment is needed 
    480                   zdeltah   (ji,jk) = - zh_i(ji) 
    481                   dh_i_bott (ji   ) = dh_i_bott(ji) + zdeltah(ji,jk) 
    482                   zinnermelt(ji   ) = 1._wp 
    483                ELSE                                  ! normal ablation 
    484                   zdeltah  (ji,jk) = - zqfont_bo(ji) / q_i_b(ji,jk) 
    485                   zqfont_bo(ji   ) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk) 
    486                   zdeltah  (ji,jk) = MAX(zdeltah(ji,jk), - zh_i(ji) ) 
    487                   dh_i_bott(ji   ) = dh_i_bott(ji) + zdeltah(ji,jk) 
    488                   zdq_i    (ji   ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 
     498 
     499               ztmelts = - tmut * s_i_b(ji,jk) + rtt  ! Melting point of layer jk (K) 
     500 
     501               IF( t_i_b(ji,jk) >= ztmelts ) THEN !!! Internal melting 
     502 
     503                  zdeltah   (ji,jk) = - zh_i(ji)  ! internal melting occurs when the internal temperature is above freezing      
     504                                                  ! this should normally not happen, but sometimes, heat diffusion leads to this 
     505 
     506                  dh_i_bott (ji)    = dh_i_bott(ji) + zdeltah(ji,jk) 
     507 
     508                  zinnermelt(ji)    = 1._wp 
     509 
     510                  zQm               = 0.          ! Not sure which specific enthalpy we should use here (MV HC 2014) 
     511                                                  ! If that happens, heat is probably not well counted 
     512                                                  ! Put zero by default, but more bug analysis should be done to investigate this case 
     513 
     514               ELSE                               !!! Basal melting 
     515 
     516                  zEi               = - q_i_b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0) 
     517 
     518                  zEw               = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0) 
     519 
     520                  zdE               = zEi - zEw              ! Specific enthalpy difference   (J/kg, <0) 
     521 
     522                  zfmdt             = - zqfont_bo(ji) / zdE  ! Mass flux x time step (kg/m2, >0) 
     523 
     524                  zdeltah(ji,jk)    = - zfmdt / rhoic        ! Gross thickness change 
     525 
     526                  zqfont_bo(ji)     = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * rhoic * ( - zdE ) ! Update heat available 
     527 
     528                  zdeltah(ji,jk)    = MAX( zdeltah(ji,jk), - zh_i(ji) ) ! Update thickness change 
     529 
     530                  dh_i_bott(ji)     = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt 
     531 
     532                  zfmdt             = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step 
     533 
     534                  zQm               = MAX ( zfmdt , 0.0 ) * zEw         ! Heat exchanged with ocean 
     535 
     536                  zdq_i(ji)         = zdq_i(ji) + zdeltah(ji,jk) * rhoic * zdE * r1_rdtice ! for heat conservation  
     537 
    489538               ENDIF 
    490                ! clem: contribution to salt flux 
     539 
     540               ! Contribution to salt flux 
    491541               sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji)    & 
    492542                    &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice 
     543 
     544               ! Contribution to energy flux to the ocean (J/m2) 
     545               rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm * a_i_b(ji) 
     546 
    493547            ENDIF 
    494548         END DO ! ji 
     
    571625         dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji) 
    572626 
    573         ! remaining heat 
     627         ! remaining heat 
    574628         zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) +  zqfont_bo(ji) )  
    575629 
     
    661715      DO ji = kideb, kiut 
    662716         ht_i_b(ji) = zhgnew(ji) 
    663       END DO  ! ji 
     717      END DO ! ji  
     718       
    664719      ! 
    665720      !------------------------------------------------------------------------------| 
     
    681736         !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew     - ht_s_b(ji) ) * rhosn  
    682737 
    683          !        Equivalent salt flux (1) Snow-ice formation component 
    684          !        ----------------------------------------------------- 
    685          ii = MOD( npb(ji) - 1, jpi ) + 1 
    686          ij =    ( npb(ji) - 1 ) / jpi + 1 
    687  
    688          IF( num_sal == 2 ) THEN   ;   zsm_snowice = sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic 
    689          ELSE                      ;   zsm_snowice = sm_i_b(ji)    
    690          ENDIF 
     738         ! Salinity of snow ice 
     739         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
     740 
     741         zs_snic = zswitch_sal          * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic     & 
     742                 + ( 1. - zswitch_sal ) * sm_i_b(ji) 
     743 
    691744         ! entrapment during snow ice formation 
    692745         ! clem: new salinity difference stored (to be used in limthd_ent.F90) 
     
    694747            i_ice_switch = MAX( 0._wp , SIGN( 1._wp , zhgnew(ji) - epsi10 ) ) 
    695748            ! salinity dif due to snow-ice formation 
    696             dsm_i_si_1d(ji) = ( zsm_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch      
     749            dsm_i_si_1d(ji) = ( zs_snic - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch      
    697750            ! salinity dif due to bottom growth  
    698751            IF (  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji)  < 0._wp ) THEN 
     
    709762 
    710763         ! diagnostic ( snow ice growth ) 
    711          ii = MOD( npb(ji) - 1, jpi ) + 1 
    712          ij =    ( npb(ji) - 1 ) / jpi + 1 
     764         ii = MOD( npb(ji) - 1, jpi ) + 1 !MV HC 2014 useless 
     765         ij =    ( npb(ji) - 1 ) / jpi + 1 ! MV HC 2014 useless 
    713766         diag_sni_gr(ii,ij)  = diag_sni_gr(ii,ij) + dh_snowice(ji)*a_i_b(ji) * r1_rdtice 
    714          ! 
    715          ! salt flux 
    716          sfx_thd_1d(ji) = sfx_thd_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice 
    717          !-------------------------------- 
    718          ! Update mass fluxes (clem) 
    719          !-------------------------------- 
     767 
     768         ! Contribution to energy flux to the ocean [J/m2], >0 
     769         zfmdt          = ( rhosn - rhoic ) * MAX( dh_snowice(ji), 0.0 )    ! <0 
     770         zsstK          = sst_m(ii,ij) + rt0                                 
     771         zEw            = rcp * ( zsstK - rt0 ) 
     772         zQm            = zfmdt * zEw  
     773         rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm * a_i_b(ji) 
     774 
     775         ! Contribution to salt flux 
     776         sfx_thd_1d(ji) = sfx_thd_1d(ji) + sss_m(ii,ij) * a_i_b(ji) * zfmdt * r1_rdtice  
     777           
     778         ! Contribution to mass fluxes 
    720779         rdm_ice_1d(ji) = rdm_ice_1d(ji) + ( a_i_b(ji) * ht_i_b(ji) - zviold(ji) ) * rhoic  
    721780         rdm_snw_1d(ji) = rdm_snw_1d(ji) + ( a_i_b(ji) * ht_s_b(ji) - zvsold(ji) ) * rhosn  
  • branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90

    r4332 r4506  
    2222   USE domain         ! 
    2323   USE phycst         ! physical constants 
     24   USE sbc_oce        ! Surface boundary condition: ocean fields 
    2425   USE ice            ! LIM variables 
    2526   USE par_ice        ! LIM parameters 
     
    8687         ztmelts        ,   &  ! ice melting point 
    8788         zqsnic         ,   &  ! enthalpy of snow ice layer 
     89         zsstK          ,   &  ! SST in Kelvin 
    8890         zhsnow         ,   &  ! temporary snow thickness variable 
    8991         zswitch        ,   &  ! dummy switch argument 
     
    540542      DO ji = kideb, kiut 
    541543         ! energy of the flooding seawater 
    542          zqsnic = rau0 * rcp * ( rtt - t_bo_b(ji) ) * dh_snowice(ji) * & 
    543             (rhoic - rhosn) / rhoic * REAL(snicswi(ji)) ! generally positive 
     544         ii = MOD( npb(ji) - 1, jpi ) + 1 
     545         ij =    ( npb(ji) - 1 ) / jpi + 1 
     546         zsstK  = sst_m(ii,ij) + rt0 
     547         zqsnic = ( rhosn - rhoic ) * dh_snowice(ji) * rcp * ( zsstK - rt0 ) * REAL ( snicswi(ji) ) ! MV HC 2014 
     548          
    544549         ! Heat conservation diagnostic 
    545550         qt_i_in(ji,jl) = qt_i_in(ji,jl) + zqsnic  
    546  
    547          qldif_1d(ji)   = qldif_1d(ji) + zqsnic * a_i_b(ji) 
    548551 
    549552         ! enthalpy of the newly formed snow-ice layer 
  • branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

    r4332 r4506  
    8181      LOGICAL  ::   iterate_frazil   ! iterate frazil ice collection thickness 
    8282      CHARACTER (len = 15) :: fieldid 
    83       ! 
     83 
     84      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean) 
     85      REAL(wp) ::   zEi          ! sea ice specific enthalpy (J/kg) 
     86      REAL(wp) ::   zEw          ! seawater specific enthalpy (J/kg) 
     87      REAL(wp) ::   zfmdt        ! mass flux x time step (kg/m2, >0 towards ocean) 
     88        
    8489      INTEGER , POINTER, DIMENSION(:) ::   zcatac      ! indexes of categories where new ice grows 
    8590      REAL(wp), POINTER, DIMENSION(:) ::   zswinew     ! switch for new ice or not 
     
    323328         CALL tab_2d_1d( nbpac, sfx_thd_1d(1:nbpac)     , sfx_thd, jpi, jpj, npac(1:nbpac) ) 
    324329         CALL tab_2d_1d( nbpac, rdm_ice_1d(1:nbpac)     , rdm_ice, jpi, jpj, npac(1:nbpac) ) 
     330         CALL tab_2d_1d( nbpac, rdq_ice_1d(1:nbpac)     , rdq_ice, jpi, jpj, npac(1:nbpac) ) 
    325331         CALL tab_2d_1d( nbpac, hicol_b   (1:nbpac)     , hicol  , jpi, jpj, npac(1:nbpac) ) 
    326332         CALL tab_2d_1d( nbpac, zvrel_ac  (1:nbpac)     , zvrel  , jpi, jpj, npac(1:nbpac) ) 
     
    365371               &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )   & 
    366372               &                       - rcp  *         ( ztmelts - rtt )  ) 
     373            ! MV HC 2014 comment I dont see why this line below is here... ? 
     374            ! This implies that ze_newice gets to rhoic*Lfus if it was negative, but this should never happen 
    367375            ze_newice(ji) =   MAX( ze_newice(ji) , 0._wp )    & 
    368376               &          +   MAX(  0.0 , SIGN( 1.0 , - ze_newice(ji) )  ) * rhoic * lfus 
     
    386394         !------------------- 
    387395         DO ji = 1, nbpac 
    388             zv_newice(ji) = - zqbgow(ji) / ze_newice(ji) 
     396 
     397            zEi           = - ze_newice(ji) / rhoic                ! specific enthalpy of forming ice [J/kg] 
     398 
     399            zEw           = rcp * ( t_bo_b(ji) - rt0 )             ! specific enthalpy of seawater at t_bo_b [J/kg] 
     400 
     401            zdE           = zEi - zEw                              ! specific enthalpy difference [J/kg] 
     402 
     403            zfmdt         = - zqbgow(ji) / zdE                     ! Fm.dt [kg/m2] (<0)  
     404 
     405            zv_newice(ji) = - zfmdt / rhoic 
     406 
     407            zQm           = zfmdt * zEw                            ! heat to the ocean >0 associated with mass flux   
     408 
     409            ! Contribution to energy flux to the ocean [J/m2], >0   
     410            rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm  
    389411 
    390412            ! A fraction zfrazb of frazil ice is accreted at the ice bottom 
     
    539561         END DO 
    540562!!gm ???  why the previous do loop  if ocerwriten by the following one ? 
     563!! MV HC 2014 just because it'snot the same index and the expression is different 
    541564         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 
    542565            DO ji = 1, nbpac 
     
    623646         CALL tab_1d_2d( nbpac, sfx_thd, npac(1:nbpac), sfx_thd_1d(1:nbpac), jpi, jpj ) 
    624647         CALL tab_1d_2d( nbpac, rdm_ice, npac(1:nbpac), rdm_ice_1d(1:nbpac), jpi, jpj ) 
     648         CALL tab_1d_2d( nbpac, rdq_ice, npac(1:nbpac), rdq_ice_1d(1:nbpac), jpi, jpj ) 
    625649         ! 
    626650      ENDIF ! nbpac > 0 
  • branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90

    r4332 r4506  
    448448         CALL prt_ctl(tab2d_1=fmmec  , clinfo1= ' lim_update1 : fmmec : ', tab2d_2=fhmec     , clinfo2= ' fhmec     : ') 
    449449         CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' lim_update1 : sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ') 
    450          CALL prt_ctl(tab2d_1=fhbri  , clinfo1= ' lim_update1 : fhbri : ', tab2d_2=fheat_mec , clinfo2= ' fheat_mec : ') 
     450         CALL prt_ctl(tab2d_1=fhbri  , clinfo1= ' lim_update1 : fhbri : ')  
    451451 
    452452         CALL prt_ctl_info(' ') 
  • branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90

    r4332 r4506  
    598598         CALL prt_ctl(tab2d_1=fmmec  , clinfo1= ' lim_update2 : fmmec : ', tab2d_2=fhmec     , clinfo2= ' fhmec     : ') 
    599599         CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' lim_update2 : sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ') 
    600          CALL prt_ctl(tab2d_1=fhbri  , clinfo1= ' lim_update2 : fhbri : ', tab2d_2=fheat_mec , clinfo2= ' fheat_mec : ') 
     600         CALL prt_ctl(tab2d_1=fhbri  , clinfo1= ' lim_update2 : fhbri : ') 
    601601 
    602602         CALL prt_ctl_info(' ') 
  • branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90

    r4045 r4506  
    6666   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fbif_1d       !: <==> the 2D  fbif 
    6767   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdm_ice_1d    !: <==> the 2D  rdm_ice 
     68   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdq_ice_1d    !: <==> the 2D  rdq_ice 
    6869   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdm_snw_1d    !: <==> the 2D  rdm_snw 
     70   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdq_snw_1d    !: <==> the 2D  rdq_snw 
    6971   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qlbbq_1d      !: <==> the 2D  qlbsbq 
    7072   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dmgwi_1d      !: <==> the 2D  dmgwi 
     
    164166      ALLOCATE( sprecip_1d (jpij) , frld_1d    (jpij) , at_i_b     (jpij) ,     & 
    165167         &      fbif_1d    (jpij) , rdm_ice_1d (jpij) , rdm_snw_1d (jpij) ,     & 
     168         &                          rdq_ice_1d (jpij) , rdq_snw_1d (jpij) ,     & !MV HC 2014 
    166169         &      qlbbq_1d   (jpij) , dmgwi_1d   (jpij) , dvsbq_1d   (jpij) ,     & 
    167170         &      dvbbq_1d   (jpij) , dvlbq_1d   (jpij) , dvnbq_1d   (jpij) ,     & 
  • branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90

    r3625 r4506  
    5454   REAL(wp), PUBLIC ::   r1_rau0                     !: = 1. / rau0                   [m3/kg] 
    5555   REAL(wp), PUBLIC ::   rauw     = 1000._wp         !: volumic mass of pure water    [m3/kg] 
    56    REAL(wp), PUBLIC ::   rcp      =    4.e3_wp       !: ocean specific heat           [J/Kelvin] 
    57    REAL(wp), PUBLIC ::   r1_rcp                      !: = 1. / rcp                    [Kelvin/J] 
     56   REAL(wp), PUBLIC ::   rcp      =    4.e3_wp       !: ocean specific heat           [J/kg/K] 
     57   REAL(wp), PUBLIC ::   r1_rcp                      !: = 1. / rcp                    [kg.K/J] 
    5858   REAL(wp), PUBLIC ::   r1_rau0_rcp                 !: = 1. / ( rau0 * rcp ) 
    5959 
     
    6969#if defined key_lim3 || defined key_cice 
    7070   REAL(wp), PUBLIC ::   rhoic    =  917._wp         !: volumic mass of sea ice                               [kg/m3] 
    71    REAL(wp), PUBLIC ::   rcdic    =    2.034396_wp   !: thermal conductivity of fresh ice 
    72    REAL(wp), PUBLIC ::   rcdsn    =    0.31_wp       !: thermal conductivity of snow 
    73    REAL(wp), PUBLIC ::   cpic     = 2067.0_wp        !: specific heat for ice  
     71   REAL(wp), PUBLIC ::   rcdic    =    2.034396_wp   !: thermal conductivity of fresh ice                     [W/m/K] 
     72   REAL(wp), PUBLIC ::   rcdsn    =    0.31_wp       !: thermal conductivity of snow                          [W/m/K]  
     73   REAL(wp), PUBLIC ::   cpic     = 2067.0_wp        !: specific heat for ice                                 [J/kg/K] 
    7474   REAL(wp), PUBLIC ::   lsub     =    2.834e+6_wp   !: pure ice latent heat of sublimation                   [J/kg] 
    7575   REAL(wp), PUBLIC ::   lfus     =    0.334e+6_wp   !: latent heat of fusion of fresh ice                    [J/kg] 
    76    REAL(wp), PUBLIC ::   tmut     =    0.054_wp      !: decrease of seawater meltpoint with salinity 
     76   REAL(wp), PUBLIC ::   tmut     =    0.054_wp      !: decrease of seawater meltpoint with salinity          [degC/ppt] 
    7777   REAL(wp), PUBLIC ::   xlsn                        !: = lfus*rhosn (volumetric latent heat fusion of snow)  [J/m3] 
    7878#else 
  • branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r4358 r4506  
    136136      REAL(wp) ::   zcoef   ! local scalar 
    137137      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice_os, zalb_ice_cs  ! albedo of the ice under overcast/clear sky 
    138       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice      ! mean albedo of ice (for coupled) 
     138      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice      ! mean albedo of ice  
    139139 
    140140      REAL(wp), POINTER, DIMENSION(:,:) :: zalb_ice_all    ! Mean albedo over all categories 
     
    152152      IF( nn_timing == 1 )  CALL timing_start('sbc_ice_lim') 
    153153 
    154       CALL wrk_alloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs ) 
     154      CALL wrk_alloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs, zalb_ice ) 
    155155 
    156156#if defined key_coupled 
    157       IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) CALL wrk_alloc( jpi,jpj,jpl, zalb_ice) 
    158157      IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) & 
    159158         &   CALL wrk_alloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 
     
    168167         ! 
    169168         IF( ln_nicep ) THEN      ! control print at a given point 
    170             jiindx = 177   ;   jjindx = 112 
     169            jiindx = 14    ;   jjindx =  25 
    171170            IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx 
    172171         ENDIF 
     
    310309         sfx    (:,:) = 0._wp   ;   sfx_thd  (:,:) = 0._wp 
    311310         sfx_bri(:,:) = 0._wp   ;   sfx_mec  (:,:) = 0._wp   ;   sfx_res  (:,:) = 0._wp 
    312          fhbri  (:,:) = 0._wp   ;   fheat_mec(:,:) = 0._wp   ;   fheat_res(:,:) = 0._wp 
     311         fhbri  (:,:) = 0._wp   ;    
     312         fheat_res(:,:) = 0._wp  
    313313         fhmec  (:,:) = 0._wp   ;    
    314314         fmmec  (:,:) = 0._wp 
     
    327327         rdm_snw(:,:) = 0._wp   ! variation of snow mass per unit area 
    328328         rdm_ice(:,:) = 0._wp   ! variation of ice mass per unit area 
     329         rdq_snw(:,:) = 0._wp   ! heat flux associated with mass exchange (ice contribution) 
     330         rdq_ice(:,:) = 0._wp   ! heat flux associated with mass exchange (snow contribution) 
    329331         hicifp (:,:) = 0._wp   ! daily thermodynamic ice production.  
    330332         ! 
     
    421423       
    422424!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!! 
    423       ! 
    424       CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs ) 
     425      CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs, zalb_ice ) 
    425426 
    426427#if defined key_coupled 
    427       IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice) 
     428   
    428429      IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) & 
    429430         &    CALL wrk_dealloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 
     
    807808               WRITE(numout,*) ' fhmec      : ', fhmec    (ji,jj) 
    808809               WRITE(numout,*) ' fhbri      : ', fhbri    (ji,jj) 
    809                WRITE(numout,*) ' fheat_mec  : ', fheat_mec(ji,jj) 
    810810               WRITE(numout,*)  
    811811               WRITE(numout,*) ' sst        : ', sst_m(ji,jj)   
Note: See TracChangeset for help on using the changeset viewer.