New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 4634 for branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 – NEMO

Ignore:
Timestamp:
2014-05-12T22:46:18+02:00 (10 years ago)
Author:
clem
Message:

major changes in heat budget

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r4345 r4634  
    55   !!====================================================================== 
    66   !! History :  LIM  ! 2006-02  (M. Vancoppenolle) Original code  
    7    !!            3.2  ! 2009-07  (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_mec 
     7   !!            3.2  ! 2009-07  (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_dyn 
    88   !!            4.0  ! 2011-02  (G. Madec) dynamical allocation 
    99   !!---------------------------------------------------------------------- 
     
    143143      REAL(wp), POINTER, DIMENSION(:,:) ::   esnow_mlt       ! energy needed to melt snow in ocean (J m-2) 
    144144      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final  !  ice volume summed over categories 
    145       REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 
     145      REAL(wp) :: zchk_v_i, zchk_smv, zchk_e_i, zchk_fs, zchk_fw, zchk_ft, zchk_v_i_b, zchk_smv_b, zchk_e_i_b, zchk_fs_b, zchk_fw_b, zchk_ft_b  
    146146      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
    147       ! mass and salt flux (clem) 
    148       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold, zvsold, zsmvold   ! old ice volume... 
    149147      !!----------------------------------------------------------------------------- 
    150148      IF( nn_timing == 1 )  CALL timing_start('limitd_me') 
    151149 
    152150      CALL wrk_alloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 
    153  
    154       CALL wrk_alloc( jpi, jpj, jpl, zviold, zvsold, zsmvold )   ! clem 
    155151 
    156152      IF( numit == nstart  )   CALL lim_itd_me_init   ! Initialization (first time-step only) 
     
    165161      !- check conservation (C Rousset) 
    166162      IF (ln_limdiahsb) THEN 
    167          zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
     163         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) 
    168164         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
    169          zchk_fw_b  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 
    170          zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 
     165         zchk_e_i_b = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) 
     166         zchk_fw_b  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) 
     167         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) 
     168         zchk_ft_b  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) 
    171169      ENDIF 
    172170      !- check conservation (C Rousset) 
    173171      ! ------------------------------- 
    174  
    175       ! mass and salt flux init (clem) 
    176       zviold(:,:,:) = v_i(:,:,:) 
    177       zvsold(:,:,:) = v_s(:,:,:) 
    178       zsmvold(:,:,:) = smv_i(:,:,:) 
    179172 
    180173      !-----------------------------------------------------------------------------! 
     
    362355            ! 5) Heat, salt and freshwater fluxes 
    363356            !-----------------------------------------------------------------------------! 
    364             fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice     ! fresh water source for ocean 
    365             fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice     ! heat sink for ocean 
     357            wfx_snw(ji,jj) = wfx_snw(ji,jj) - msnow_mlt(ji,jj) * r1_rdtice     ! fresh water source for ocean 
     358            hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * unit_fac / area(ji,jj) * r1_rdtice  ! heat sink for ocean (<0, W.m-2) 
    366359 
    367360         END DO 
     
    399392      CALL lim_itd_me_zapsmall 
    400393 
    401       !-------------------------------- 
    402       ! Update mass/salt fluxes (clem) 
    403       !-------------------------------- 
    404       DO jl = 1, jpl 
    405          DO jj = 1, jpj  
    406             DO ji = 1, jpi 
    407                diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice 
    408                rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * rhoic  
    409                rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * rhosn  
    410                sfx_mec(ji,jj) = sfx_mec(ji,jj) - ( smv_i(ji,jj,jl) - zsmvold(ji,jj,jl) ) * rhoic * r1_rdtice  
    411             END DO 
    412          END DO 
    413       END DO 
    414394 
    415395      IF(ln_ctl) THEN     ! Control print 
     
    448428      !- check conservation (C Rousset) 
    449429      IF (ln_limdiahsb) THEN 
    450          zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 
    451          zchk_fw  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 
     430         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 
     431         zchk_fw  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b 
     432         zchk_ft  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b 
    452433  
    453          zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 
     434         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw  
    454435         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 
     436         zchk_e_i =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zchk_e_i_b * r1_rdtice + zchk_ft 
    455437 
    456438         zchk_vmin = glob_min(v_i) 
     
    459441        
    460442         IF(lwp) THEN 
    461             IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limitd_me) = ',(zchk_v_i * rday) 
     443            IF ( ABS( zchk_v_i   ) >  1.e-2 ) WRITE(numout,*) 'violation volume [kg/day]     (limitd_me) = ',(zchk_v_i * rday) 
    462444            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_me) = ',(zchk_smv * rday) 
     445            IF ( ABS( zchk_e_i   ) >  1.e-4 ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (limitd_me) = ',(zchk_e_i) 
    463446            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_me) = ',(zchk_vmin * 1.e-3) 
    464447            IF ( zchk_amax >  kamax+epsi10  ) WRITE(numout,*) 'violation a_i>amax            (limitd_me) = ',zchk_amax 
     
    472455      ! 
    473456      CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 
    474       ! 
    475       CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zsmvold )   ! clem 
    476457      ! 
    477458      IF( nn_timing == 1 )  CALL timing_stop('limitd_me') 
     
    908889      INTEGER ::   ij                ! horizontal index, combines i and j loops 
    909890      INTEGER ::   icells            ! number of cells with aicen > puny 
    910       REAL(wp) ::   zindb, zsrdg2   ! local scalar 
     891      REAL(wp) ::   zindb    ! local scalar 
    911892      REAL(wp) ::   hL, hR, farea, zdummy, zdummy0, ztmelts    ! left and right limits of integration 
     893      REAL(wp) ::   zsstK            ! SST in Kelvin 
    912894 
    913895      INTEGER , POINTER, DIMENSION(:) ::   indxi, indxj   ! compressed indices 
     
    917899 
    918900      REAL(wp), POINTER, DIMENSION(:,:,:) ::   aicen_init, vicen_init   ! ice  area    & volume before ridging 
    919       REAL(wp), POINTER, DIMENSION(:,:,:) ::   vsnon_init, esnon_init   ! snow volume  & energy before ridging 
     901      REAL(wp), POINTER, DIMENSION(:,:,:) ::   vsnwn_init, esnwn_init   ! snow volume  & energy before ridging 
    920902      REAL(wp), POINTER, DIMENSION(:,:,:) ::   smv_i_init, oa_i_init    ! ice salinity & age    before ridging 
    921903 
     
    952934      CALL wrk_alloc( jpi, jpj,             vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw ) 
    953935      CALL wrk_alloc( jpi, jpj,             afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 
    954       CALL wrk_alloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnon_init, esnon_init, smv_i_init, oa_i_init ) 
     936      CALL wrk_alloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 
    955937      CALL wrk_alloc( jpi, jpj, jkmax,      eirft, erdg1, erdg2, ersw ) 
    956938      CALL wrk_alloc( jpi, jpj, jkmax, jpl, eicen_init ) 
     
    1008990         aicen_init(:,:,jl) = a_i(:,:,jl) 
    1009991         vicen_init(:,:,jl) = v_i(:,:,jl) 
    1010          vsnon_init(:,:,jl) = v_s(:,:,jl) 
     992         vsnwn_init(:,:,jl) = v_s(:,:,jl) 
    1011993         ! 
    1012994         smv_i_init(:,:,jl) = smv_i(:,:,jl) 
     
    1014996      END DO !jl 
    1015997 
    1016       esnon_init(:,:,:) = e_s(:,:,1,:) 
     998      esnwn_init(:,:,:) = e_s(:,:,1,:) 
    1017999 
    10181000      DO jl = 1, jpl   
     
    10951077            vsw  (ji,jj) = vrdg1(ji,jj) * ridge_por 
    10961078 
    1097             vsrdg(ji,jj) = vsnon_init(ji,jj,jl1) * afrac(ji,jj) 
    1098             esrdg(ji,jj) = esnon_init(ji,jj,jl1) * afrac(ji,jj) 
     1079            vsrdg(ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj) 
     1080            esrdg(ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj) 
    10991081            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) 
     1082            srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) !! MV HC 2014 this line seems useless 
    11011083 
    11021084            ! rafting volumes, heat contents ... 
    11031085            virft(ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj) 
    1104             vsrft(ji,jj) = vsnon_init(ji,jj,jl1) * afrft(ji,jj) 
    1105             esrft(ji,jj) = esnon_init(ji,jj,jl1) * afrft(ji,jj) 
     1086            vsrft(ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj) 
     1087            esrft(ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj) 
    11061088            smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj)  
    11071089 
     
    11201102            ! Salinity 
    11211103            !------------- 
    1122             smsw(ji,jj)  = sss_m(ji,jj) * vsw(ji,jj) * rhoic / rau0       ! salt content of seawater frozen in voids 
    1123  
    1124             zsrdg2       = srdg1(ji,jj) + smsw(ji,jj)                     ! salt content of new ridge 
    1125  
    1126             srdg2(ji,jj) = MIN( s_i_max * vrdg2(ji,jj) , zsrdg2 )         ! impose a maximum salinity 
     1104            smsw(ji,jj)  = vsw(ji,jj) * sss_m(ji,jj)                      ! salt content of seawater frozen in voids !! MV HC2014 
     1105            srdg2(ji,jj) = srdg1(ji,jj) + smsw(ji,jj)                     ! salt content of new ridge 
     1106 
     1107            !srdg2(ji,jj) = MIN( s_i_max * vrdg2(ji,jj) , zsrdg2 )         ! impose a maximum salinity 
    11271108             
    1128             !                                                             ! excess of salt is flushed into the ocean 
    1129             !sfx_mec(ji,jj) = sfx_mec(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice 
    1130  
    1131             !rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic    ! gurvan: increase in ice volume du to seawater frozen in voids              
     1109            sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ji,jj) * rhoic * r1_rdtice 
     1110            wfx_dyn(ji,jj) = wfx_dyn(ji,jj) + vsw (ji,jj) * rhoic * r1_rdtice   ! gurvan: increase in ice volume du to seawater frozen in voids              
     1111            ! MV HC 2014 this previous line seems ok, i'm not sure at this moment of the sign convention 
    11321112 
    11331113            !------------------------------------             
     
    11581138               &                                + rhosn*vsrft(ji,jj)*(1.0-fsnowrft) 
    11591139 
    1160             esnow_mlt(ji,jj) = esnow_mlt(ji,jj) + esrdg(ji,jj)*(1.0-fsnowrdg)         &   !rafting included 
    1161                &                                + esrft(ji,jj)*(1.0-fsnowrft)           
     1140            ! in 1e-9 Joules (same as e_s) 
     1141            esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0-fsnowrdg)         &   !rafting included 
     1142               &                                - esrft(ji,jj)*(1.0-fsnowrft)           
    11621143 
    11631144            !----------------------------------------------------------------- 
     
    11871168               eirft(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj) 
    11881169               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 
     1170                
     1171                
     1172               ! enthalpy of the trapped seawater (J/m2, >0) 
     1173               ! clem: if sst>0, then ersw <0 (is that possible?) 
     1174               zsstK  = sst_m(ji,jj) + rt0 
     1175               ersw(ji,jj,jk)   = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 ) / REAL( nlay_i ) 
     1176 
     1177               ! heat flux to the ocean 
     1178               hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ji,jj,jk) * r1_rdtice  ! > 0 [W.m-2] ocean->ice flux  
    12031179 
    12041180               ! Correct dimensions to avoid big values 
    1205                ersw(ji,jj,jk)   = ersw(ji,jj,jk) * 1.e-09 
    1206  
    1207                ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 
    1208                ersw (ji,jj,jk)  = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / REAL( nlay_i ) 
     1181               ersw(ji,jj,jk)   = ersw(ji,jj,jk) / unit_fac 
     1182 
     1183               ! Mutliply by ice volume, and divide by number of layers to get heat content in 1.e9 J 
     1184               ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean  
     1185               !! MV HC 2014 
     1186               ersw (ji,jj,jk)  = ersw(ji,jj,jk) * area(ji,jj) 
    12091187 
    12101188               erdg2(ji,jj,jk)  = erdg1(ji,jj,jk) + ersw(ji,jj,jk) 
     1189 
    12111190            END DO ! ij 
    12121191         END DO !jk 
     
    13611340      CALL wrk_dealloc( jpi, jpj,             vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw ) 
    13621341      CALL wrk_dealloc( jpi, jpj,             afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 
    1363       CALL wrk_dealloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnon_init, esnon_init, smv_i_init, oa_i_init ) 
     1342      CALL wrk_dealloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 
    13641343      CALL wrk_dealloc( jpi, jpj, jkmax,      eirft, erdg1, erdg2, ersw ) 
    13651344      CALL wrk_dealloc( jpi, jpj, jkmax, jpl, eicen_init ) 
     
    14551434 
    14561435      REAL(wp), POINTER, DIMENSION(:,:) ::   zmask   ! 2D workspace 
    1457       REAL(wp)                          ::   zmask_glo 
     1436      REAL(wp)                          ::   zmask_glo, zsal, zvi, zvs, zei, zes 
    14581437!!gm      REAL(wp) ::   xtmp      ! temporary variable 
    14591438      !!------------------------------------------------------------------- 
     
    14711450         DO jj = 1, jpj 
    14721451            DO ji = 1, jpi 
    1473                IF(  ( a_i(ji,jj,jl) >= -epsi10 .AND. a_i(ji,jj,jl) <  0._wp   ) .OR.   & 
    1474                   & ( a_i(ji,jj,jl) >  0._wp   .AND. a_i(ji,jj,jl) <= epsi10  ) .OR.   & 
    1475                   & ( v_i(ji,jj,jl) == 0._wp   .AND. a_i(ji,jj,jl) >  0._wp   ) .OR.   & 
    1476                   & ( v_i(ji,jj,jl) >  0._wp   .AND. v_i(ji,jj,jl) <= epsi10  ) )   zmask(ji,jj) = 1._wp 
     1452!               IF(  ( a_i(ji,jj,jl) >= -epsi10 .AND. a_i(ji,jj,jl) <  0._wp   ) .OR.   & 
     1453!                  & ( a_i(ji,jj,jl) >  0._wp   .AND. a_i(ji,jj,jl) <= epsi10  ) .OR.   & 
     1454!                  & ( v_i(ji,jj,jl) == 0._wp   .AND. a_i(ji,jj,jl) >  0._wp   ) .OR.   & 
     1455!                  & ( v_i(ji,jj,jl) >  0._wp   .AND. v_i(ji,jj,jl) <= epsi10  ) )   zmask(ji,jj) = 1._wp 
     1456               IF(  ( a_i(ji,jj,jl) >= -epsi10 .AND. a_i(ji,jj,jl) <= epsi10  ) .OR.   & 
     1457                  & ( v_i(ji,jj,jl) >= 0._wp   .AND. v_i(ji,jj,jl) <= epsi10  ) )   zmask(ji,jj) = 1._wp 
    14771458            END DO 
    14781459         END DO 
     
    14871468            DO jj = 1 , jpj 
    14881469               DO ji = 1 , jpi 
    1489 !!gm                  xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) * r1_rdtice 
    1490 !!gm                  xtmp = xtmp * unit_fac 
    1491                   ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp 
     1470                  zei  = e_i(ji,jj,jk,jl) 
    14921471                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) ) 
     1472                  ! update exchanges with ocean 
     1473                  hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
    14931474               END DO 
    14941475            END DO 
     
    14971478         DO jj = 1 , jpj 
    14981479            DO ji = 1 , jpi 
    1499  
     1480                
     1481               zsal = smv_i(ji,jj,jl) 
     1482               zvi  = v_i(ji,jj,jl) 
     1483               zvs  = v_s(ji,jj,jl) 
     1484               zes  = e_s(ji,jj,1,jl) 
    15001485               !----------------------------------------------------------------- 
    15011486               ! Zap snow energy and use ocean heat to melt snow 
     
    15071492               ! fluxes are positive to the ocean 
    15081493               ! here the flux has to be negative for the ocean 
    1509 !!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice 
    1510                !           fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp 
    1511  
    1512 !!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice !RB   ??????? 
    1513  
    15141494               t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) ) 
    15151495 
     
    15321512               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
    15331513               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
    1534                ! 
     1514               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) ) 
     1515               ! additional condition 
     1516               IF( v_s(ji,jj,jl) <= epsi10 ) THEN 
     1517                  v_s(ji,jj,jl)   = 0._wp 
     1518                  e_s(ji,jj,1,jl) = 0._wp 
     1519               ENDIF 
     1520               ! update exchanges with ocean 
     1521               sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 
     1522               wfx_res(ji,jj)  = wfx_res(ji,jj) + ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice 
     1523               wfx_snw(ji,jj)  = wfx_snw(ji,jj) + ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice 
     1524               hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
    15351525            END DO 
    15361526         END DO 
Note: See TracChangeset for help on using the changeset viewer.