Changeset 8562


Ignore:
Timestamp:
2017-09-25T21:11:19+02:00 (3 years ago)
Author:
clem
Message:

cosmetics only

Location:
branches/2017/dev_r8183_ICEMODEL/NEMOGCM
Files:
13 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/field_def_nemo-lim.xml

    r8306 r8562  
    129129         <field id="albedo"       long_name="Mean albedo over sea ice and ocean"                                                                                              unit=""             /> 
    130130 
    131        <field id="iceamp"       long_name="melt pond fraction"                                           standard_name="sea_ice_meltpond_fraction"                          unit="%"            />  
     131    <field id="iceamp"       long_name="melt pond fraction"                                           standard_name="sea_ice_meltpond_fraction"                          unit="%"            />  
    132132         <field id="icevmp"       long_name="melt pond volume"                                             standard_name="sea_ice_meltpond_volume"                            unit="m"            />  
    133133 
     
    154154         <field id="tau_icebfr"   long_name="ice friction on ocean bottom for landfast ice"                unit="N/2"      /> 
    155155 
    156          <field id="icetrp"       long_name="ice volume transport"                                         unit="m/day"          /> 
    157          <field id="snwtrp"       long_name="snw volume transport"                                         unit="m/day"          /> 
    158          <field id="saltrp"       long_name="salt content transport"                                       unit="1e-3*kg/m2/day" /> 
     156         <field id="icetrp"       long_name="ice mass transport"                                           unit="kg/m2/s"          /> 
     157         <field id="snwtrp"       long_name="snw mass transport"                                           unit="kg/m2/s"          /> 
     158         <field id="saltrp"       long_name="salt     transport"                                           unit="1e-3*kg/m2/s" /> 
    159159         <field id="deitrp"       long_name="advected ice enthalpy"                                        unit="W/m2"           /> 
    160160         <field id="destrp"       long_name="advected snw enthalpy"                                        unit="W/m2"           /> 
    161161 
    162          <field id="sfxbri"       long_name="brine salt flux"                                              unit="1e-3*kg/m2/day" /> 
    163          <field id="sfxdyn"       long_name="salt flux from ridging rafting"                               unit="1e-3*kg/m2/day" /> 
    164          <field id="sfxres"       long_name="salt flux from lipupdate (resultant)"                         unit="1e-3*kg/m2/day" /> 
    165          <field id="sfxbog"       long_name="salt flux from bot growth"                                    unit="1e-3*kg/m2/day" /> 
    166          <field id="sfxbom"       long_name="salt flux from bot melt"                                      unit="1e-3*kg/m2/day" /> 
    167          <field id="sfxsum"       long_name="salt flux from surf melt"                                     unit="1e-3*kg/m2/day" /> 
    168          <field id="sfxlam"       long_name="salt flux from lateral melt"                                  unit="1e-3*kg/m2/day" /> 
    169          <field id="sfxsni"       long_name="salt flux from snow-ice formation"                            unit="1e-3*kg/m2/day" /> 
    170          <field id="sfxopw"       long_name="salt flux from open water ice formation"                      unit="1e-3*kg/m2/day" /> 
    171          <field id="sfxsub"       long_name="salt flux from sublimation"                                   unit="1e-3*kg/m2/day" /> 
    172          <field id="sfx"          long_name="Salt flux from sea ice"                                       unit="1e-3*kg/m2/day" /> 
    173  
    174          <field id="vfxbog"       long_name="daily bottom thermo ice prod."                                unit="m/day"   /> 
    175          <field id="vfxdyn"       long_name="daily  dynamic ice prod."                                     unit="m/day"   /> 
    176          <field id="vfxopw"       long_name="daily lateral thermo ice prod."                               unit="m/day"   /> 
    177          <field id="vfxsni"       long_name="daily snowice ice prod."                                      unit="m/day"   /> 
    178          <field id="vfxsum"       long_name="surface melt"                                                 unit="m/day"   /> 
    179          <field id="vfxlam"       long_name="lateral melt"                                                 unit="m/day"   /> 
    180          <field id="vfxbom"       long_name="bottom melt"                                                  unit="m/day"   /> 
    181          <field id="vfxres"       long_name="daily resultant ice prod./melting from limupdate"             unit="m/day"   /> 
    182          <field id="vfxice"       long_name="ice melt/growth"                                              unit="m/day"   /> 
    183          <field id="vfxsnw"       long_name="snw melt/growth"                                              unit="m/day"   /> 
    184          <field id="vfxsub"       long_name="snw sublimation"                                              unit="m/day"   /> 
    185          <field id="vfxsub_err"   long_name="excess of snw sublimation sent to ocean"                      unit="m/day"   /> 
    186          <field id="vfxspr"       long_name="snw precipitation on ice"                                     unit="m/day"   /> 
    187          <field id="vfxthin"      long_name="daily thermo ice prod. for thin ice(20cm) + open water"       unit="m/day"   /> 
     162         <field id="sfxbri"       long_name="salt flux from brines"                                        unit="1e-3*kg/m2/s" /> 
     163         <field id="sfxdyn"       long_name="salt flux from ridging rafting"                               unit="1e-3*kg/m2/s" /> 
     164         <field id="sfxres"       long_name="salt flux from lipupdate (resultant)"                         unit="1e-3*kg/m2/s" /> 
     165         <field id="sfxbog"       long_name="salt flux from bot growth"                                    unit="1e-3*kg/m2/s" /> 
     166         <field id="sfxbom"       long_name="salt flux from bot melt"                                      unit="1e-3*kg/m2/s" /> 
     167         <field id="sfxsum"       long_name="salt flux from surf melt"                                     unit="1e-3*kg/m2/s" /> 
     168         <field id="sfxlam"       long_name="salt flux from lateral melt"                                  unit="1e-3*kg/m2/s" /> 
     169         <field id="sfxsni"       long_name="salt flux from snow-ice formation"                            unit="1e-3*kg/m2/s" /> 
     170         <field id="sfxopw"       long_name="salt flux from open water ice formation"                      unit="1e-3*kg/m2/s" /> 
     171         <field id="sfxsub"       long_name="salt flux from sublimation"                                   unit="1e-3*kg/m2/s" /> 
     172         <field id="sfx"          long_name="Salt flux from sea ice"                                       unit="1e-3*kg/m2/s" /> 
     173 
     174         <field id="vfxbog"       long_name="bottom thermo ice prod."                                      unit="kg/m2/s"   /> 
     175         <field id="vfxdyn"       long_name="dynamic ice prod."                                            unit="kg/m2/s"   /> 
     176         <field id="vfxopw"       long_name="lateral thermo ice prod."                                     unit="kg/m2/s"   /> 
     177         <field id="vfxsni"       long_name="snowice ice prod."                                            unit="kg/m2/s"   /> 
     178         <field id="vfxsum"       long_name="surface melt"                                                 unit="kg/m2/s"   /> 
     179         <field id="vfxlam"       long_name="lateral melt"                                                 unit="kg/m2/s"   /> 
     180         <field id="vfxbom"       long_name="bottom melt"                                                  unit="kg/m2/s"   /> 
     181         <field id="vfxres"       long_name="resultant ice prod./melting"                                  unit="kg/m2/s"   /> 
     182         <field id="vfxice"       long_name="ice melt/growth"                                              unit="kg/m2/s"   /> 
     183         <field id="vfxsnw"       long_name="snw melt/growth"                                              unit="kg/m2/s"   /> 
     184         <field id="vfxsub"       long_name="snw sublimation"                                              unit="kg/m2/s"   /> 
     185         <field id="vfxsub_err"   long_name="excess of snw sublimation sent to ocean"                      unit="kg/m2/s"   /> 
     186         <field id="vfxspr"       long_name="snw precipitation on ice"                                     unit="kg/m2/s"   /> 
     187         <field id="vfxthin"      long_name="thermo ice prod. for thin ice(20cm) + open water"             unit="kg/m2/s"   /> 
    188188 
    189189         <field id="afxtot"       long_name="area tendency (total)"                                        unit="s-1"   /> 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r8534 r8562  
    172172   !                                     !!** ice-dynamics namelist (namdyn) ** 
    173173   REAL(wp), PUBLIC ::   rn_ishlat        !: lateral boundary condition for sea-ice 
    174    REAL(wp), PUBLIC ::   rn_cio           !: drag coefficient for oceanic stress 
    175174   LOGICAL , PUBLIC ::   ln_landfast      !: landfast ice parameterization (T or F)  
    176175   REAL(wp), PUBLIC ::   rn_gamma         !:    fraction of ocean depth that ice must reach to initiate landfast ice 
     
    184183   REAL(wp), PUBLIC ::   rn_relast        !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp)  
    185184   ! 
    186    !                                     !!** ice-thermodynamics namelist (namthd) ** 
     185   !                                     !!** ice-surface forcing namelist (namforcing) ** 
    187186                                          ! -- icethd_dh -- ! 
    188187   REAL(wp), PUBLIC ::   rn_blow_s        !: coef. for partitioning of snowfall between leads and sea ice 
    189188                                          ! -- icethd -- ! 
     189   REAL(wp), PUBLIC ::   rn_cio           !: drag coefficient for oceanic stress 
    190190   INTEGER , PUBLIC ::   nn_iceflx        !: Redistribute heat flux over ice categories 
    191191   !                                      !   =-1  Do nothing (needs N(cat) fluxes) 
     
    235235   !                                     !!** define arrays 
    236236   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   u_oce,v_oce !: surface ocean velocity used in ice dynamics 
    237    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hicol       !: ice collection thickness accreted in leads 
     237   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_i_new    !: ice collection thickness accreted in leads 
    238238   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   strength    !: ice strength 
    239239   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   stress1_i, stress2_i, stress12_i   !: 1st, 2nd & diagonal stress tensor element 
     
    437437      ! stay within Fortran's max-line length limit. 
    438438      ii = 1 
    439       ALLOCATE( u_oce   (jpi,jpj) , v_oce    (jpi,jpj) , hicol    (jpi,jpj) ,                        & 
     439      ALLOCATE( u_oce   (jpi,jpj) , v_oce    (jpi,jpj) , ht_i_new (jpi,jpj) ,                        & 
    440440         &      strength(jpi,jpj) , stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) ,  & 
    441441         &      delta_i (jpi,jpj) , divu_i   (jpi,jpj) , shear_i  (jpi,jpj) , STAT=ierr(ii) ) 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90

    r8559 r8562  
    208208      ii = ii + 1 
    209209      ALLOCATE( t_s_1d  (jpij,nlay_s)     , t_i_1d (jpij,nlay_i)     , s_i_1d(jpij,nlay_i) ,  &             
    210          &      e_i_1d  (jpij,nlay_i+1)   , e_s_1d (jpij,nlay_s)     ,                        & 
     210         &      e_i_1d  (jpij,nlay_i   , e_s_1d (jpij,nlay_s)     ,                        & 
    211211         &      eh_i_old(jpij,0:nlay_i+1) , h_i_old(jpij,0:nlay_i+1) , STAT=ierr(ii) ) 
    212212      ! 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedia.F90

    r8534 r8562  
    7676      IF( kt == nit000 .AND. lwp ) THEN 
    7777         WRITE(numout,*) 
    78          WRITE(numout,*)'icedia : outpout ice diagnostics (integrated over the domain)' 
     78         WRITE(numout,*)'icedia: output ice diagnostics (integrated over the domain)' 
    7979         WRITE(numout,*)'~~~~~~' 
    8080      ENDIF 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90

    r8534 r8562  
    8787            WHERE( ht_i(:,:,jl) > ht_n(:,:) * rn_gamma )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 
    8888         END DO 
     89         IF( iom_use('tau_icebfr') )   CALL iom_put( 'tau_icebfr', tau_icebfr )   
    8990      ENDIF 
    9091 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv.F90

    r8534 r8562  
    105105      diag_trp_vi (:,:) = SUM(     v_i  (:,:,:)          - v_i_b  (:,:,:)                  , dim=3 ) * r1_rdtice 
    106106      diag_trp_vs (:,:) = SUM(     v_s  (:,:,:)          - v_s_b  (:,:,:)                  , dim=3 ) * r1_rdtice 
    107       IF( iom_use('icetrp') )   CALL iom_put( "icetrp" , diag_trp_vi * rday  )         ! ice volume transport 
    108       IF( iom_use('snwtrp') )   CALL iom_put( "snwtrp" , diag_trp_vs * rday  )         ! snw volume transport 
    109       IF( iom_use('saltrp') )   CALL iom_put( "saltrp" , diag_trp_smv * rday * rhoic ) ! salt content transport 
    110       IF( iom_use('deitrp') )   CALL iom_put( "deitrp" , diag_trp_ei         )         ! advected ice enthalpy (W/m2) 
    111       IF( iom_use('destrp') )   CALL iom_put( "destrp" , diag_trp_es         )         ! advected snw enthalpy (W/m2) 
     107      IF( iom_use('icetrp') )   CALL iom_put( "icetrp" , diag_trp_vi )          ! ice volume transport 
     108      IF( iom_use('snwtrp') )   CALL iom_put( "snwtrp" , diag_trp_vs )          ! snw volume transport 
     109      IF( iom_use('saltrp') )   CALL iom_put( "saltrp" , diag_trp_smv * rhoic ) ! salt content transport 
     110      IF( iom_use('deitrp') )   CALL iom_put( "deitrp" , diag_trp_ei )          ! advected ice enthalpy (W/m2) 
     111      IF( iom_use('destrp') )   CALL iom_put( "destrp" , diag_trp_es )          ! advected snw enthalpy (W/m2) 
    112112 
    113113      IF( lrst_ice ) THEN                       !* write Prather fields in the restart file 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90

    r8534 r8562  
    128128         !    utau_ice, vtau_ice = surface ice stress [N/m2] 
    129129         !------------------------------------------------! 
    130          CALL ice_forcing_tau( kt, ksbc, utau_ice, vtau_ice ) 
    131                                        
     130                                        CALL ice_forcing_tau( kt, ksbc, utau_ice, vtau_ice )           
    132131         !-------------------------------------! 
    133132         ! --- ice dynamics and advection  --- ! 
     
    136135         CALL ice_rst_opn( kt )     ! Open Ice restart file (if necessary)  
    137136         ! 
    138          IF( ln_icedyn .AND. .NOT.lk_c1d )   CALL ice_dyn( kt )            ! -- Ice dynamics 
     137         IF( ln_icedyn .AND. .NOT.lk_c1d )   & 
     138            &                           CALL ice_dyn( kt )            ! -- Ice dynamics 
    139139         ! 
    140140         !                          !==  lateral boundary conditions  ==! 
    141141#if defined key_agrif 
    142          IF( .NOT. Agrif_Root()     )        CALL agrif_interp_lim3('T')   ! -- AGRIF  
     142         IF( .NOT. Agrif_Root()     )   CALL agrif_interp_lim3('T')   ! -- AGRIF  
    143143#endif 
    144          IF( ln_icethd .AND. ln_bdy )        CALL bdy_ice( kt )            ! -- bdy ice thermo 
     144         IF( ln_icethd .AND. ln_bdy )   CALL bdy_ice( kt )            ! -- bdy ice thermo 
    145145         ! 
    146146         ! 
    147147         !                          !==  previous lead fraction and ice volume for flux calculations 
    148          ! 
    149148         CALL ice_var_glo2eqv            ! ht_i and ht_s for ice albedo calculation 
    150149         CALL ice_var_agg(1)             ! at_i for coupling  
     
    165164         !    qprec_ice, qevap_ice 
    166165         !------------------------------------------------------! 
    167                                     CALL ice_forcing_flx( kt, ksbc ) 
    168  
     166                                        CALL ice_forcing_flx( kt, ksbc ) 
    169167         !----------------------------! 
    170168         ! --- ice thermodynamics --- ! 
    171169         !----------------------------! 
    172          IF( ln_icethd )            CALL ice_thd( kt )          ! -- Ice thermodynamics       
    173  
    174          ! MV MP 2016 
    175          IF ( ln_pnd )              CALL lim_mp( kt )           ! -- Melt ponds 
    176          ! END MV MP 2016 
    177  
    178          IF( ln_icethd )            CALL ice_cor( kt , 2 )      ! -- Corrections 
    179          ! --- 
     170         IF( ln_icethd )                CALL ice_thd( kt )            ! -- Ice thermodynamics       
     171         ! 
     172         IF ( ln_pnd )                  CALL lim_mp( kt )             ! -- Melt ponds 
     173         ! 
     174         IF( ln_icethd )                CALL ice_cor( kt , 2 )        ! -- Corrections 
    180175# if defined key_agrif 
    181          IF( .NOT. Agrif_Root() )   CALL agrif_update_lim3( kt ) 
     176         IF( .NOT. Agrif_Root() )       CALL agrif_update_lim3( kt ) 
    182177# endif 
    183                                     CALL ice_var_glo2eqv        ! necessary calls (at least for coupling) 
    184                                     CALL ice_var_agg( 2 )       ! necessary calls (at least for coupling) 
    185                                     ! 
     178         CALL ice_var_glo2eqv        ! necessary calls (at least for coupling) 
     179         CALL ice_var_agg( 2 )       ! necessary calls (at least for coupling) 
     180                                     ! 
    186181!! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work 
    187182!!       moreover it should only be called at the update frequency (as in agrif_lim3_update.F90) 
    188183!!# if defined key_agrif 
    189 !!         IF( .NOT. Agrif_Root() )   CALL Agrif_ChildGrid_To_ParentGrid() 
     184!!         IF( .NOT. Agrif_Root() )     CALL Agrif_ChildGrid_To_ParentGrid() 
    190185!!# endif 
    191                                     CALL ice_update_flx( kt )   ! -- Update ocean surface mass, heat and salt fluxes 
     186                                        CALL ice_update_flx( kt )     ! -- Update ocean surface mass, heat and salt fluxes 
    192187!!# if defined key_agrif 
    193 !!         IF( .NOT. Agrif_Root() )   CALL Agrif_ParentGrid_To_ChildGrid() 
     188!!         IF( .NOT. Agrif_Root() )     CALL Agrif_ParentGrid_To_ChildGrid() 
    194189!!# endif 
    195          IF( ln_icediahsb )         CALL ice_dia( kt )          ! -- Diagnostics and outputs  
    196          ! 
    197                                     CALL ice_wri( 1 )           ! -- Ice outputs  
    198          ! 
    199          ! 
    200          IF( lrst_ice )             CALL ice_rst_write( kt )    ! -- Ice restart file  
    201          ! 
    202          IF( ln_icectl )            CALL ice_ctl( kt )          ! -- alerts in case of model crash 
     190         IF( ln_icediahsb )             CALL ice_dia( kt )            ! -- Diagnostics outputs  
     191         ! 
     192                                        CALL ice_wri( kt )            ! -- Ice outputs  
     193         ! 
     194         IF( lrst_ice )                 CALL ice_rst_write( kt )      ! -- Ice restart file  
     195         ! 
     196         IF( ln_icectl )                CALL ice_ctl( kt )            ! -- alerts in case of model crash 
    203197         ! 
    204198      ENDIF   ! End sea-ice time step only 
     
    207201      ! --- Ocean time step --- ! 
    208202      !-------------------------! 
    209       ! Update surface ocean stresses (only in ice-dynamic case) otherwise the atm.-ocean stresses are used everywhere 
    210       !    using before instantaneous surf. currents 
    211       IF( ln_icedyn )               CALL ice_update_tau( kt, ub(:,:,1), vb(:,:,1) ) 
     203      IF( ln_icedyn )                   CALL ice_update_tau( kt, ub(:,:,1), vb(:,:,1) )   ! -- update surface ocean stresses 
    212204!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!! 
    213205      ! 
    214       IF( nn_timing == 1 ) CALL timing_stop('ice_stp') 
     206      IF( nn_timing == 1 )   CALL timing_stop('ice_stp') 
    215207      ! 
    216208   END SUBROUTINE ice_stp 
     
    246238      CALL ice_itd_init                ! ice thickness distribution initialization 
    247239      ! 
    248       ! MV MP 2016 
    249240      CALL lim_mp_init                 ! set melt ponds parameters (clem: important to be located here) 
    250       ! END MV MP 2016 
     241      ! 
    251242      !                                ! Initial sea-ice state 
    252243      IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90

    r8559 r8562  
    116116      END SELECT 
    117117 
    118       DO ji = 1, nidx 
    119          icount (ji,:) = 0 
    120          zdh_s_mel(ji) = 0._wp 
    121          e_i_1d(ji,nlay_i+1) = 0._wp ! Initialize enthalpy at nlay_i+1 
    122       END DO 
    123  
    124118      ! initialize layer thicknesses and enthalpies 
    125119      h_i_old (1:nidx,0:nlay_i+1) = 0._wp 
     
    227221      ! If heat still available (zq_su > 0), then melt more snow 
    228222      zdeltah(1:nidx,:) = 0._wp 
     223      zdh_s_mel(1:nidx) = 0._wp 
    229224      DO jk = 1, nlay_s 
    230225         DO ji = 1, nidx 
     
    443438 
    444439               dh_i_bott(ji)      = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) ) 
    445  
    446                e_i_1d(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0) 
    447440                
    448441            END DO 
     
    475468 
    476469            ! update heat content (J.m-2) and layer thickness 
    477             eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * e_i_1d(ji,nlay_i+1) 
     470            eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * (-zEi * rhoic) 
    478471            h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji) 
    479472 
     
    654647 
    655648      ! --- ensure that a_i = 0 where ht_i = 0 --- 
    656       DO ji = 1, nidx 
    657          IF( ht_i_1d(ji) == 0._wp )   a_i_1d(ji) = 0._wp 
    658       END DO          
     649      WHERE( ht_i_1d(1:nidx) == 0._wp )   a_i_1d(1:nidx) = 0._wp 
    659650      ! 
    660651   END SUBROUTINE ice_thd_dh 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_do.F90

    r8559 r8562  
    123123      ! 3) Collection thickness of ice formed in leads and polynyas 
    124124      !------------------------------------------------------------------------------!     
    125       ! hicol is the thickness of new ice formed in open water 
    126       ! hicol can be either prescribed (frazswi = 0) or computed (frazswi = 1) 
     125      ! ht_i_new is the thickness of new ice formed in open water 
     126      ! ht_i_new can be either prescribed (frazswi = 0) or computed (frazswi = 1) 
    127127      ! Frazil ice forms in open water, is transported by wind 
    128128      ! accumulates at the edge of the consolidated ice edge 
     
    136136 
    137137      ! Default new ice thickness 
    138       WHERE( qlead(:,:) < 0._wp )   ;   hicol(:,:) = rn_hinew 
    139       ELSEWHERE                     ;   hicol(:,:) = 0._wp 
     138      WHERE( qlead(:,:) < 0._wp )   ;   ht_i_new(:,:) = rn_hinew 
     139      ELSEWHERE                     ;   ht_i_new(:,:) = 0._wp 
    140140      END WHERE 
    141141 
     
    145145         ! Physical constants 
    146146         !-------------------- 
    147          hicol(:,:) = 0._wp 
     147         ht_i_new(:,:) = 0._wp 
    148148 
    149149         zhicrit = 0.04 ! frazil ice thickness 
     
    192192                  ! Iterative procedure 
    193193                  !--------------------- 
    194                   hicol(ji,jj) = zhicrit +   ( zhicrit + 0.1 )    & 
     194                  ht_i_new(ji,jj) = zhicrit +   ( zhicrit + 0.1 )    & 
    195195                     &                   / ( ( zhicrit + 0.1 ) * ( zhicrit + 0.1 ) -  zhicrit * zhicrit ) * ztwogp * zvrel2 
    196196 
    197197                  iter = 1 
    198198                  DO WHILE ( iter < 20 )  
    199                      zf  = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj) * hicol(ji,jj) - zhicrit * zhicrit ) -   & 
    200                         &    hicol(ji,jj) * zhicrit * ztwogp * zvrel2 
    201                      zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0 * hicol(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2 
    202  
    203                      hicol(ji,jj) = hicol(ji,jj) - zf / MAX( zfp, epsi20 ) 
     199                     zf  = ( ht_i_new(ji,jj) - zhicrit ) * ( ht_i_new(ji,jj) * ht_i_new(ji,jj) - zhicrit * zhicrit ) -   & 
     200                        &    ht_i_new(ji,jj) * zhicrit * ztwogp * zvrel2 
     201                     zfp = ( ht_i_new(ji,jj) - zhicrit ) * ( 3.0 * ht_i_new(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2 
     202 
     203                     ht_i_new(ji,jj) = ht_i_new(ji,jj) - zf / MAX( zfp, epsi20 ) 
    204204                     iter = iter + 1 
    205205                  END DO 
     
    210210         END DO  
    211211         !  
    212          CALL lbc_lnk_multi( zvrel, 'T', 1., hicol, 'T', 1.  ) 
     212         CALL lbc_lnk_multi( zvrel, 'T', 1., ht_i_new, 'T', 1.  ) 
    213213 
    214214      ENDIF ! End of computation of frazil ice collection thickness 
     
    252252         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx) , sfx_opw     ) 
    253253         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx) , wfx_opw     ) 
    254          CALL tab_2d_1d( nidx, idxice(1:nidx), zh_newice (1:nidx) , hicol       ) 
     254         CALL tab_2d_1d( nidx, idxice(1:nidx), zh_newice (1:nidx) , ht_i_new    ) 
    255255         CALL tab_2d_1d( nidx, idxice(1:nidx), zvrel_1d  (1:nidx) , zvrel       ) 
    256256 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_ent.F90

    r8534 r8562  
    7878      !  1) Cumulative integral of old enthalpy * thickness and layers interfaces 
    7979      !-------------------------------------------------------------------------- 
    80       zeh_cum0(:,0:nlay_i+2) = 0._wp  
    81       zh_cum0 (:,0:nlay_i+2) = 0._wp 
     80      zeh_cum0(1:nidx,0) = 0._wp  
     81      zh_cum0 (1:nidx,0) = 0._wp 
    8282      DO jk0 = 1, nlay_i+2 
    8383         DO ji = 1, nidx 
     
    9696 
    9797      ! new layers interfaces 
    98       zh_cum1(:,0:nlay_i) = 0._wp 
     98      zh_cum1(1:nidx,0) = 0._wp 
    9999      DO jk1 = 1, nlay_i 
    100100         DO ji = 1, nidx 
     
    103103      END DO 
    104104 
    105       zeh_cum1(:,0:nlay_i) = 0._wp  
     105      zeh_cum1(1:nidx,0) = 0._wp  
    106106      ! new cumulative q*h => linear interpolation 
    107107      DO jk0 = 1, nlay_i+1 
     
    117117      END DO 
    118118      ! to ensure that total heat content is strictly conserved, set: 
    119       zeh_cum1(:,nlay_i) = zeh_cum0(:,nlay_i+2)  
     119      zeh_cum1(1:nidx,nlay_i) = zeh_cum0(1:nidx,nlay_i+2)  
    120120 
    121121      ! new enthalpies 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90

    r8534 r8562  
    8484      !!------------------------------------------------------------------- 
    8585      INTEGER ::   ji, jk         ! spatial loop index 
    86       INTEGER ::   numeq          ! current reference number of equation 
    87       INTEGER ::   minnumeqmin, maxnumeqmax 
     86      INTEGER ::   jm             ! current reference number of equation 
     87      INTEGER ::   jm_mint, jm_maxt 
    8888      INTEGER ::   iconv          ! number of iterations in iterative procedure 
    8989      INTEGER ::   iconv_max = 50 ! max number of iterations in iterative procedure 
    9090       
    91       INTEGER, DIMENSION(jpij) ::   numeqmin   ! reference number of top equation 
    92       INTEGER, DIMENSION(jpij) ::   numeqmax   ! reference number of bottom equation 
     91      INTEGER, DIMENSION(jpij) ::   jm_min    ! reference number of top equation 
     92      INTEGER, DIMENSION(jpij) ::   jm_max    ! reference number of bottom equation 
    9393       
    9494      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system 
     
    113113      REAL(wp), DIMENSION(jpij)     ::   zfsw        ! solar radiation absorbed at the surface 
    114114      REAL(wp), DIMENSION(jpij)     ::   zqns_ice_b  ! solar radiation absorbed at the surface 
    115       REAL(wp), DIMENSION(jpij)     ::   zf          ! surface flux function 
     115      REAL(wp), DIMENSION(jpij)     ::   zfnet       ! surface flux function 
    116116      REAL(wp), DIMENSION(jpij)     ::   zdqns_ice_b ! derivative of the surface flux function 
    117117      REAL(wp), DIMENSION(jpij)     ::   zftrice     ! solar radiation transmitted through the ice 
     
    176176      ! 2) Radiation 
    177177      !------------- 
    178       ! 
    179178      z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0 
    180179      DO ji = 1, nidx 
     
    224223      iconv    =  0          ! number of iterations 
    225224      zdti_max =  1000._wp   ! maximal value of error on all points 
    226       !                                                          !----------------------------! 
     225      !                                                          !============================! 
    227226      DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max )   ! Iterative procedure begins ! 
    228          !                                                       !----------------------------! 
    229          ! 
     227         !                                                       !============================! 
    230228         iconv = iconv + 1 
    231229         ! 
     
    272270         ! Used in mono-category case only to simulate an ITD implicitly 
    273271         ! Fichefet and Morales Maqueda, JGR 1997 
    274  
    275272         zghe(1:nidx) = 1._wp 
    276           
     273         ! 
    277274         SELECT CASE ( nn_monocat ) 
    278275 
     
    281278            zepsilon = 0.1_wp 
    282279            DO ji = 1, nidx 
    283     
    284                ! Mean sea ice thermal conductivity 
    285                zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp ) 
    286  
    287                ! Effective thickness he (zhe) 
    288                zhe  = ( rn_cnd_s * ht_i_1d(ji) + zcnd_i * ht_s_1d(ji) ) / ( rn_cnd_s + zcnd_i ) 
    289  
    290                ! G(he) 
     280               zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp )                             ! Mean sea ice thermal conductivity 
     281               zhe = ( rn_cnd_s * ht_i_1d(ji) + zcnd_i * ht_s_1d(ji) ) / ( rn_cnd_s + zcnd_i )   ! Effective thickness he (zhe) 
    291282               IF( zhe >=  zepsilon * 0.5_wp * EXP(1._wp) ) THEN 
    292                   zghe(ji) =  MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) ) 
     283                  zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) )    ! G(he) 
    293284               ENDIF 
    294     
    295285            END DO 
    296286 
     
    303293         DO jk = 0, nlay_s-1 
    304294            DO ji = 1, nidx 
    305                zkappa_s(ji,jk)  = zghe(ji) * rn_cnd_s * z1_h_s(ji) 
     295               zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 
    306296            END DO 
    307297         END DO 
     
    322312         END DO 
    323313         DO ji = 1, nidx   ! Snow-ice interface 
    324             zkappa_i(ji,0)     = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 
     314            zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 
    325315         END DO 
    326316         ! 
     
    337327         DO jk = 1, nlay_s 
    338328            DO ji = 1, nidx 
    339                zeta_s(ji,jk)  = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji) 
     329               zeta_s(ji,jk) = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji) 
    340330            END DO 
    341331         END DO 
     
    352342 
    353343         DO ji = 1, nidx 
    354             zf(ji) = zfsw(ji) + qns_ice_1d(ji) ! incoming = net absorbed solar radiation + non solar total flux (LWup, LWdw, SH, LH) 
     344            zfnet(ji) = zfsw(ji) + qns_ice_1d(ji) ! incoming = net absorbed solar radiation + non solar total flux (LWup, LWdw, SH, LH) 
    355345         END DO 
    356346         ! 
     
    359349         !---------------------------- 
    360350         !!layer denotes the number of the layer in the snow or in the ice 
    361          !!numeq denotes the reference number of the equation in the tridiagonal 
     351         !!jm denotes the reference number of the equation in the tridiagonal 
    362352         !!system, terms of tridiagonal system are indexed as following : 
    363353         !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 
    364354 
    365355         !!ice interior terms (top equation has the same form as the others) 
    366  
    367          DO numeq=1,nlay_i+3 
    368             DO ji = 1, nidx 
    369                ztrid(ji,numeq,1) = 0. 
    370                ztrid(ji,numeq,2) = 0. 
    371                ztrid(ji,numeq,3) = 0. 
    372                zindterm(ji,numeq)= 0. 
    373                zindtbis(ji,numeq)= 0. 
    374                zdiagbis(ji,numeq)= 0. 
    375             ENDDO 
     356         ztrid   (1:nidx,:,:) = 0._wp 
     357         zindterm(1:nidx,:)   = 0._wp 
     358         zindtbis(1:nidx,:)   = 0._wp 
     359         zdiagbis(1:nidx,:)   = 0._wp 
     360 
     361         DO jm = nlay_s + 2, nlay_s + nlay_i  
     362            DO ji = 1, nidx 
     363               jk = jm - nlay_s - 1 
     364               ztrid(ji,jm,1)  = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) 
     365               ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 
     366               ztrid(ji,jm,3)  = - zeta_i(ji,jk) * zkappa_i(ji,jk) 
     367               zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 
     368            END DO 
    376369         ENDDO 
    377370 
    378          DO numeq = nlay_s + 2, nlay_s + nlay_i  
    379             DO ji = 1, nidx 
    380                jk                 = numeq - nlay_s - 1 
    381                ztrid(ji,numeq,1)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk-1) 
    382                ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 
    383                ztrid(ji,numeq,3)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk) 
    384                zindterm(ji,numeq) =  ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 
    385             END DO 
    386          ENDDO 
    387  
    388          numeq =  nlay_s + nlay_i + 1 
     371         jm =  nlay_s + nlay_i + 1 
    389372         DO ji = 1, nidx 
    390373            !!ice bottom term 
    391             ztrid(ji,numeq,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)    
    392             ztrid(ji,numeq,2)  = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) ) 
    393             ztrid(ji,numeq,3)  = 0.0 
    394             zindterm(ji,numeq) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  & 
    395                &                  ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) )  
     374            ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)    
     375            ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) ) 
     376            ztrid(ji,jm,3)  = 0.0 
     377            zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  & 
     378               &            ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) )  
    396379         ENDDO 
    397380 
     
    401384            IF ( ht_s_1d(ji) > 0.0 ) THEN   !  snow-covered cells ! 
    402385               !                            !---------------------! 
    403                ! 
    404                !!snow interior terms (bottom equation has the same form as the others) 
    405                DO numeq = 3, nlay_s + 1 
    406                   jk                  =  numeq - 1 
    407                   ztrid(ji,numeq,1)   =  - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 
    408                   ztrid(ji,numeq,2)   =  1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 
    409                   ztrid(ji,numeq,3)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk) 
    410                   zindterm(ji,numeq)  =  ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 
     386               ! snow interior terms (bottom equation has the same form as the others) 
     387               DO jm = 3, nlay_s + 1 
     388                  jk = jm - 1 
     389                  ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 
     390                  ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 
     391                  ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk) 
     392                  zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 
    411393               END DO 
    412394 
    413                !!case of only one layer in the ice (ice equation is altered) 
     395               ! case of only one layer in the ice (ice equation is altered) 
    414396               IF ( nlay_i == 1 ) THEN 
    415                   ztrid(ji,nlay_s+2,3)    = 0.0 
    416                   zindterm(ji,nlay_s+2)   = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji)  
     397                  ztrid(ji,nlay_s+2,3)  = 0.0 
     398                  zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji)  
    417399               ENDIF 
    418400 
    419401               IF ( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting 
    420402 
    421                   numeqmin(ji)    = 1 
    422                   numeqmax(ji)    = nlay_i + nlay_s + 1 
    423  
    424                   !!surface equation 
     403                  jm_min(ji) = 1 
     404                  jm_max(ji) = nlay_i + nlay_s + 1 
     405 
     406                  ! surface equation 
    425407                  ztrid(ji,1,1)  = 0.0 
    426408                  ztrid(ji,1,2)  = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0) 
    427409                  ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0) 
    428                   zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zf(ji) 
    429  
    430                   !!first layer of snow equation 
    431                   ztrid(ji,2,1)  =  - zkappa_s(ji,0) * zg1s * zeta_s(ji,1) 
    432                   ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 
    433                   ztrid(ji,2,3)  =  - zeta_s(ji,1)* zkappa_s(ji,1) 
    434                   zindterm(ji,2) =  ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1) 
     410                  zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji) 
     411 
     412                  ! first layer of snow equation 
     413                  ztrid(ji,2,1)  = - zkappa_s(ji,0) * zg1s * zeta_s(ji,1) 
     414                  ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 
     415                  ztrid(ji,2,3)  = - zeta_s(ji,1)* zkappa_s(ji,1) 
     416                  zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1) 
    435417 
    436418               ELSE                            !--  case 2 : surface is melting 
    437419                  ! 
    438                   numeqmin(ji)    = 2 
    439                   numeqmax(ji)    = nlay_i + nlay_s + 1 
    440  
    441                   !!first layer of snow equation 
    442                   ztrid(ji,2,1)  =  0.0 
    443                   ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 
    444                   ztrid(ji,2,3)  =  - zeta_s(ji,1)*zkappa_s(ji,1)  
     420                  jm_min(ji) = 2 
     421                  jm_max(ji) = nlay_i + nlay_s + 1 
     422 
     423                  ! first layer of snow equation 
     424                  ztrid(ji,2,1)  = 0.0 
     425                  ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 
     426                  ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1)  
    445427                  zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  & 
    446                      &             ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) )  
     428                     &           ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) )  
    447429               ENDIF 
    448430            !                               !---------------------! 
     
    452434               IF ( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting 
    453435                  ! 
    454                   numeqmin(ji)      =  nlay_s + 1 
    455                   numeqmax(ji)      =  nlay_i + nlay_s + 1 
    456  
    457                   !!surface equation    
    458                   ztrid(ji,numeqmin(ji),1)   =  0.0 
    459                   ztrid(ji,numeqmin(ji),2)   =  zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1     
    460                   ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1 
    461                   zindterm(ji,numeqmin(ji))  =  zdqns_ice_b(ji)*t_su_1d(ji) - zf(ji) 
    462  
    463                   !!first layer of ice equation 
    464                   ztrid(ji,numeqmin(ji)+1,1) =  - zkappa_i(ji,0) * zg1 * zeta_i(ji,1) 
    465                   ztrid(ji,numeqmin(ji)+1,2) =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 
    466                   ztrid(ji,numeqmin(ji)+1,3) =  - zeta_i(ji,1) * zkappa_i(ji,1)   
    467                   zindterm(ji,numeqmin(ji)+1)=  ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)   
    468  
    469                   !!case of only one layer in the ice (surface & ice equations are altered) 
    470  
     436                  jm_min(ji) = nlay_s + 1 
     437                  jm_max(ji) = nlay_i + nlay_s + 1 
     438 
     439                  ! surface equation    
     440                  ztrid(ji,jm_min(ji),1)  = 0.0 
     441                  ztrid(ji,jm_min(ji),2)  = zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1     
     442                  ztrid(ji,jm_min(ji),3)  = zkappa_i(ji,0)*zg1 
     443                  zindterm(ji,jm_min(ji)) = zdqns_ice_b(ji)*t_su_1d(ji) - zfnet(ji) 
     444 
     445                  ! first layer of ice equation 
     446                  ztrid(ji,jm_min(ji)+1,1)  = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1) 
     447                  ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 
     448                  ztrid(ji,jm_min(ji)+1,3)  = - zeta_i(ji,1) * zkappa_i(ji,1)   
     449                  zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)   
     450 
     451                  ! case of only one layer in the ice (surface & ice equations are altered) 
    471452                  IF ( nlay_i == 1 ) THEN 
    472                      ztrid(ji,numeqmin(ji),1)    =  0.0 
    473                      ztrid(ji,numeqmin(ji),2)    =  zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0 
    474                      ztrid(ji,numeqmin(ji),3)    =  zkappa_i(ji,0) * 2.0 
    475                      ztrid(ji,numeqmin(ji)+1,1)  =  -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1) 
    476                      ztrid(ji,numeqmin(ji)+1,2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 
    477                      ztrid(ji,numeqmin(ji)+1,3)  =  0.0 
    478  
    479                      zindterm(ji,numeqmin(ji)+1) =  ztiold(ji,1) + zeta_i(ji,1) *  & 
    480                         &                          ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) 
     453                     ztrid(ji,jm_min(ji),1)    = 0.0 
     454                     ztrid(ji,jm_min(ji),2)    = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0 
     455                     ztrid(ji,jm_min(ji),3)    = zkappa_i(ji,0) * 2.0 
     456                     ztrid(ji,jm_min(ji)+1,1)  = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1) 
     457                     ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 
     458                     ztrid(ji,jm_min(ji)+1,3)  = 0.0 
     459                     zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) *  & 
     460                        &                      ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) 
    481461                  ENDIF 
    482462 
    483463               ELSE                            !--  case 2 : surface is melting 
    484464 
    485                   numeqmin(ji)    =  nlay_s + 2 
    486                   numeqmax(ji)    =  nlay_i + nlay_s + 1 
    487  
    488                   !!first layer of ice equation 
    489                   ztrid(ji,numeqmin(ji),1)      = 0.0 
    490                   ztrid(ji,numeqmin(ji),2)      = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )   
    491                   ztrid(ji,numeqmin(ji),3)      = - zeta_i(ji,1) * zkappa_i(ji,1) 
    492                   zindterm(ji,numeqmin(ji))     = ztiold(ji,1) + zeta_i(ji,1) *  & 
    493                      &                             ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) )  
    494  
    495                   !!case of only one layer in the ice (surface & ice equations are altered) 
     465                  jm_min(ji)    =  nlay_s + 2 
     466                  jm_max(ji)    =  nlay_i + nlay_s + 1 
     467 
     468                  ! first layer of ice equation 
     469                  ztrid(ji,jm_min(ji),1)  = 0.0 
     470                  ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )   
     471                  ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1) 
     472                  zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) *  & 
     473                     &                    ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) )  
     474 
     475                  ! case of only one layer in the ice (surface & ice equations are altered) 
    496476                  IF ( nlay_i == 1 ) THEN 
    497                      ztrid(ji,numeqmin(ji),1)  = 0.0 
    498                      ztrid(ji,numeqmin(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 
    499                      ztrid(ji,numeqmin(ji),3)  = 0.0 
    500                      zindterm(ji,numeqmin(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  & 
    501                         &                       + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0 
     477                     ztrid(ji,jm_min(ji),1)  = 0.0 
     478                     ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 
     479                     ztrid(ji,jm_min(ji),3)  = 0.0 
     480                     zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  & 
     481                        &                    + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0 
    502482                  ENDIF 
    503483 
    504484               ENDIF 
    505485            ENDIF 
    506  
     486            ! 
     487            zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji)) 
     488            zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2) 
     489            ! 
    507490         END DO 
    508491         ! 
     
    511494         !------------------------------ 
    512495         ! Solve the tridiagonal system with Gauss elimination method. 
    513          ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984.    
    514  
    515          maxnumeqmax = 0 
    516          minnumeqmin = nlay_i+5 
    517  
    518          DO ji = 1, nidx 
    519             zindtbis(ji,numeqmin(ji)) =  zindterm(ji,numeqmin(ji)) 
    520             zdiagbis(ji,numeqmin(ji)) =  ztrid(ji,numeqmin(ji),2) 
    521             minnumeqmin               =  MIN(numeqmin(ji),minnumeqmin) 
    522             maxnumeqmax               =  MAX(numeqmax(ji),maxnumeqmax) 
    523          END DO 
    524  
    525          DO jk = minnumeqmin+1, maxnumeqmax 
    526             DO ji = 1, nidx 
    527                numeq               =  min(max(numeqmin(ji)+1,jk),numeqmax(ji)) 
    528                zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2)  - ztrid(ji,numeq,1) * ztrid(ji,numeq-1,3)  / zdiagbis(ji,numeq-1) 
    529                zindtbis(ji,numeq)  =  zindterm(ji,numeq) - ztrid(ji,numeq,1) * zindtbis(ji,numeq-1) / zdiagbis(ji,numeq-1) 
     496         ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 
     497         jm_maxt = 0 
     498         jm_mint = nlay_i+5 
     499         DO ji = 1, nidx 
     500            jm_mint = MIN(jm_min(ji),jm_mint) 
     501            jm_maxt = MAX(jm_max(ji),jm_maxt) 
     502         END DO 
     503 
     504         DO jk = jm_mint+1, jm_maxt 
     505            DO ji = 1, nidx 
     506               jm = min(max(jm_min(ji)+1,jk),jm_max(ji)) 
     507               zdiagbis(ji,jm) = ztrid(ji,jm,2)  - ztrid(ji,jm,1) * ztrid(ji,jm-1,3)  / zdiagbis(ji,jm-1) 
     508               zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 
    530509            END DO 
    531510         END DO 
     
    533512         DO ji = 1, nidx 
    534513            ! ice temperatures 
    535             t_i_1d(ji,nlay_i) =  zindtbis(ji,numeqmax(ji)) / zdiagbis(ji,numeqmax(ji)) 
    536          END DO 
    537  
    538          DO numeq = nlay_i + nlay_s, nlay_s + 2, -1 
    539             DO ji = 1, nidx 
    540                jk    =  numeq - nlay_s - 1 
    541                t_i_1d(ji,jk)  =  ( zindtbis(ji,numeq) - ztrid(ji,numeq,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,numeq) 
     514            t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 
     515         END DO 
     516 
     517         DO jm = nlay_i + nlay_s, nlay_s + 2, -1 
     518            DO ji = 1, nidx 
     519               jk    =  jm - nlay_s - 1 
     520               t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
    542521            END DO 
    543522         END DO 
     
    546525            ! snow temperatures       
    547526            IF( ht_s_1d(ji) > 0._wp ) THEN 
    548                t_s_1d(ji,nlay_s) =  ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  & 
    549                   &                 / zdiagbis(ji,nlay_s+1) 
     527               t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  & 
     528                  &                / zdiagbis(ji,nlay_s+1) 
    550529            ENDIF 
    551530            ! surface temperature 
    552531            ztsub(ji) = t_su_1d(ji) 
    553532            IF( t_su_1d(ji) < rt0 ) THEN 
    554                t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) *  & 
    555                   &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 
     533               t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  & 
     534                  &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 
    556535            ENDIF 
    557536         END DO 
     
    599578      DO ji = 1, nidx 
    600579         !                                ! surface ice conduction flux 
    601          fc_su(ji)       =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   & 
    602             &               - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji)) 
     580         fc_su(ji)   =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   & 
     581            &           - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji)) 
    603582         !                                ! bottom ice conduction flux 
    604          fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 
     583         fc_bo_i(ji) =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 
    605584      END DO 
    606585 
     
    637616      DO ji = 1, nidx 
    638617         !--- Conduction fluxes (positive downwards) 
    639          diag_fc_bo_1d(ji)   = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji) 
    640          diag_fc_su_1d(ji)   = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji) 
     618         diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji) 
     619         diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji) 
    641620 
    642621         !--- Snow-ice interfacial temperature (diagnostic SIMIP) 
    643622         zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) 
    644623         IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN 
    645             t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) + & 
     624            t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   & 
    646625               &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac 
    647626         ELSE 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90

    r8534 r8562  
    2929!!gm  It should be probably better to pass these variable in argument of the routine,  
    3030!!gm  rather than having this long list in USE. This will also highlight what is updated, and what is just use. 
    31    USE sbc_ice , ONLY : emp_oce, qns_oce, qsr_oce , qemp_oce ,                             & 
    32       &                 emp_ice, qsr_ice, qemp_ice, qevap_ice, alb_ice, tn_ice, cldf_ice,  & 
    33       &                 snwice_mass, snwice_mass_b, snwice_fmass 
    34    USE sbc_oce , ONLY : nn_fsbc, ln_ice_embd, sfx, fr_i, qsr_tot, qns, qsr, fmmflx, emp, taum, utau, vtau 
     31   USE sbc_ice 
     32   USE sbc_oce 
    3533!!gm end 
    3634   USE sbccpl         ! Surface boundary condition: coupled interface 
     
    107105      REAL(wp) ::   zqmass           ! Heat flux associated with mass exchange ice->ocean (W.m-2) 
    108106      REAL(wp) ::   zqsr             ! New solar flux received by the ocean 
     107      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d                  ! 2D workspace 
    109108      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_cs, zalb_os     ! 3D workspace 
    110109      !!--------------------------------------------------------------------- 
     
    208207      ENDIF 
    209208      ! 
     209      ! output all fluxes 
     210      !------------------ 
     211      IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * ( 1._wp - at_i_b(:,:) ) )                      !     solar flux at ocean surface 
     212      IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + qemp_oce(:,:) )      ! non-solar flux at ocean surface 
     213      IF( iom_use('qsr_ice') )   CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux at ice surface 
     214      IF( iom_use('qns_ice') )   CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) ! non-solar flux at ice surface 
     215      IF( iom_use('qtr_ice') )   CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux transmitted thru ice 
     216      IF( iom_use('qt_oce' ) )   CALL iom_put( "qt_oce"  , ( qsr_oce(:,:) + qns_oce(:,:) ) * ( 1._wp - at_i_b(:,:) ) + qemp_oce(:,:) ) 
     217      IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) )   & 
     218         &                                                      * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) 
     219!!gm I don't understand the variable below.... why not multiplied by a_i_b or (1-a_i_b) ???  
     220      IF( iom_use('qemp_oce') )  CALL iom_put( "qemp_oce" , qemp_oce(:,:) )   
     221      IF( iom_use('qemp_ice') )  CALL iom_put( "qemp_ice" , qemp_ice(:,:) )   
     222      IF( iom_use('emp_oce' ) )  CALL iom_put( "emp_oce"  , emp_oce (:,:) )   ! emp over ocean (taking into account the snow blown away from the ice) 
     223      IF( iom_use('emp_ice' ) )  CALL iom_put( "emp_ice"  , emp_ice (:,:) )   ! emp over ice   (taking into account the snow blown away from the ice) 
     224 
     225      CALL iom_put( "snowpre"     , sprecip * 86400.    )        ! snow precipitation [m/day] 
     226 
     227      CALL iom_put( "sfxbog"      , sfx_bog             )        ! salt flux from bottom growth 
     228      CALL iom_put( "sfxbom"      , sfx_bom             )        ! salt flux from bottom melting 
     229      CALL iom_put( "sfxsum"      , sfx_sum             )        ! salt flux from surface melting 
     230      CALL iom_put( "sfxlam"      , sfx_lam             )        ! salt flux from lateral melting 
     231      CALL iom_put( "sfxsni"      , sfx_sni             )        ! salt flux from snow ice formation 
     232      CALL iom_put( "sfxopw"      , sfx_opw             )        ! salt flux from open water formation 
     233      CALL iom_put( "sfxdyn"      , sfx_dyn             )        ! salt flux from ridging rafting 
     234      CALL iom_put( "sfxres"      , sfx_res             )        ! salt flux from corrections (resultant) 
     235      CALL iom_put( "sfxbri"      , sfx_bri             )        ! salt flux from brines 
     236      CALL iom_put( "sfxsub"      , sfx_sub             )        ! salt flux from sublimation 
     237      CALL iom_put( "sfx"         , sfx                 )        ! total salt flux 
     238 
     239      CALL iom_put( "vfxres"     , wfx_res              )        ! prod./melting due to corrections  
     240      CALL iom_put( "vfxopw"     , wfx_opw              )        ! lateral thermodynamic ice production 
     241      CALL iom_put( "vfxsni"     , wfx_sni              )        ! snowice ice production 
     242      CALL iom_put( "vfxbog"     , wfx_bog              )        ! bottom thermodynamic ice production 
     243      CALL iom_put( "vfxdyn"     , wfx_dyn              )        ! dynamic ice production (rid/raft) 
     244      CALL iom_put( "vfxsum"     , wfx_sum              )        ! surface melt  
     245      CALL iom_put( "vfxbom"     , wfx_bom              )        ! bottom melt  
     246      CALL iom_put( "vfxlam"     , wfx_lam              )        ! lateral melt  
     247      CALL iom_put( "vfxice"     , wfx_ice              )        ! total ice growth/melt  
     248      IF ( ln_pnd )   CALL iom_put( "vfxpnd", wfx_pnd   )        ! melt pond water flux 
     249 
     250      IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations   
     251         WHERE( htm_i(:,:) < 0.2 .AND. htm_i(:,:) > 0. ) ; z2d = wfx_bog 
     252         ELSEWHERE                                       ; z2d = 0._wp 
     253         END WHERE 
     254         CALL iom_put( "vfxthin", ( wfx_opw + z2d )     ) 
     255      ENDIF 
     256 
     257      CALL iom_put( "vfxspr"     , wfx_spr              )        ! precip (snow) 
     258      CALL iom_put( "vfxsnw"     , wfx_snw              )        ! total snw growth/melt  
     259      CALL iom_put( "vfxsub"     , wfx_sub              )        ! sublimation (snow/ice)  
     260      CALL iom_put( "vfxsub_err" , wfx_err_sub          )        ! "excess" of sublimation sent to ocean       
     261  
     262      CALL iom_put ('hfxthd'     , hfx_thd(:,:)         )   !   
     263      CALL iom_put ('hfxdyn'     , hfx_dyn(:,:)         )   !   
     264      CALL iom_put ('hfxres'     , hfx_res(:,:)         )   !   
     265      CALL iom_put ('hfxout'     , hfx_out(:,:)         )   !   
     266      CALL iom_put ('hfxin'      , hfx_in(:,:)          )   !   
     267      CALL iom_put ('hfxsnw'     , hfx_snw(:,:)         )   !   
     268      CALL iom_put ('hfxsub'     , hfx_sub(:,:)         )   !   
     269      CALL iom_put ('hfxerr'     , hfx_err_dif(:,:)     )   !   
     270      CALL iom_put ('hfxerr_rem' , hfx_err_rem(:,:)     )   !   
     271       
     272      CALL iom_put ('hfxsum'     , hfx_sum(:,:)         )   !   
     273      CALL iom_put ('hfxbom'     , hfx_bom(:,:)         )   !   
     274      CALL iom_put ('hfxbog'     , hfx_bog(:,:)         )   !   
     275      CALL iom_put ('hfxdif'     , hfx_dif(:,:)         )   !   
     276      CALL iom_put ('hfxopw'     , hfx_opw(:,:)         )   !   
     277      CALL iom_put ('hfxtur'     , fhtur(:,:) * at_i_b(:,:) ) ! turbulent heat flux at ice base  
     278      CALL iom_put ('hfxdhc'     , diag_heat(:,:)       )   ! Heat content variation in snow and ice  
     279      CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip  
     280      ! 
    210281      ! controls 
     282      !--------- 
    211283      IF( ln_icediachk .AND. .NOT. ln_bdy)   CALL ice_cons_final('iceupdate')                                       ! conservation 
    212284      IF( ln_icectl                      )   CALL ice_prt       (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints 
     
    282354      !                                      !==  every ocean time-step  ==! 
    283355      ! 
    284       DO jj = 2, jpjm1                                !* update the stress WITHOUT a ice-ocean rotation angle 
     356      DO jj = 2, jpjm1                                !* update the stress WITHOUT an ice-ocean rotation angle 
    285357         DO ji = fs_2, fs_jpim1   ! Vect. Opt. 
    286358            zat_u  = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5_wp   ! ice area at u and V-points 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90

    r8534 r8562  
    4040CONTAINS 
    4141 
    42    SUBROUTINE ice_wri( kindic ) 
     42   SUBROUTINE ice_wri( kt ) 
    4343      !!------------------------------------------------------------------- 
    4444      !!  This routine computes the average of some variables and write it 
     
    5050      !!  modif : 03/06/98 
    5151      !!------------------------------------------------------------------- 
    52       INTEGER, INTENT(in) ::   kindic   ! if kindic < 0 there has been an error somewhere 
     52      INTEGER, INTENT(in) ::   kt   ! time-step 
    5353      ! 
    5454      INTEGER  ::  ji, jj, jk, jl  ! dummy loop indices 
    55       REAL(wp) ::  z2da, z2db, ztmp, zrho1, zrho2, zmiss_val 
     55      REAL(wp) ::  z2da, z2db, zrho1, zrho2, zmiss_val 
    5656      REAL(wp), DIMENSION(jpi,jpj)     ::  z2d, zswi, zmiss !  2D workspace 
    5757      REAL(wp), DIMENSION(jpi,jpj)     ::  zfb              ! ice freeboard 
     
    9595      ! Standard outputs 
    9696      !---------------------------------------- 
    97       ! fluxes  
    98       IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * ( 1._wp - at_i_b(:,:) ) )                      !     solar flux at ocean surface 
    99       IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + qemp_oce(:,:) )      ! non-solar flux at ocean surface 
    100       IF( iom_use('qsr_ice') )   CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux at ice surface 
    101       IF( iom_use('qns_ice') )   CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) ! non-solar flux at ice surface 
    102       IF( iom_use('qtr_ice') )   CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux transmitted thru ice 
    103       IF( iom_use('qt_oce' ) )   CALL iom_put( "qt_oce"  , ( qsr_oce(:,:) + qns_oce(:,:) ) * ( 1._wp - at_i_b(:,:) ) + qemp_oce(:,:) ) 
    104       IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) )   & 
    105          &                                                      * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) 
    106 !!gm I don't understand the variable below.... why not multiplied by a_i_b or (1-a_i_b) ???  
    107       IF( iom_use('qemp_oce') )  CALL iom_put( "qemp_oce" , qemp_oce(:,:) )   
    108       IF( iom_use('qemp_ice') )  CALL iom_put( "qemp_ice" , qemp_ice(:,:) )   
    109       IF( iom_use('emp_oce' ) )  CALL iom_put( "emp_oce"  , emp_oce (:,:) )   ! emp over ocean (taking into account the snow blown away from the ice) 
    110       IF( iom_use('emp_ice' ) )  CALL iom_put( "emp_ice"  , emp_ice (:,:) )   ! emp over ice   (taking into account the snow blown away from the ice) 
    111  
    11297      ! velocity 
    11398      IF( iom_use('uice_ipa') )  CALL iom_put( "uice_ipa" , u_ice         )   ! ice velocity u component 
     
    126111         IF( iom_use('icevel_mv') )   CALL iom_put( "icevel_mv"    , z2d(:,:) * zswi(:,:) + zmiss(:,:) )   ! ice velocity module (missing value) 
    127112      ENDIF 
    128  
    129       IF( iom_use('tau_icebfr') )     CALL iom_put( "tau_icebfr"  , tau_icebfr             )  ! ice friction with ocean bottom (landfast ice)   
    130113      ! 
    131114      IF( iom_use('miceage')  )       CALL iom_put( "miceage"     , om_i * zswi * zamask15 )  ! mean ice age 
    132115      IF( iom_use('micet')    )       CALL iom_put( "micet"       , ( tm_i  - rt0 ) * zswi )  ! ice mean    temperature 
    133116      IF( iom_use('icest')    )       CALL iom_put( "icest"       , ( tm_su - rt0 ) * zswi )  ! ice surface temperature 
    134       IF( iom_use('icecolf')  )       CALL iom_put( "icecolf"     , hicol                  )  ! frazil ice collection thickness 
     117      IF( iom_use('icecolf')  )       CALL iom_put( "icecolf"     , ht_i_new               )  ! new ice thickness formed in the leads 
    135118      ! 
    136119      CALL iom_put( "isst"        , sst_m               )        ! sea surface temperature 
     
    142125      CALL iom_put( "isnowhc"     , et_s  * zswi        )        ! snow total heat content 
    143126      CALL iom_put( "ibrinv"      , bvm_i * zswi * 100. )        ! brine volume 
    144       CALL iom_put( "snowpre"     , sprecip * 86400.    )        ! snow precipitation [m/day] 
    145127      CALL iom_put( "micesalt"    , smt_i   * zswi      )        ! mean ice salinity 
    146128      CALL iom_put( "snowvol"     , vt_s    * zswi      )        ! snow volume 
    147129       
    148       CALL iom_put( "sfxbog"      , sfx_bog * rday      )        ! salt flux from bottom growth 
    149       CALL iom_put( "sfxbom"      , sfx_bom * rday      )        ! salt flux from bottom melting 
    150       CALL iom_put( "sfxsum"      , sfx_sum * rday      )        ! salt flux from surface melting 
    151       CALL iom_put( "sfxlam"      , sfx_lam * rday      )        ! salt flux from lateral melting 
    152       CALL iom_put( "sfxsni"      , sfx_sni * rday      )        ! salt flux from snow ice formation 
    153       CALL iom_put( "sfxopw"      , sfx_opw * rday      )        ! salt flux from open water formation 
    154       CALL iom_put( "sfxdyn"      , sfx_dyn * rday      )        ! salt flux from ridging rafting 
    155       CALL iom_put( "sfxres"      , sfx_res * rday      )        ! salt flux from corrections (resultant) 
    156       CALL iom_put( "sfxbri"      , sfx_bri * rday      )        ! salt flux from brines 
    157       CALL iom_put( "sfxsub"      , sfx_sub * rday      )        ! salt flux from sublimation 
    158       CALL iom_put( "sfx"         , sfx     * rday      )        ! total salt flux 
    159  
    160       ztmp = rday / rhoic 
    161       CALL iom_put( "vfxres"     , wfx_res * ztmp       )        ! daily prod./melting due to corrections  
    162       CALL iom_put( "vfxopw"     , wfx_opw * ztmp       )        ! daily lateral thermodynamic ice production 
    163       CALL iom_put( "vfxsni"     , wfx_sni * ztmp       )        ! daily snowice ice production 
    164       CALL iom_put( "vfxbog"     , wfx_bog * ztmp       )        ! daily bottom thermodynamic ice production 
    165       CALL iom_put( "vfxdyn"     , wfx_dyn * ztmp       )        ! daily dynamic ice production (rid/raft) 
    166       CALL iom_put( "vfxsum"     , wfx_sum * ztmp       )        ! surface melt  
    167       CALL iom_put( "vfxbom"     , wfx_bom * ztmp       )        ! bottom melt  
    168       CALL iom_put( "vfxlam"     , wfx_lam * ztmp       )        ! lateral melt  
    169       CALL iom_put( "vfxice"     , wfx_ice * ztmp       )        ! total ice growth/melt  
    170  
    171       IF ( ln_pnd ) & 
    172          CALL iom_put( "vfxpnd"  , wfx_pnd * ztmp       )        ! melt pond water flux 
    173  
    174       IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations   
    175          WHERE( htm_i(:,:) < 0.2 .AND. htm_i(:,:) > 0. ) ; z2d = wfx_bog 
    176          ELSEWHERE                                       ; z2d = 0._wp 
    177          END WHERE 
    178          CALL iom_put( "vfxthin", ( wfx_opw + z2d ) * ztmp ) 
    179       ENDIF 
    180  
    181       ztmp = rday * r1_rhosn 
    182       CALL iom_put( "vfxspr"     , wfx_spr * ztmp       )        ! precip (snow) 
    183       CALL iom_put( "vfxsnw"     , wfx_snw * ztmp       )        ! total snw growth/melt  
    184       CALL iom_put( "vfxsub"     , wfx_sub * ztmp       )        ! sublimation (snow/ice)  
    185       CALL iom_put( "vfxsub_err" , wfx_err_sub * ztmp   )        ! "excess" of sublimation sent to ocean       
    186   
    187       CALL iom_put ('hfxthd'     , hfx_thd(:,:)         )   !   
    188       CALL iom_put ('hfxdyn'     , hfx_dyn(:,:)         )   !   
    189       CALL iom_put ('hfxres'     , hfx_res(:,:)         )   !   
    190       CALL iom_put ('hfxout'     , hfx_out(:,:)         )   !   
    191       CALL iom_put ('hfxin'      , hfx_in(:,:)          )   !   
    192       CALL iom_put ('hfxsnw'     , hfx_snw(:,:)         )   !   
    193       CALL iom_put ('hfxsub'     , hfx_sub(:,:)         )   !   
    194       CALL iom_put ('hfxerr'     , hfx_err_dif(:,:)     )   !   
    195       CALL iom_put ('hfxerr_rem' , hfx_err_rem(:,:)     )   !   
    196        
    197       CALL iom_put ('hfxsum'     , hfx_sum(:,:)         )   !   
    198       CALL iom_put ('hfxbom'     , hfx_bom(:,:)         )   !   
    199       CALL iom_put ('hfxbog'     , hfx_bog(:,:)         )   !   
    200       CALL iom_put ('hfxdif'     , hfx_dif(:,:)         )   !   
    201       CALL iom_put ('hfxopw'     , hfx_opw(:,:)         )   !   
    202       CALL iom_put ('hfxtur'     , fhtur(:,:) * at_i_b(:,:) ) ! turbulent heat flux at ice base  
    203       CALL iom_put ('hfxdhc'     , diag_heat(:,:)       )   ! Heat content variation in snow and ice  
    204       CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip  
    205  
    206 !!gm ====>>>>>  THIS should be moved where at_ip, vt_ip are computed fro the last time in the time-step  (limmpd) 
    207       ! MV MP 2016 
    208130      IF ( ln_pnd ) THEN 
    209131         CALL iom_put( "iceamp"  , at_ip  * zswi        )   ! melt pond total fraction 
    210132         CALL iom_put( "icevmp"  , vt_ip  * zswi        )   ! melt pond total volume per unit area 
    211133      ENDIF 
    212       ! END MV MP 2016 
    213 !!gm  <<<<<<======= end 
    214134 
    215135      !---------------------------------- 
     
    225145      IF ( iom_use('brinevol_cat') )  CALL iom_put( "brinevol_cat", bv_i * 100. * zswi2 )          ! brine volume 
    226146 
    227       ! MV MP 2016 
    228147      IF ( ln_pnd ) THEN 
    229148         IF ( iom_use('iceamp_cat') )  CALL iom_put( "iceamp_cat"     , a_ip       * zswi2   )       ! melt pond frac for categories 
     
    232151         IF ( iom_use('iceafp_cat') )  CALL iom_put( "iceafp_cat"     , a_ip_frac  * zswi2   )       ! melt pond frac for categories 
    233152      ENDIF 
    234       ! END MV MP 2016 
    235153 
    236154      !-------------------------------- 
     
    360278      !! History :   4.0  !  2013-06  (C. Rousset) 
    361279      !!---------------------------------------------------------------------- 
    362       INTEGER, INTENT( in )   ::   kt               ! ocean time-step index) 
     280      INTEGER, INTENT( in )   ::   kt               ! ocean time-step index 
    363281      INTEGER, INTENT( in )   ::   kid , kh_i 
    364282      INTEGER                 ::   nz_i, jl 
     
    476394      CALL histwrite( kid, "sntemcat", kt, tm_su - rt0 , jpi*jpj*jpl, (/1/) )     
    477395 
    478       ! Close the file 
    479       ! ----------------- 
    480 !!gm I don't understand why the file is not closed ! 
    481       !CALL histclo( kid ) 
     396      !! The file is closed in dia_wri_state (ocean routine) 
     397      !! CALL histclo( kid ) 
    482398      ! 
    483399    END SUBROUTINE ice_wri_state 
Note: See TracChangeset for help on using the changeset viewer.