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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

icethd.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90 @ 8626

Last change on this file since 8626 was 8626, checked in by clem, 7 years ago
rebin thickness after bdy in case itd is not the same between the input and simulation + prepare some changes in salinity
File size: 34.6 KB

Rev	Line
[8422]	1	MODULE icethd
	2	!!======================================================================
	3	!! * MODULE icethd *
[8534]	4	!! sea-ice : master routine for thermodynamics
[8422]	5	!!======================================================================
	6	!! History : LIM ! 2000-01 (M.A. Morales Maqueda, H. Goosse, T. Fichefet) LIM-1
	7	!! 2.0 ! 2002-07 (C. Ethe, G. Madec) LIM-2 (F90 rewriting)
	8	!! 3.0 ! 2005-11 (M. Vancoppenolle) LIM-3 : Multi-layer thermodynamics + salinity variations
	9	!! - ! 2007-04 (M. Vancoppenolle) add ice_thd_glohec, ice_thd_con_dh and ice_thd_con_dif
	10	!! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in wfx_snw
	11	!! 3.3 ! 2010-11 (G. Madec) corrected snow melting heat (due to factor betas)
	12	!! 4.0 ! 2011-02 (G. Madec) dynamical allocation
	13	!! - ! 2012-05 (C. Rousset) add penetration solar flux
	14	!!----------------------------------------------------------------------
	15	#if defined key_lim3
	16	!!----------------------------------------------------------------------
[8534]	17	!! 'key_lim3' ESIM sea-ice model
[8422]	18	!!----------------------------------------------------------------------
[8534]	19	!! ice_thd : thermodynamics of sea ice
	20	!! ice_thd_init : initialisation of sea-ice thermodynamics
[8422]	21	!!----------------------------------------------------------------------
	22	USE phycst ! physical constants
	23	USE dom_oce ! ocean space and time domain variables
[8534]	24	USE ice ! sea-ice: variables
[8486]	25	!!gm list trop longue ==>>> why not passage en argument d'appel ?
[8422]	26	USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, qns_tot, qsr_tot, sprecip, ln_cpl
	27	USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, &
[8512]	28	& fr1_i0, fr2_i0
[8534]	29	USE ice1D ! sea-ice: thermodynamics variables
	30	USE icethd_zdf ! sea-ice: vertical heat diffusion
	31	USE icethd_dh ! sea-ice: ice-snow growth and melt
	32	USE icethd_da ! sea-ice: lateral melting
	33	USE icethd_sal ! sea-ice: salinity
	34	USE icethd_ent ! sea-ice: enthalpy redistribution
	35	USE icethd_do ! sea-ice: growth in open water
[8598]	36	USE icethd_pnd ! sea-ice: melt ponds
[8534]	37	USE iceitd ! sea-ice: remapping thickness distribution
	38	USE icetab ! sea-ice: 1D <==> 2D transformation
	39	USE icevar ! sea-ice: operations
	40	USE icectl ! sea-ice: control print
[8422]	41	!
	42	USE in_out_manager ! I/O manager
	43	USE lib_mpp ! MPP library
[8534]	44	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
	45	USE lbclnk ! lateral boundary conditions (or mpp links)
[8422]	46	USE timing ! Timing
	47
	48	IMPLICIT NONE
	49	PRIVATE
	50
	51	PUBLIC ice_thd ! called by limstp module
	52	PUBLIC ice_thd_init ! called by ice_init
	53
[8531]	54	!! namelist (namthd)
	55	LOGICAL :: ln_icedH ! activate ice thickness change from growing/melting (T) or not (F)
	56	LOGICAL :: ln_icedA ! activate lateral melting param. (T) or not (F)
	57	LOGICAL :: ln_icedO ! activate ice growth in open-water (T) or not (F)
	58	LOGICAL :: ln_icedS ! activate gravity drainage and flushing (T) or not (F)
	59
[8422]	60	!! * Substitutions
	61	# include "vectopt_loop_substitute.h90"
	62	!!----------------------------------------------------------------------
[8486]	63	!! NEMO/ICE 4.0 , NEMO Consortium (2017)
[8422]	64	!! $Id: icethd.F90 8420 2017-08-08 12:18:46Z clem $
	65	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	66	!!----------------------------------------------------------------------
	67	CONTAINS
	68
	69	SUBROUTINE ice_thd( kt )
	70	!!-------------------------------------------------------------------
	71	!! * ROUTINE ice_thd *
	72	!!
	73	!! ** Purpose : This routine manages ice thermodynamics
	74	!!
	75	!! ** Action : - Initialisation of some variables
	76	!! - Some preliminary computation (oceanic heat flux
	77	!! at the ice base, snow acc.,heat budget of the leads)
	78	!! - selection of the icy points and put them in an array
[8531]	79	!! - call ice_thd_zdf for vertical heat diffusion
[8422]	80	!! - call ice_thd_dh for vertical ice growth and melt
	81	!! - call ice_thd_ent for enthalpy remapping
	82	!! - call ice_thd_sal for ice desalination
	83	!! - call ice_thd_temp to retrieve temperature from ice enthalpy
[8531]	84	!! - call ice_thd_lam for extra lateral ice melt if monocat
	85	!! - call ice_thd_da for lateral ice melt
[8422]	86	!! - back to the geographic grid
[8534]	87	!! - call ice_thd_rem for remapping thickness distribution
	88	!! - call ice_thd_do for ice growth in leads
[8598]	89	!!-------------------------------------------------------------------
[8422]	90	INTEGER, INTENT(in) :: kt ! number of iteration
	91	!
	92	INTEGER :: ji, jj, jk, jl ! dummy loop indices
	93	REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg
	94	REAL(wp), PARAMETER :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04)
	95	REAL(wp), PARAMETER :: zch = 0.0057_wp ! heat transfer coefficient
	96	REAL(wp), DIMENSION(jpi,jpj) :: zu_io, zv_io, zfric ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)
	97	!
	98	!!-------------------------------------------------------------------
[8517]	99	! controls
	100	IF( nn_timing == 1 ) CALL timing_start('icethd') ! timing
	101	IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
[8422]	102
	103	IF( kt == nit000 .AND. lwp ) THEN
[8426]	104	WRITE(numout,*)
[8512]	105	WRITE(numout,*) 'ice_thd: sea-ice thermodynamics'
	106	WRITE(numout,*) '~~~~~~~'
[8422]	107	ENDIF
	108
[8424]	109	CALL ice_var_glo2eqv
[8422]	110
	111	!---------------------------------------------!
	112	! computation of friction velocity at T points
	113	!---------------------------------------------!
[8512]	114	IF( ln_icedyn ) THEN
[8422]	115	zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)
	116	zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
	117	DO jj = 2, jpjm1
	118	DO ji = fs_2, fs_jpim1
	119	zfric(ji,jj) = rn_cio * ( 0.5_wp * &
	120	& ( zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj) &
	121	& + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) ) * tmask(ji,jj,1)
	122	END DO
	123	END DO
	124	ELSE ! if no ice dynamics => transmit directly the atmospheric stress to the ocean
	125	DO jj = 2, jpjm1
	126	DO ji = fs_2, fs_jpim1
	127	zfric(ji,jj) = r1_rau0 * SQRT( 0.5_wp * &
	128	& ( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) &
	129	& + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1)
	130	END DO
	131	END DO
	132	ENDIF
	133	CALL lbc_lnk( zfric, 'T', 1. )
	134	!
	135	ftr_ice(:,:,:) = 0._wp ! initialization (part of solar radiation transmitted through the ice)
	136
	137	!--------------------------------------------------------------------!
	138	! Partial computation of forcing for the thermodynamic sea ice model
	139	!--------------------------------------------------------------------!
	140	DO jj = 1, jpj
	141	DO ji = 1, jpi
	142	rswitch = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice
	143	!
	144	! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
	145	! ! practically no "direct lateral ablation"
	146	!
	147	! ! net downward heat flux from the ice to the ocean, expressed as a function of ocean
	148	! ! temperature and turbulent mixing (McPhee, 1992)
	149	!
	150	! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- !
	151	zqld = tmask(ji,jj,1) * rdt_ice * &
	152	& ( ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) * frq_m(ji,jj) + &
	153	& ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) )
	154
	155	! --- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- !
	156	! includes supercooling potential energy (>0) or "above-freezing" energy (<0)
	157	zqfr = tmask(ji,jj,1) * rau0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) )
	158
	159	! --- Above-freezing sensible heat content (J/m2 grid)
	160	zqfr_neg = tmask(ji,jj,1) * rau0 * rcp * e3t_m(ji,jj) * MIN( ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ), 0._wp )
	161
	162	! --- Sensible ocean-to-ice heat flux (W/m2)
	163	zfric_u = MAX( SQRT( zfric(ji,jj) ), zfric_umin )
	164	fhtur(ji,jj) = rswitch * rau0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2
	165
	166	fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - zqfr_neg * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )
	167	! upper bound for fhtur: the heat retrieved from the ocean must be smaller than the heat necessary to reach
	168	! the freezing point, so that we do not have SST < T_freeze
	169	! This implies: - ( fhtur(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0
	170
	171	!-- Energy Budget of the leads (J.m-2), source of lateral accretion. Must be < 0 to form ice
	172	qlead(ji,jj) = MIN( 0._wp , zqld - ( fhtur(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr )
	173
	174	! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting
	175	IF( zqld > 0._wp ) THEN
	176	fhld (ji,jj) = rswitch * zqld * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90
	177	qlead(ji,jj) = 0._wp
	178	ELSE
	179	fhld (ji,jj) = 0._wp
	180	ENDIF
	181	!
	182	! Net heat flux on top of the ice-ocean [W.m-2]
	183	! ---------------------------------------------
	184	hfx_in(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj)
	185	END DO
	186	END DO
	187
	188	! In case we bypass open-water ice formation
[8514]	189	IF( .NOT. ln_icedO ) qlead(:,:) = 0._wp
[8422]	190	! In case we bypass growing/melting from top and bottom: we suppose ice is impermeable => ocean is isolated from atmosphere
[8514]	191	IF( .NOT. ln_icedH ) THEN
[8483]	192	hfx_in(:,:) = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
	193	fhtur (:,:) = 0._wp
	194	fhld (:,:) = 0._wp
	195	ENDIF
[8422]	196
	197	! ---------------------------------------------------------------------
	198	! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2]
	199	! ---------------------------------------------------------------------
	200	! First step here : non solar + precip - qlead - qturb
	201	! Second step in icethd_dh : heat remaining if total melt (zq_rema)
	202	! Third step in iceupdate.F90 : heat from ice-ocean mass exchange (zf_mass) + solar
[8498]	203	hfx_out(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:) & ! Non solar heat flux received by the ocean
	204	& - qlead(:,:) * r1_rdtice & ! heat flux taken from the ocean where there is open water ice formation
	205	& - at_i (:,:) * fhtur(:,:) & ! heat flux taken by turbulence
	206	& - at_i (:,:) * fhld(:,:) ! heat flux taken during bottom growth/melt
	207	! (fhld should be 0 while bott growth)
[8422]	208	!-------------------------------------------------------------------------------------------!
	209	! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories
	210	!-------------------------------------------------------------------------------------------!
	211	DO jl = 1, jpl
	212
	213	! select ice covered grid points
[8565]	214	npti = 0 ; nptidx(:) = 0
[8422]	215	DO jj = 1, jpj
	216	DO ji = 1, jpi
	217	IF ( a_i(ji,jj,jl) > epsi10 ) THEN
[8565]	218	npti = npti + 1
	219	nptidx(npti) = (jj - 1) * jpi + ji
[8422]	220	ENDIF
	221	END DO
	222	END DO
	223
[8565]	224	IF( lk_mpp ) CALL mpp_ini_ice( npti , numout )
[8422]	225
[8565]	226	IF( npti > 0 ) THEN ! If there is no ice, do nothing.
[8422]	227	!
	228	CALL ice_thd_1d2d( jl, 1 ) ! --- Move to 1D arrays --- !
	229	! ! --- & Change units of e_i, e_s from J/m2 to J/m3 --- !
	230	!
[8565]	231	s_i_new (1:npti) = 0._wp ; dh_s_tot (1:npti) = 0._wp ! --- some init --- ! (important to have them here)
	232	dh_i_surf (1:npti) = 0._wp ; dh_i_bott(1:npti) = 0._wp
[8623]	233	dh_snowice(1:npti) = 0._wp ; dh_i_sub (1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp
[8422]	234	!
[8514]	235	IF( ln_icedH ) THEN ! --- growing/melting --- !
[8531]	236	CALL ice_thd_zdf ! Ice/Snow Temperature profile
[8486]	237	CALL ice_thd_dh ! Ice/Snow thickness
[8623]	238	CALL ice_thd_pnd ! Melt ponds formation
[8565]	239	CALL ice_thd_ent( e_i_1d(1:npti,:) ) ! Ice enthalpy remapping
[8486]	240	ENDIF
[8422]	241	!
[8531]	242	CALL ice_thd_sal( ln_icedS ) ! --- Ice salinity --- !
[8422]	243	!
	244	CALL ice_thd_temp ! --- temperature update --- !
	245	!
[8486]	246	!!gm please create a new logical (l_thd_lam or a better explicit name) set one for all in icestp.F90 module
[8514]	247	!!gm l_thd_lam = ln_icedH .AND. ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 )
[8486]	248	!!gm by the way, the different options associated with nn_monocat =1 to 4 are quite impossible to identify
	249	!!gm more comment to add when ready the namelist, with an explicit print in the ocean.output
[8514]	250	IF( ln_icedH ) THEN
[8422]	251	IF ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 ) THEN
	252	CALL ice_thd_lam ! --- extra lateral melting if monocat --- !
	253	END IF
	254	END IF
	255	!
[8514]	256	IF( ln_icedA ) CALL ice_thd_da ! --- lateral melting --- !
[8422]	257	!
	258	CALL ice_thd_1d2d( jl, 2 ) ! --- Change units of e_i, e_s from J/m3 to J/m2 --- !
	259	! ! --- & Move to 2D arrays --- !
	260	!
	261	IF( lk_mpp ) CALL mpp_comm_free( ncomm_ice ) !RB necessary ??
	262	ENDIF
	263	!
	264	END DO
	265	! update ice age (in case a_i changed, i.e. becomes 0 or lateral melting)
	266	oa_i(:,:,:) = o_i(:,:,:) * a_i(:,:,:)
	267
[8517]	268	IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
[8422]	269	!
[8486]	270	CALL ice_var_zapsmall ! --- remove very small ice concentration (<1e-10) --- !
	271	! ! & make sure at_i=SUM(a_i) & ato_i=1 where at_i=0
[8422]	272	!
[8486]	273	IF( jpl > 1 ) CALL ice_itd_rem( kt ) ! --- Transport ice between thickness categories --- !
[8422]	274	!
[8531]	275	IF( ln_icedO ) CALL ice_thd_do ! --- frazil ice growing in leads --- !
[8422]	276	!
[8517]	277	! controls
	278	IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ') ! prints
	279	IF( ln_ctl ) CALL ice_prt3D ('icethd') ! prints
	280	IF( nn_timing == 1 ) CALL timing_stop('icethd') ! timing
[8422]	281	!
	282	END SUBROUTINE ice_thd
	283
	284
	285	SUBROUTINE ice_thd_temp
	286	!!-----------------------------------------------------------------------
	287	!! * ROUTINE ice_thd_temp *
	288	!!
	289	!! ** Purpose : Computes sea ice temperature (Kelvin) from enthalpy
	290	!!
	291	!! ** Method : Formula (Bitz and Lipscomb, 1999)
	292	!!-------------------------------------------------------------------
	293	INTEGER :: ji, jk ! dummy loop indices
[8522]	294	REAL(wp) :: ztmelts, zbbb, zccc ! local scalar
[8422]	295	!!-------------------------------------------------------------------
	296	! Recover ice temperature
	297	DO jk = 1, nlay_i
[8565]	298	DO ji = 1, npti
[8564]	299	ztmelts = -tmut * sz_i_1d(ji,jk)
[8422]	300	! Conversion q(S,T) -> T (second order equation)
[8498]	301	zbbb = ( rcp - cpic ) * ztmelts + e_i_1d(ji,jk) * r1_rhoic - lfus
	302	zccc = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts, 0._wp ) )
[8522]	303	t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_cpic
[8422]	304
	305	! mask temperature
[8563]	306	rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) )
[8498]	307	t_i_1d(ji,jk) = rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0
[8422]	308	END DO
	309	END DO
	310	!
	311	END SUBROUTINE ice_thd_temp
	312
	313
	314	SUBROUTINE ice_thd_lam
	315	!!-----------------------------------------------------------------------
	316	!! * ROUTINE ice_thd_lam *
	317	!!
	318	!! ** Purpose : Lateral melting in case monocategory
	319	!! ( dA = A/2h dh )
	320	!!-----------------------------------------------------------------------
	321	INTEGER :: ji ! dummy loop indices
	322	REAL(wp) :: zhi_bef ! ice thickness before thermo
	323	REAL(wp) :: zdh_mel, zda_mel ! net melting
	324	REAL(wp) :: zvi, zvs ! ice/snow volumes
	325	!!-----------------------------------------------------------------------
	326	!
[8565]	327	DO ji = 1, npti
[8422]	328	zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) + dh_i_sub(ji) )
	329	IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp ) THEN
[8563]	330	zvi = a_i_1d(ji) * h_i_1d(ji)
	331	zvs = a_i_1d(ji) * h_s_1d(ji)
[8422]	332	! lateral melting = concentration change
[8563]	333	zhi_bef = h_i_1d(ji) - zdh_mel
[8422]	334	rswitch = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) )
	335	zda_mel = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) )
	336	a_i_1d(ji) = MAX( epsi20, a_i_1d(ji) + zda_mel )
	337	! adjust thickness
[8563]	338	h_i_1d(ji) = zvi / a_i_1d(ji)
	339	h_s_1d(ji) = zvs / a_i_1d(ji)
[8422]	340	! retrieve total concentration
	341	at_i_1d(ji) = a_i_1d(ji)
	342	END IF
	343	END DO
	344	!
	345	END SUBROUTINE ice_thd_lam
	346
	347
[8486]	348	SUBROUTINE ice_thd_1d2d( kl, kn )
[8422]	349	!!-----------------------------------------------------------------------
	350	!! * ROUTINE ice_thd_1d2d *
	351	!!
	352	!! ** Purpose : move arrays from 1d to 2d and the reverse
	353	!!-----------------------------------------------------------------------
[8486]	354	INTEGER, INTENT(in) :: kl ! index of the ice category
	355	INTEGER, INTENT(in) :: kn ! 1= from 2D to 1D ; 2= from 1D to 2D
[8422]	356	!
[8486]	357	INTEGER :: jk ! dummy loop indices
[8422]	358	!!-----------------------------------------------------------------------
	359	!
	360	SELECT CASE( kn )
[8486]	361	! !---------------------!
	362	CASE( 1 ) !== from 2D to 1D ==!
	363	! !---------------------!
[8565]	364	CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i )
	365	CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl) )
[8623]	366	CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,kl) )
	367	CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti), h_s (:,:,kl) )
[8565]	368	CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl) )
[8623]	369	CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti), s_i (:,:,kl) )
[8422]	370	DO jk = 1, nlay_s
[8623]	371	CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl) )
	372	CALL tab_2d_1d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl) )
[8422]	373	END DO
	374	DO jk = 1, nlay_i
[8623]	375	CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,kl) )
	376	CALL tab_2d_1d( npti, nptidx(1:npti), e_i_1d (1:npti,jk), e_i (:,:,jk,kl) )
	377	CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl) )
[8422]	378	END DO
[8623]	379	CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,kl) )
	380	CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,kl) )
	381	CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) )
[8422]	382	!
[8565]	383	CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d(1:npti), qprec_ice )
	384	CALL tab_2d_1d( npti, nptidx(1:npti), qsr_ice_1d (1:npti), qsr_ice (:,:,kl) )
	385	CALL tab_2d_1d( npti, nptidx(1:npti), fr1_i0_1d (1:npti), fr1_i0 )
	386	CALL tab_2d_1d( npti, nptidx(1:npti), fr2_i0_1d (1:npti), fr2_i0 )
	387	CALL tab_2d_1d( npti, nptidx(1:npti), qns_ice_1d (1:npti), qns_ice (:,:,kl) )
	388	CALL tab_2d_1d( npti, nptidx(1:npti), ftr_ice_1d (1:npti), ftr_ice (:,:,kl) )
	389	CALL tab_2d_1d( npti, nptidx(1:npti), evap_ice_1d (1:npti), evap_ice(:,:,kl) )
	390	CALL tab_2d_1d( npti, nptidx(1:npti), dqns_ice_1d (1:npti), dqns_ice(:,:,kl) )
	391	CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d (1:npti), t_bo )
	392	CALL tab_2d_1d( npti, nptidx(1:npti), sprecip_1d (1:npti), sprecip )
	393	CALL tab_2d_1d( npti, nptidx(1:npti), fhtur_1d (1:npti), fhtur )
	394	CALL tab_2d_1d( npti, nptidx(1:npti), fhld_1d (1:npti), fhld )
[8422]	395	!
[8565]	396	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni )
	397	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum )
	398	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sub_1d (1:npti), wfx_sub )
	399	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sub_1d(1:npti), wfx_snw_sub )
	400	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_ice_sub_1d(1:npti), wfx_ice_sub )
	401	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_err_sub_1d(1:npti), wfx_err_sub )
[8422]	402	!
[8565]	403	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_bog_1d (1:npti), wfx_bog )
	404	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_bom_1d (1:npti), wfx_bom )
	405	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum )
	406	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sni_1d (1:npti), wfx_sni )
	407	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_res_1d (1:npti), wfx_res )
	408	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr )
	409	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam )
[8623]	410	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd )
[8422]	411	!
[8565]	412	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog )
	413	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bom_1d (1:npti), sfx_bom )
	414	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sum_1d (1:npti), sfx_sum )
	415	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sni_1d (1:npti), sfx_sni )
	416	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri )
	417	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_res_1d (1:npti), sfx_res )
	418	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sub_1d (1:npti), sfx_sub )
	419	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_lam_1d (1:npti), sfx_lam )
[8422]	420	!
[8565]	421	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d (1:npti), hfx_thd )
	422	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_spr_1d (1:npti), hfx_spr )
	423	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_sum_1d (1:npti), hfx_sum )
	424	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_bom_1d (1:npti), hfx_bom )
	425	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_bog_1d (1:npti), hfx_bog )
	426	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dif_1d (1:npti), hfx_dif )
	427	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_opw_1d (1:npti), hfx_opw )
	428	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_snw_1d (1:npti), hfx_snw )
	429	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_sub_1d (1:npti), hfx_sub )
	430	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_res_1d (1:npti), hfx_res )
	431	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif )
	432	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_rem_1d(1:npti), hfx_err_rem )
	433	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_out_1d (1:npti), hfx_out )
[8422]	434	!
	435	! SIMIP diagnostics
[8565]	436	CALL tab_2d_1d( npti, nptidx(1:npti), diag_fc_bo_1d(1:npti), diag_fc_bo )
	437	CALL tab_2d_1d( npti, nptidx(1:npti), diag_fc_su_1d(1:npti), diag_fc_su )
[8422]	438	! ocean surface fields
[8565]	439	CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m )
	440	CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m )
[8422]	441
	442	! --- Change units of e_i, e_s from J/m2 to J/m3 --- !
	443	DO jk = 1, nlay_i
[8565]	444	WHERE( h_i_1d(1:npti)>0._wp ) e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) / (h_i_1d(1:npti) * a_i_1d(1:npti)) * nlay_i
[8486]	445	END DO
[8422]	446	DO jk = 1, nlay_s
[8565]	447	WHERE( h_s_1d(1:npti)>0._wp ) e_s_1d(1:npti,jk) = e_s_1d(1:npti,jk) / (h_s_1d(1:npti) * a_i_1d(1:npti)) * nlay_s
[8486]	448	END DO
[8422]	449	!
[8486]	450	! !---------------------!
	451	CASE( 2 ) !== from 1D to 2D ==!
	452	! !---------------------!
[8422]	453	! --- Change units of e_i, e_s from J/m3 to J/m2 --- !
	454	DO jk = 1, nlay_i
[8565]	455	e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) * h_i_1d(1:npti) * a_i_1d(1:npti) * r1_nlay_i
[8486]	456	END DO
[8422]	457	DO jk = 1, nlay_s
[8565]	458	e_s_1d(1:npti,jk) = e_s_1d(1:npti,jk) * h_s_1d(1:npti) * a_i_1d(1:npti) * r1_nlay_s
[8486]	459	END DO
[8422]	460	!
[8506]	461	! Change thickness to volume (replaces routine ice_var_eqv2glo)
[8623]	462	v_i_1d (1:npti) = h_i_1d (1:npti) * a_i_1d (1:npti)
	463	v_s_1d (1:npti) = h_s_1d (1:npti) * a_i_1d (1:npti)
	464	sv_i_1d(1:npti) = s_i_1d (1:npti) * v_i_1d (1:npti)
	465	v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti)
[8422]	466
[8565]	467	CALL tab_1d_2d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i )
	468	CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl) )
[8623]	469	CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,kl) )
	470	CALL tab_1d_2d( npti, nptidx(1:npti), h_s_1d (1:npti), h_s (:,:,kl) )
[8565]	471	CALL tab_1d_2d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl) )
[8623]	472	CALL tab_1d_2d( npti, nptidx(1:npti), s_i_1d (1:npti), s_i (:,:,kl) )
[8422]	473	DO jk = 1, nlay_s
[8623]	474	CALL tab_1d_2d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl) )
	475	CALL tab_1d_2d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl) )
[8422]	476	END DO
	477	DO jk = 1, nlay_i
[8623]	478	CALL tab_1d_2d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,kl) )
	479	CALL tab_1d_2d( npti, nptidx(1:npti), e_i_1d (1:npti,jk), e_i (:,:,jk,kl) )
	480	CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl) )
[8422]	481	END DO
[8623]	482	CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,kl) )
	483	CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,kl) )
	484	CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) )
[8422]	485	!
[8565]	486	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni )
	487	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum )
	488	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sub_1d (1:npti), wfx_sub )
	489	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sub_1d(1:npti), wfx_snw_sub )
	490	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_ice_sub_1d(1:npti), wfx_ice_sub )
	491	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_err_sub_1d(1:npti), wfx_err_sub )
[8422]	492	!
[8565]	493	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_bog_1d (1:npti), wfx_bog )
	494	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_bom_1d (1:npti), wfx_bom )
	495	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum )
	496	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sni_1d (1:npti), wfx_sni )
	497	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_res_1d (1:npti), wfx_res )
	498	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr )
	499	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam )
[8623]	500	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd )
[8422]	501	!
[8565]	502	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog )
	503	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bom_1d (1:npti), sfx_bom )
	504	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sum_1d (1:npti), sfx_sum )
	505	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sni_1d (1:npti), sfx_sni )
	506	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri )
	507	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_res_1d (1:npti), sfx_res )
	508	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sub_1d (1:npti), sfx_sub )
	509	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_lam_1d (1:npti), sfx_lam )
[8422]	510	!
[8565]	511	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_thd_1d (1:npti), hfx_thd )
	512	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_spr_1d (1:npti), hfx_spr )
	513	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_sum_1d (1:npti), hfx_sum )
	514	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_bom_1d (1:npti), hfx_bom )
	515	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_bog_1d (1:npti), hfx_bog )
	516	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dif_1d (1:npti), hfx_dif )
	517	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_opw_1d (1:npti), hfx_opw )
	518	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_snw_1d (1:npti), hfx_snw )
	519	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_sub_1d (1:npti), hfx_sub )
	520	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_res_1d (1:npti), hfx_res )
	521	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif )
	522	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_rem_1d(1:npti), hfx_err_rem )
	523	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_out_1d (1:npti), hfx_out )
[8422]	524	!
[8565]	525	CALL tab_1d_2d( npti, nptidx(1:npti), qns_ice_1d(1:npti), qns_ice(:,:,kl) )
	526	CALL tab_1d_2d( npti, nptidx(1:npti), ftr_ice_1d(1:npti), ftr_ice(:,:,kl) )
[8422]	527	!
	528	! SIMIP diagnostics
[8565]	529	CALL tab_1d_2d( npti, nptidx(1:npti), t_si_1d (1:npti), t_si(:,:,kl) )
	530	CALL tab_1d_2d( npti, nptidx(1:npti), diag_fc_bo_1d(1:npti), diag_fc_bo )
	531	CALL tab_1d_2d( npti, nptidx(1:npti), diag_fc_su_1d(1:npti), diag_fc_su )
[8422]	532	! extensive variables
[8565]	533	CALL tab_1d_2d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i (:,:,kl) )
	534	CALL tab_1d_2d( npti, nptidx(1:npti), v_s_1d (1:npti), v_s (:,:,kl) )
	535	CALL tab_1d_2d( npti, nptidx(1:npti), sv_i_1d(1:npti), sv_i(:,:,kl) )
[8623]	536	CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,kl) )
[8486]	537	!
[8422]	538	END SELECT
	539	!
	540	END SUBROUTINE ice_thd_1d2d
	541
	542
	543	SUBROUTINE ice_thd_init
[8598]	544	!!-------------------------------------------------------------------
[8422]	545	!! * ROUTINE ice_thd_init *
	546	!!
[8531]	547	!! ** Purpose : Physical constants and parameters associated with
	548	!! ice thermodynamics
[8422]	549	!!
[8531]	550	!! ** Method : Read the namthd namelist and check the parameters
	551	!! called at the first timestep (nit000)
[8422]	552	!!
[8531]	553	!! ** input : Namelist namthd
[8422]	554	!!-------------------------------------------------------------------
[8486]	555	INTEGER :: ios ! Local integer output status for namelist read
	556	!!
[8531]	557	NAMELIST/namthd/ ln_icedH, ln_icedA, ln_icedO, ln_icedS
[8422]	558	!!-------------------------------------------------------------------
	559	!
[8531]	560	REWIND( numnam_ice_ref ) ! Namelist namthd in reference namelist : Ice thermodynamics
	561	READ ( numnam_ice_ref, namthd, IOSTAT = ios, ERR = 901)
	562	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd in reference namelist', lwp )
[8422]	563
[8531]	564	REWIND( numnam_ice_cfg ) ! Namelist namthd in configuration namelist : Ice thermodynamics
	565	READ ( numnam_ice_cfg, namthd, IOSTAT = ios, ERR = 902 )
	566	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd in configuration namelist', lwp )
	567	IF(lwm) WRITE ( numoni, namthd )
[8422]	568	!
	569	!
	570	IF(lwp) THEN ! control print
[8531]	571	WRITE(numout,*) 'ice_thd_init: Ice Thermodynamics'
	572	WRITE(numout,*) '~~~~~~~~~~~~'
	573	WRITE(numout,*) ' Namelist namthd:'
	574	WRITE(numout,*) ' activate ice thick change from top/bot (T) or not (F) ln_icedH = ', ln_icedH
	575	WRITE(numout,*) ' activate lateral melting (T) or not (F) ln_icedA = ', ln_icedA
	576	WRITE(numout,*) ' activate ice growth in open-water (T) or not (F) ln_icedO = ', ln_icedO
	577	WRITE(numout,*) ' activate gravity drainage and flushing (T) or not (F) ln_icedS = ', ln_icedS
	578	ENDIF
[8422]	579	!
[8531]	580	CALL ice_thd_zdf_init ! set ice heat diffusion parameters
	581	IF( ln_icedA ) CALL ice_thd_da_init ! set ice lateral melting parameters
	582	IF( ln_icedO ) CALL ice_thd_do_init ! set ice growth in open water parameters
	583	CALL ice_thd_sal_init ! set ice salinity parameters
[8598]	584	CALL ice_thd_pnd_init ! set melt ponds parameters
[8531]	585	!
[8422]	586	END SUBROUTINE ice_thd_init
	587
	588	#else
	589	!!----------------------------------------------------------------------
[8534]	590	!! Default option Dummy module NO ESIM sea-ice model
[8422]	591	!!----------------------------------------------------------------------
	592	#endif
	593
	594	!!======================================================================
	595	END MODULE icethd

Note: See TracBrowser for help on using the repository browser.

Download in other formats: