Changeset 8562 for branches/2017/dev_r8183_ICEMODEL
- Timestamp:
- 2017-09-25T21:11:19+02:00 (7 years ago)
- 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 129 129 <field id="albedo" long_name="Mean albedo over sea ice and ocean" unit="" /> 130 130 131 131 <field id="iceamp" long_name="melt pond fraction" standard_name="sea_ice_meltpond_fraction" unit="%" /> 132 132 <field id="icevmp" long_name="melt pond volume" standard_name="sea_ice_meltpond_volume" unit="m" /> 133 133 … … 154 154 <field id="tau_icebfr" long_name="ice friction on ocean bottom for landfast ice" unit="N/2" /> 155 155 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" /> 159 159 <field id="deitrp" long_name="advected ice enthalpy" unit="W/m2" /> 160 160 <field id="destrp" long_name="advected snw enthalpy" unit="W/m2" /> 161 161 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="d aily 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" /> 188 188 189 189 <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 172 172 ! !!** ice-dynamics namelist (namdyn) ** 173 173 REAL(wp), PUBLIC :: rn_ishlat !: lateral boundary condition for sea-ice 174 REAL(wp), PUBLIC :: rn_cio !: drag coefficient for oceanic stress175 174 LOGICAL , PUBLIC :: ln_landfast !: landfast ice parameterization (T or F) 176 175 REAL(wp), PUBLIC :: rn_gamma !: fraction of ocean depth that ice must reach to initiate landfast ice … … 184 183 REAL(wp), PUBLIC :: rn_relast !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 185 184 ! 186 ! !!** ice- thermodynamics namelist (namthd) **185 ! !!** ice-surface forcing namelist (namforcing) ** 187 186 ! -- icethd_dh -- ! 188 187 REAL(wp), PUBLIC :: rn_blow_s !: coef. for partitioning of snowfall between leads and sea ice 189 188 ! -- icethd -- ! 189 REAL(wp), PUBLIC :: rn_cio !: drag coefficient for oceanic stress 190 190 INTEGER , PUBLIC :: nn_iceflx !: Redistribute heat flux over ice categories 191 191 ! ! =-1 Do nothing (needs N(cat) fluxes) … … 235 235 ! !!** define arrays 236 236 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: u_oce,v_oce !: surface ocean velocity used in ice dynamics 237 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: h icol!: ice collection thickness accreted in leads237 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_i_new !: ice collection thickness accreted in leads 238 238 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: strength !: ice strength 239 239 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: stress1_i, stress2_i, stress12_i !: 1st, 2nd & diagonal stress tensor element … … 437 437 ! stay within Fortran's max-line length limit. 438 438 ii = 1 439 ALLOCATE( u_oce (jpi,jpj) , v_oce (jpi,jpj) , h icol(jpi,jpj) , &439 ALLOCATE( u_oce (jpi,jpj) , v_oce (jpi,jpj) , ht_i_new (jpi,jpj) , & 440 440 & strength(jpi,jpj) , stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) , & 441 441 & 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 208 208 ii = ii + 1 209 209 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) , & 211 211 & eh_i_old(jpij,0:nlay_i+1) , h_i_old(jpij,0:nlay_i+1) , STAT=ierr(ii) ) 212 212 ! -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedia.F90
r8534 r8562 76 76 IF( kt == nit000 .AND. lwp ) THEN 77 77 WRITE(numout,*) 78 WRITE(numout,*)'icedia : outpout ice diagnostics (integrated over the domain)'78 WRITE(numout,*)'icedia: output ice diagnostics (integrated over the domain)' 79 79 WRITE(numout,*)'~~~~~~' 80 80 ENDIF -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90
r8534 r8562 87 87 WHERE( ht_i(:,:,jl) > ht_n(:,:) * rn_gamma ) tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 88 88 END DO 89 IF( iom_use('tau_icebfr') ) CALL iom_put( 'tau_icebfr', tau_icebfr ) 89 90 ENDIF 90 91 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv.F90
r8534 r8562 105 105 diag_trp_vi (:,:) = SUM( v_i (:,:,:) - v_i_b (:,:,:) , dim=3 ) * r1_rdtice 106 106 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 transport108 IF( iom_use('snwtrp') ) CALL iom_put( "snwtrp" , diag_trp_vs * rday )! snw volume transport109 IF( iom_use('saltrp') ) CALL iom_put( "saltrp" , diag_trp_smv * r day * rhoic ) ! salt content transport110 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) 112 112 113 113 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 128 128 ! utau_ice, vtau_ice = surface ice stress [N/m2] 129 129 !------------------------------------------------! 130 CALL ice_forcing_tau( kt, ksbc, utau_ice, vtau_ice ) 131 130 CALL ice_forcing_tau( kt, ksbc, utau_ice, vtau_ice ) 132 131 !-------------------------------------! 133 132 ! --- ice dynamics and advection --- ! … … 136 135 CALL ice_rst_opn( kt ) ! Open Ice restart file (if necessary) 137 136 ! 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 139 139 ! 140 140 ! !== lateral boundary conditions ==! 141 141 #if defined key_agrif 142 IF( .NOT. Agrif_Root() ) 142 IF( .NOT. Agrif_Root() ) CALL agrif_interp_lim3('T') ! -- AGRIF 143 143 #endif 144 IF( ln_icethd .AND. ln_bdy ) 144 IF( ln_icethd .AND. ln_bdy ) CALL bdy_ice( kt ) ! -- bdy ice thermo 145 145 ! 146 146 ! 147 147 ! !== previous lead fraction and ice volume for flux calculations 148 !149 148 CALL ice_var_glo2eqv ! ht_i and ht_s for ice albedo calculation 150 149 CALL ice_var_agg(1) ! at_i for coupling … … 165 164 ! qprec_ice, qevap_ice 166 165 !------------------------------------------------------! 167 CALL ice_forcing_flx( kt, ksbc ) 168 166 CALL ice_forcing_flx( kt, ksbc ) 169 167 !----------------------------! 170 168 ! --- ice thermodynamics --- ! 171 169 !----------------------------! 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 180 175 # if defined key_agrif 181 IF( .NOT. Agrif_Root() ) CALL agrif_update_lim3( kt )176 IF( .NOT. Agrif_Root() ) CALL agrif_update_lim3( kt ) 182 177 # endif 183 184 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 ! 186 181 !! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work 187 182 !! moreover it should only be called at the update frequency (as in agrif_lim3_update.F90) 188 183 !!# if defined key_agrif 189 !! IF( .NOT. Agrif_Root() ) CALL Agrif_ChildGrid_To_ParentGrid()184 !! IF( .NOT. Agrif_Root() ) CALL Agrif_ChildGrid_To_ParentGrid() 190 185 !!# endif 191 CALL ice_update_flx( kt )! -- Update ocean surface mass, heat and salt fluxes186 CALL ice_update_flx( kt ) ! -- Update ocean surface mass, heat and salt fluxes 192 187 !!# if defined key_agrif 193 !! IF( .NOT. Agrif_Root() ) CALL Agrif_ParentGrid_To_ChildGrid()188 !! IF( .NOT. Agrif_Root() ) CALL Agrif_ParentGrid_To_ChildGrid() 194 189 !!# 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 203 197 ! 204 198 ENDIF ! End sea-ice time step only … … 207 201 ! --- Ocean time step --- ! 208 202 !-------------------------! 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 212 204 !!gm remark, the ocean-ice stress is not saved in ice diag call above ..... find a solution!!! 213 205 ! 214 IF( nn_timing == 1 ) CALL timing_stop('ice_stp')206 IF( nn_timing == 1 ) CALL timing_stop('ice_stp') 215 207 ! 216 208 END SUBROUTINE ice_stp … … 246 238 CALL ice_itd_init ! ice thickness distribution initialization 247 239 ! 248 ! MV MP 2016249 240 CALL lim_mp_init ! set melt ponds parameters (clem: important to be located here) 250 ! END MV MP 2016241 ! 251 242 ! ! Initial sea-ice state 252 243 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 116 116 END SELECT 117 117 118 DO ji = 1, nidx119 icount (ji,:) = 0120 zdh_s_mel(ji) = 0._wp121 e_i_1d(ji,nlay_i+1) = 0._wp ! Initialize enthalpy at nlay_i+1122 END DO123 124 118 ! initialize layer thicknesses and enthalpies 125 119 h_i_old (1:nidx,0:nlay_i+1) = 0._wp … … 227 221 ! If heat still available (zq_su > 0), then melt more snow 228 222 zdeltah(1:nidx,:) = 0._wp 223 zdh_s_mel(1:nidx) = 0._wp 229 224 DO jk = 1, nlay_s 230 225 DO ji = 1, nidx … … 443 438 444 439 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)447 440 448 441 END DO … … 475 468 476 469 ! 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) 478 471 h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji) 479 472 … … 654 647 655 648 ! --- 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 659 650 ! 660 651 END SUBROUTINE ice_thd_dh -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_do.F90
r8559 r8562 123 123 ! 3) Collection thickness of ice formed in leads and polynyas 124 124 !------------------------------------------------------------------------------! 125 ! h icolis the thickness of new ice formed in open water126 ! h icolcan 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) 127 127 ! Frazil ice forms in open water, is transported by wind 128 128 ! accumulates at the edge of the consolidated ice edge … … 136 136 137 137 ! Default new ice thickness 138 WHERE( qlead(:,:) < 0._wp ) ; h icol(:,:) = rn_hinew139 ELSEWHERE ; h icol(:,:) = 0._wp138 WHERE( qlead(:,:) < 0._wp ) ; ht_i_new(:,:) = rn_hinew 139 ELSEWHERE ; ht_i_new(:,:) = 0._wp 140 140 END WHERE 141 141 … … 145 145 ! Physical constants 146 146 !-------------------- 147 h icol(:,:) = 0._wp147 ht_i_new(:,:) = 0._wp 148 148 149 149 zhicrit = 0.04 ! frazil ice thickness … … 192 192 ! Iterative procedure 193 193 !--------------------- 194 h icol(ji,jj) = zhicrit + ( zhicrit + 0.1 ) &194 ht_i_new(ji,jj) = zhicrit + ( zhicrit + 0.1 ) & 195 195 & / ( ( zhicrit + 0.1 ) * ( zhicrit + 0.1 ) - zhicrit * zhicrit ) * ztwogp * zvrel2 196 196 197 197 iter = 1 198 198 DO WHILE ( iter < 20 ) 199 zf = ( h icol(ji,jj) - zhicrit ) * ( hicol(ji,jj) * hicol(ji,jj) - zhicrit * zhicrit ) - &200 & h icol(ji,jj) * zhicrit * ztwogp * zvrel2201 zfp = ( h icol(ji,jj) - zhicrit ) * ( 3.0 * hicol(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2202 203 h icol(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 ) 204 204 iter = iter + 1 205 205 END DO … … 210 210 END DO 211 211 ! 212 CALL lbc_lnk_multi( zvrel, 'T', 1., h icol, 'T', 1. )212 CALL lbc_lnk_multi( zvrel, 'T', 1., ht_i_new, 'T', 1. ) 213 213 214 214 ENDIF ! End of computation of frazil ice collection thickness … … 252 252 CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx) , sfx_opw ) 253 253 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) , h icol)254 CALL tab_2d_1d( nidx, idxice(1:nidx), zh_newice (1:nidx) , ht_i_new ) 255 255 CALL tab_2d_1d( nidx, idxice(1:nidx), zvrel_1d (1:nidx) , zvrel ) 256 256 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_ent.F90
r8534 r8562 78 78 ! 1) Cumulative integral of old enthalpy * thickness and layers interfaces 79 79 !-------------------------------------------------------------------------- 80 zeh_cum0( :,0:nlay_i+2) = 0._wp81 zh_cum0 ( :,0:nlay_i+2) = 0._wp80 zeh_cum0(1:nidx,0) = 0._wp 81 zh_cum0 (1:nidx,0) = 0._wp 82 82 DO jk0 = 1, nlay_i+2 83 83 DO ji = 1, nidx … … 96 96 97 97 ! new layers interfaces 98 zh_cum1( :,0:nlay_i) = 0._wp98 zh_cum1(1:nidx,0) = 0._wp 99 99 DO jk1 = 1, nlay_i 100 100 DO ji = 1, nidx … … 103 103 END DO 104 104 105 zeh_cum1( :,0:nlay_i) = 0._wp105 zeh_cum1(1:nidx,0) = 0._wp 106 106 ! new cumulative q*h => linear interpolation 107 107 DO jk0 = 1, nlay_i+1 … … 117 117 END DO 118 118 ! 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) 120 120 121 121 ! new enthalpies -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90
r8534 r8562 84 84 !!------------------------------------------------------------------- 85 85 INTEGER :: ji, jk ! spatial loop index 86 INTEGER :: numeq! current reference number of equation87 INTEGER :: minnumeqmin, maxnumeqmax86 INTEGER :: jm ! current reference number of equation 87 INTEGER :: jm_mint, jm_maxt 88 88 INTEGER :: iconv ! number of iterations in iterative procedure 89 89 INTEGER :: iconv_max = 50 ! max number of iterations in iterative procedure 90 90 91 INTEGER, DIMENSION(jpij) :: numeqmin! reference number of top equation92 INTEGER, DIMENSION(jpij) :: numeqmax! reference number of bottom equation91 INTEGER, DIMENSION(jpij) :: jm_min ! reference number of top equation 92 INTEGER, DIMENSION(jpij) :: jm_max ! reference number of bottom equation 93 93 94 94 REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system … … 113 113 REAL(wp), DIMENSION(jpij) :: zfsw ! solar radiation absorbed at the surface 114 114 REAL(wp), DIMENSION(jpij) :: zqns_ice_b ! solar radiation absorbed at the surface 115 REAL(wp), DIMENSION(jpij) :: zf 115 REAL(wp), DIMENSION(jpij) :: zfnet ! surface flux function 116 116 REAL(wp), DIMENSION(jpij) :: zdqns_ice_b ! derivative of the surface flux function 117 117 REAL(wp), DIMENSION(jpij) :: zftrice ! solar radiation transmitted through the ice … … 176 176 ! 2) Radiation 177 177 !------------- 178 !179 178 z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0 180 179 DO ji = 1, nidx … … 224 223 iconv = 0 ! number of iterations 225 224 zdti_max = 1000._wp ! maximal value of error on all points 226 ! ! ----------------------------!225 ! !============================! 227 226 DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max ) ! Iterative procedure begins ! 228 ! !----------------------------! 229 ! 227 ! !============================! 230 228 iconv = iconv + 1 231 229 ! … … 272 270 ! Used in mono-category case only to simulate an ITD implicitly 273 271 ! Fichefet and Morales Maqueda, JGR 1997 274 275 272 zghe(1:nidx) = 1._wp 276 273 ! 277 274 SELECT CASE ( nn_monocat ) 278 275 … … 281 278 zepsilon = 0.1_wp 282 279 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) 291 282 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) 293 284 ENDIF 294 295 285 END DO 296 286 … … 303 293 DO jk = 0, nlay_s-1 304 294 DO ji = 1, nidx 305 zkappa_s(ji,jk) 295 zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 306 296 END DO 307 297 END DO … … 322 312 END DO 323 313 DO ji = 1, nidx ! Snow-ice interface 324 zkappa_i(ji,0) 314 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 325 315 END DO 326 316 ! … … 337 327 DO jk = 1, nlay_s 338 328 DO ji = 1, nidx 339 zeta_s(ji,jk) 329 zeta_s(ji,jk) = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji) 340 330 END DO 341 331 END DO … … 352 342 353 343 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) 355 345 END DO 356 346 ! … … 359 349 !---------------------------- 360 350 !!layer denotes the number of the layer in the snow or in the ice 361 !! numeqdenotes the reference number of the equation in the tridiagonal351 !!jm denotes the reference number of the equation in the tridiagonal 362 352 !!system, terms of tridiagonal system are indexed as following : 363 353 !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 364 354 365 355 !!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 376 369 ENDDO 377 370 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 389 372 DO ji = 1, nidx 390 373 !!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.0394 zindterm(ji, numeq) =ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * &395 & 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) ) 396 379 ENDDO 397 380 … … 401 384 IF ( ht_s_1d(ji) > 0.0 ) THEN ! snow-covered cells ! 402 385 ! !---------------------! 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) 411 393 END DO 412 394 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) 414 396 IF ( nlay_i == 1 ) THEN 415 ztrid(ji,nlay_s+2,3) =0.0416 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) 417 399 ENDIF 418 400 419 401 IF ( t_su_1d(ji) < rt0 ) THEN !-- case 1 : no surface melting 420 402 421 numeqmin(ji) =1422 numeqmax(ji) =nlay_i + nlay_s + 1423 424 ! !surface equation403 jm_min(ji) = 1 404 jm_max(ji) = nlay_i + nlay_s + 1 405 406 ! surface equation 425 407 ztrid(ji,1,1) = 0.0 426 408 ztrid(ji,1,2) = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0) 427 409 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 equation431 ztrid(ji,2,1) = 432 ztrid(ji,2,2) = 433 ztrid(ji,2,3) = 434 zindterm(ji,2) = 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) 435 417 436 418 ELSE !-- case 2 : surface is melting 437 419 ! 438 numeqmin(ji) =2439 numeqmax(ji) =nlay_i + nlay_s + 1440 441 ! !first layer of snow equation442 ztrid(ji,2,1) = 443 ztrid(ji,2,2) = 444 ztrid(ji,2,3) = 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) 445 427 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * & 446 & 428 & ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 447 429 ENDIF 448 430 ! !---------------------! … … 452 434 IF ( t_su_1d(ji) < rt0 ) THEN !-- case 1 : no surface melting 453 435 ! 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) 471 452 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) ) 481 461 ENDIF 482 462 483 463 ELSE !-- case 2 : surface is melting 484 464 485 numeqmin(ji) = nlay_s + 2486 numeqmax(ji) = nlay_i + nlay_s + 1487 488 ! !first layer of ice equation489 ztrid(ji, numeqmin(ji),1) =0.0490 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 & 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) 496 476 IF ( nlay_i == 1 ) THEN 497 ztrid(ji, numeqmin(ji),1) =0.0498 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.0500 zindterm(ji, numeqmin(ji)) =ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) &501 & 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 502 482 ENDIF 503 483 504 484 ENDIF 505 485 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 ! 507 490 END DO 508 491 ! … … 511 494 !------------------------------ 512 495 ! 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) 530 509 END DO 531 510 END DO … … 533 512 DO ji = 1, nidx 534 513 ! 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, -1539 DO ji = 1, nidx 540 jk = numeq- nlay_s - 1541 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) 542 521 END DO 543 522 END DO … … 546 525 ! snow temperatures 547 526 IF( ht_s_1d(ji) > 0._wp ) THEN 548 t_s_1d(ji,nlay_s) = 549 & 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) 550 529 ENDIF 551 530 ! surface temperature 552 531 ztsub(ji) = t_su_1d(ji) 553 532 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)) 556 535 ENDIF 557 536 END DO … … 599 578 DO ji = 1, nidx 600 579 ! ! surface ice conduction flux 601 fc_su(ji) 602 & 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)) 603 582 ! ! bottom ice conduction flux 604 fc_bo_i(ji) 583 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 605 584 END DO 606 585 … … 637 616 DO ji = 1, nidx 638 617 !--- Conduction fluxes (positive downwards) 639 diag_fc_bo_1d(ji) 640 diag_fc_su_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) 641 620 642 621 !--- Snow-ice interfacial temperature (diagnostic SIMIP) 643 622 zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) 644 623 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) + & 646 625 & ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac 647 626 ELSE -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90
r8534 r8562 29 29 !!gm It should be probably better to pass these variable in argument of the routine, 30 30 !!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 35 33 !!gm end 36 34 USE sbccpl ! Surface boundary condition: coupled interface … … 107 105 REAL(wp) :: zqmass ! Heat flux associated with mass exchange ice->ocean (W.m-2) 108 106 REAL(wp) :: zqsr ! New solar flux received by the ocean 107 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace 109 108 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zalb_cs, zalb_os ! 3D workspace 110 109 !!--------------------------------------------------------------------- … … 208 207 ENDIF 209 208 ! 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 ! 210 281 ! controls 282 !--------- 211 283 IF( ln_icediachk .AND. .NOT. ln_bdy) CALL ice_cons_final('iceupdate') ! conservation 212 284 IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints … … 282 354 ! !== every ocean time-step ==! 283 355 ! 284 DO jj = 2, jpjm1 !* update the stress WITHOUT a ice-ocean rotation angle356 DO jj = 2, jpjm1 !* update the stress WITHOUT an ice-ocean rotation angle 285 357 DO ji = fs_2, fs_jpim1 ! Vect. Opt. 286 358 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 40 40 CONTAINS 41 41 42 SUBROUTINE ice_wri( k indic)42 SUBROUTINE ice_wri( kt ) 43 43 !!------------------------------------------------------------------- 44 44 !! This routine computes the average of some variables and write it … … 50 50 !! modif : 03/06/98 51 51 !!------------------------------------------------------------------- 52 INTEGER, INTENT(in) :: k indic ! if kindic < 0 there has been an error somewhere52 INTEGER, INTENT(in) :: kt ! time-step 53 53 ! 54 54 INTEGER :: ji, jj, jk, jl ! dummy loop indices 55 REAL(wp) :: z2da, z2db, z tmp, zrho1, zrho2, zmiss_val55 REAL(wp) :: z2da, z2db, zrho1, zrho2, zmiss_val 56 56 REAL(wp), DIMENSION(jpi,jpj) :: z2d, zswi, zmiss ! 2D workspace 57 57 REAL(wp), DIMENSION(jpi,jpj) :: zfb ! ice freeboard … … 95 95 ! Standard outputs 96 96 !---------------------------------------- 97 ! fluxes98 IF( iom_use('qsr_oce') ) CALL iom_put( "qsr_oce" , qsr_oce(:,:) * ( 1._wp - at_i_b(:,:) ) ) ! solar flux at ocean surface99 IF( iom_use('qns_oce') ) CALL iom_put( "qns_oce" , qns_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + qemp_oce(:,:) ) ! non-solar flux at ocean surface100 IF( iom_use('qsr_ice') ) CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) ! solar flux at ice surface101 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 surface102 IF( iom_use('qtr_ice') ) CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) ! solar flux transmitted thru ice103 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 112 97 ! velocity 113 98 IF( iom_use('uice_ipa') ) CALL iom_put( "uice_ipa" , u_ice ) ! ice velocity u component … … 126 111 IF( iom_use('icevel_mv') ) CALL iom_put( "icevel_mv" , z2d(:,:) * zswi(:,:) + zmiss(:,:) ) ! ice velocity module (missing value) 127 112 ENDIF 128 129 IF( iom_use('tau_icebfr') ) CALL iom_put( "tau_icebfr" , tau_icebfr ) ! ice friction with ocean bottom (landfast ice)130 113 ! 131 114 IF( iom_use('miceage') ) CALL iom_put( "miceage" , om_i * zswi * zamask15 ) ! mean ice age 132 115 IF( iom_use('micet') ) CALL iom_put( "micet" , ( tm_i - rt0 ) * zswi ) ! ice mean temperature 133 116 IF( iom_use('icest') ) CALL iom_put( "icest" , ( tm_su - rt0 ) * zswi ) ! ice surface temperature 134 IF( iom_use('icecolf') ) CALL iom_put( "icecolf" , h icol ) ! frazil ice collection thickness117 IF( iom_use('icecolf') ) CALL iom_put( "icecolf" , ht_i_new ) ! new ice thickness formed in the leads 135 118 ! 136 119 CALL iom_put( "isst" , sst_m ) ! sea surface temperature … … 142 125 CALL iom_put( "isnowhc" , et_s * zswi ) ! snow total heat content 143 126 CALL iom_put( "ibrinv" , bvm_i * zswi * 100. ) ! brine volume 144 CALL iom_put( "snowpre" , sprecip * 86400. ) ! snow precipitation [m/day]145 127 CALL iom_put( "micesalt" , smt_i * zswi ) ! mean ice salinity 146 128 CALL iom_put( "snowvol" , vt_s * zswi ) ! snow volume 147 129 148 CALL iom_put( "sfxbog" , sfx_bog * rday ) ! salt flux from bottom growth149 CALL iom_put( "sfxbom" , sfx_bom * rday ) ! salt flux from bottom melting150 CALL iom_put( "sfxsum" , sfx_sum * rday ) ! salt flux from surface melting151 CALL iom_put( "sfxlam" , sfx_lam * rday ) ! salt flux from lateral melting152 CALL iom_put( "sfxsni" , sfx_sni * rday ) ! salt flux from snow ice formation153 CALL iom_put( "sfxopw" , sfx_opw * rday ) ! salt flux from open water formation154 CALL iom_put( "sfxdyn" , sfx_dyn * rday ) ! salt flux from ridging rafting155 CALL iom_put( "sfxres" , sfx_res * rday ) ! salt flux from corrections (resultant)156 CALL iom_put( "sfxbri" , sfx_bri * rday ) ! salt flux from brines157 CALL iom_put( "sfxsub" , sfx_sub * rday ) ! salt flux from sublimation158 CALL iom_put( "sfx" , sfx * rday ) ! total salt flux159 160 ztmp = rday / rhoic161 CALL iom_put( "vfxres" , wfx_res * ztmp ) ! daily prod./melting due to corrections162 CALL iom_put( "vfxopw" , wfx_opw * ztmp ) ! daily lateral thermodynamic ice production163 CALL iom_put( "vfxsni" , wfx_sni * ztmp ) ! daily snowice ice production164 CALL iom_put( "vfxbog" , wfx_bog * ztmp ) ! daily bottom thermodynamic ice production165 CALL iom_put( "vfxdyn" , wfx_dyn * ztmp ) ! daily dynamic ice production (rid/raft)166 CALL iom_put( "vfxsum" , wfx_sum * ztmp ) ! surface melt167 CALL iom_put( "vfxbom" , wfx_bom * ztmp ) ! bottom melt168 CALL iom_put( "vfxlam" , wfx_lam * ztmp ) ! lateral melt169 CALL iom_put( "vfxice" , wfx_ice * ztmp ) ! total ice growth/melt170 171 IF ( ln_pnd ) &172 CALL iom_put( "vfxpnd" , wfx_pnd * ztmp ) ! melt pond water flux173 174 IF ( iom_use( "vfxthin" ) ) THEN ! ice production for open water + thin ice (<20cm) => comparable to observations175 WHERE( htm_i(:,:) < 0.2 .AND. htm_i(:,:) > 0. ) ; z2d = wfx_bog176 ELSEWHERE ; z2d = 0._wp177 END WHERE178 CALL iom_put( "vfxthin", ( wfx_opw + z2d ) * ztmp )179 ENDIF180 181 ztmp = rday * r1_rhosn182 CALL iom_put( "vfxspr" , wfx_spr * ztmp ) ! precip (snow)183 CALL iom_put( "vfxsnw" , wfx_snw * ztmp ) ! total snw growth/melt184 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 ocean186 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 base203 CALL iom_put ('hfxdhc' , diag_heat(:,:) ) ! Heat content variation in snow and ice204 CALL iom_put ('hfxspr' , hfx_spr(:,:) ) ! Heat content of snow precip205 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 2016208 130 IF ( ln_pnd ) THEN 209 131 CALL iom_put( "iceamp" , at_ip * zswi ) ! melt pond total fraction 210 132 CALL iom_put( "icevmp" , vt_ip * zswi ) ! melt pond total volume per unit area 211 133 ENDIF 212 ! END MV MP 2016213 !!gm <<<<<<======= end214 134 215 135 !---------------------------------- … … 225 145 IF ( iom_use('brinevol_cat') ) CALL iom_put( "brinevol_cat", bv_i * 100. * zswi2 ) ! brine volume 226 146 227 ! MV MP 2016228 147 IF ( ln_pnd ) THEN 229 148 IF ( iom_use('iceamp_cat') ) CALL iom_put( "iceamp_cat" , a_ip * zswi2 ) ! melt pond frac for categories … … 232 151 IF ( iom_use('iceafp_cat') ) CALL iom_put( "iceafp_cat" , a_ip_frac * zswi2 ) ! melt pond frac for categories 233 152 ENDIF 234 ! END MV MP 2016235 153 236 154 !-------------------------------- … … 360 278 !! History : 4.0 ! 2013-06 (C. Rousset) 361 279 !!---------------------------------------------------------------------- 362 INTEGER, INTENT( in ) :: kt ! ocean time-step index )280 INTEGER, INTENT( in ) :: kt ! ocean time-step index 363 281 INTEGER, INTENT( in ) :: kid , kh_i 364 282 INTEGER :: nz_i, jl … … 476 394 CALL histwrite( kid, "sntemcat", kt, tm_su - rt0 , jpi*jpj*jpl, (/1/) ) 477 395 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 ) 482 398 ! 483 399 END SUBROUTINE ice_wri_state
Note: See TracChangeset
for help on using the changeset viewer.