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.

limthd.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90 @ 4788

Last change on this file since 4788 was 4765, checked in by rblod, 10 years ago
Compilation issue, see ticket #1379
Property svn:keywords set to `Id`
File size: 33.6 KB

Rev	Line
[825]	1	MODULE limthd
	2	!!======================================================================
	3	!! * MODULE limthd *
[1572]	4	!! LIM-3 : ice thermodynamic
[825]	5	!!======================================================================
[1572]	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
[2528]	9	!! - ! 2007-04 (M. Vancoppenolle) add lim_thd_glohec, lim_thd_con_dh and lim_thd_con_dif
[4688]	10	!! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in wfx_snw
[2528]	11	!! 3.3 ! 2010-11 (G. Madec) corrected snow melting heat (due to factor betas)
[2715]	12	!! 4.0 ! 2011-02 (G. Madec) dynamical allocation
[4161]	13	!! - ! 2012-05 (C. Rousset) add penetration solar flux
[1572]	14	!!----------------------------------------------------------------------
[825]	15	#if defined key_lim3
	16	!!----------------------------------------------------------------------
[834]	17	!! 'key_lim3' LIM3 sea-ice model
[825]	18	!!----------------------------------------------------------------------
[3625]	19	!! lim_thd : thermodynamic of sea ice
	20	!! lim_thd_init : initialisation of sea-ice thermodynamic
[825]	21	!!----------------------------------------------------------------------
[3625]	22	USE phycst ! physical constants
	23	USE dom_oce ! ocean space and time domain variables
[4205]	24	USE oce , ONLY : iatte, oatte
[3625]	25	USE ice ! LIM: sea-ice variables
	26	USE par_ice ! LIM: sea-ice parameters
	27	USE sbc_oce ! Surface boundary condition: ocean fields
	28	USE sbc_ice ! Surface boundary condition: ice fields
	29	USE thd_ice ! LIM thermodynamic sea-ice variables
	30	USE dom_ice ! LIM sea-ice domain
	31	USE domvvl ! domain: variable volume level
	32	USE limthd_dif ! LIM: thermodynamics, vertical diffusion
	33	USE limthd_dh ! LIM: thermodynamics, ice and snow thickness variation
	34	USE limthd_sal ! LIM: thermodynamics, ice salinity
	35	USE limthd_ent ! LIM: thermodynamics, ice enthalpy redistribution
	36	USE limtab ! LIM: 1D <==> 2D transformation
	37	USE limvar ! LIM: sea-ice variables
	38	USE lbclnk ! lateral boundary condition - MPP links
	39	USE lib_mpp ! MPP library
	40	USE wrk_nemo ! work arrays
	41	USE in_out_manager ! I/O manager
	42	USE prtctl ! Print control
	43	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[4161]	44	USE timing ! Timing
[4688]	45	USE cpl_oasis3, ONLY : lk_cpl
	46	USE limcons ! conservation tests
[825]	47
	48	IMPLICIT NONE
	49	PRIVATE
	50
[2528]	51	PUBLIC lim_thd ! called by limstp module
	52	PUBLIC lim_thd_init ! called by iceini module
[825]	53
[4333]	54	REAL(wp) :: epsi10 = 1.e-10_wp !
[825]	55
	56	!! * Substitutions
	57	# include "domzgr_substitute.h90"
	58	# include "vectopt_loop_substitute.h90"
	59	!!----------------------------------------------------------------------
[2528]	60	!! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
[1156]	61	!! $Id$
[2528]	62	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	63	!!----------------------------------------------------------------------
	64	CONTAINS
	65
[921]	66	SUBROUTINE lim_thd( kt )
[825]	67	!!-------------------------------------------------------------------
	68	!! * ROUTINE lim_thd *
	69	!!
	70	!! ** Purpose : This routine manages the ice thermodynamic.
	71	!!
	72	!! ** Action : - Initialisation of some variables
	73	!! - Some preliminary computation (oceanic heat flux
	74	!! at the ice base, snow acc.,heat budget of the leads)
	75	!! - selection of the icy points and put them in an array
	76	!! - call lim_vert_ther for vert ice thermodynamic
	77	!! - back to the geographic grid
	78	!! - selection of points for lateral accretion
	79	!! - call lim_lat_acc for the ice accretion
	80	!! - back to the geographic grid
	81	!!
[1572]	82	!! ** References : H. Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90
	83	!!---------------------------------------------------------------------
	84	INTEGER, INTENT(in) :: kt ! number of iteration
[825]	85	!!
[4688]	86	INTEGER :: ji, jj, jk, jl ! dummy loop indices
	87	INTEGER :: nbpb ! nb of icy pts for thermo. cal.
	88	INTEGER :: ii, ij ! temporary dummy loop index
	89	REAL(wp) :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04)
	90	REAL(wp) :: zch = 0.0057_wp ! heat transfer coefficient
	91	REAL(wp) :: zinda, zindb, zareamin
	92	REAL(wp) :: zfric_u, zqld, zqfr
	93	!
	94	REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
[825]	95	!!-------------------------------------------------------------------
[4161]	96	IF( nn_timing == 1 ) CALL timing_start('limthd')
[2715]	97
[4688]	98	! conservation test
	99	IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[4161]	100
[921]	101	!------------------------------------------------------------------------------!
	102	! 1) Initialization of diagnostic variables !
	103	!------------------------------------------------------------------------------!
[825]	104
	105	!--------------------
	106	! 1.2) Heat content
	107	!--------------------
[1572]	108	! Change the units of heat content; from global units to J.m3
[825]	109	DO jl = 1, jpl
[921]	110	DO jk = 1, nlay_i
	111	DO jj = 1, jpj
	112	DO ji = 1, jpi
	113	!0 if no ice and 1 if yes
[4333]	114	zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) )
[4688]	115	!Energy of melting q(S,T) [J.m-3]
	116	e_i(ji,jj,jk,jl) = zindb * e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i )
	117	!convert units ! very important that this line is here
	118	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac
[921]	119	END DO
[825]	120	END DO
[921]	121	END DO
	122	DO jk = 1, nlay_s
	123	DO jj = 1, jpj
	124	DO ji = 1, jpi
	125	!0 if no ice and 1 if yes
[4333]	126	zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi10 ) )
[4688]	127	!Energy of melting q(S,T) [J.m-3]
	128	e_s(ji,jj,jk,jl) = zindb * e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL( nlay_s )
[921]	129	!convert units ! very important that this line is here
[4688]	130	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac
[921]	131	END DO
[825]	132	END DO
[921]	133	END DO
[825]	134	END DO
	135
[921]	136	! 2) Partial computation of forcing for the thermodynamic sea ice model. !
	137	!-----------------------------------------------------------------------------!
[825]	138
[921]	139	!CDIR NOVERRCHK
	140	DO jj = 1, jpj
	141	!CDIR NOVERRCHK
	142	DO ji = 1, jpi
[4688]	143	zinda = tms(ji,jj) * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice
[2528]	144	!
[921]	145	! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
	146	! ! practically no "direct lateral ablation"
	147	!
	148	! ! net downward heat flux from the ice to the ocean, expressed as a function of ocean
	149	! ! temperature and turbulent mixing (McPhee, 1992)
[4688]	150	!
	151	! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- !
	152	zqld = tms(ji,jj) * rdt_ice * &
	153	& ( pfrld(ji,jj) * ( qsr(ji,jj) * oatte(ji,jj) & ! solar heat + clem modif
	154	& + qns(ji,jj) ) & ! non solar heat
	155	! latent heat of precip (note that precip is included in qns but not in qns_ice)
[4765]	156	& + ( pfrld(ji,jj)*betas - pfrld(ji,jj) ) sprecip(ji,jj) &
	157	& * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus ) &
	158	& + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) &
	159	& * rcp * ( tatm_ice(ji,jj) - rtt ) )
[825]	160
[4688]	161	!-- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- !
	162	zqfr = tms(ji,jj) * rau0 * rcp * fse3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) )
	163
	164	!-- Energy Budget of the leads (J.m-2). Must be < 0 to form ice
	165	qlead(ji,jj) = MIN( 0._wp , zqld - zqfr )
	166
	167	! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting
	168	IF( at_i(ji,jj) > epsi10 .AND. zqld > 0._wp ) THEN
	169	fhld (ji,jj) = zqld * r1_rdtice / at_i(ji,jj) ! divided by at_i since this is (re)multiplied by a_i in limthd_dh.F90
	170	qlead(ji,jj) = 0._wp
	171	ENDIF
[2528]	172	!
[4688]	173	!-- Energy from the turbulent oceanic heat flux --- !
	174	!clem zfric_u = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin )
	175	zfric_u = MAX( SQRT( ust2s(ji,jj) ), zfric_umin )
	176	fhtur(ji,jj) = MAX( 0._wp, zinda * rau0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ) ! W.m-2
	177	! upper bound for fhtur: we do not want SST to drop below Tfreeze.
	178	! So we say that the heat retrieved from the ocean (fhtur+fhld) must be < to the heat necessary to reach Tfreeze (zqfr)
	179	! This is not a clean budget, so that should be corrected at some point
	180	fhtur(ji,jj) = zinda * MIN( fhtur(ji,jj), - fhld(ji,jj) - zqfr * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )
	181
	182	! -----------------------------------------
	183	! Net heat flux on top of ice-ocean [W.m-2]
	184	! -----------------------------------------
	185	! First step here : heat flux at the ocean surface + precip
	186	! Second step below : heat flux at the ice surface (after limthd_dif)
	187	hfx_in(ji,jj) = hfx_in(ji,jj) &
	188	! heat flux above the ocean
	189	& + pfrld(ji,jj) * ( qns(ji,jj) + qsr(ji,jj) ) &
	190	! latent heat of precip (note that precip is included in qns but not in qns_ice)
	191	& + ( 1._wp - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus ) &
	192	& + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt )
	193
	194	! -----------------------------------------------------------------------------
	195	! Net heat flux that is retroceded to the ocean or taken from the ocean [W.m-2]
	196	! -----------------------------------------------------------------------------
	197	! First step here : non solar + precip - qlead - qturb
	198	! Second step in limthd_dh : heat remaining if total melt (zq_rema)
	199	! Third step in limsbc : heat from ice-ocean mass exchange (zf_mass) + solar
[4765]	200	hfx_out(ji,jj) = hfx_out(ji,jj) &
[4688]	201	! Non solar heat flux received by the ocean
[4765]	202	& + pfrld(ji,jj) * qns(ji,jj) &
[4688]	203	! latent heat of precip (note that precip is included in qns but not in qns_ice)
[4765]	204	& + ( pfrld(ji,jj)*betas - pfrld(ji,jj) ) sprecip(ji,jj) &
	205	& * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus ) &
	206	& + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) &
[4688]	207	! heat flux taken from the ocean where there is open water ice formation
[4765]	208	& - qlead(ji,jj) * r1_rdtice &
[4688]	209	! heat flux taken from the ocean during bottom growth/melt (fhld should be 0 while bott growth)
[4765]	210	& - at_i(ji,jj) * fhtur(ji,jj) &
[4688]	211	& - at_i(ji,jj) * fhld(ji,jj)
	212
[825]	213	END DO
	214	END DO
	215
[921]	216	!------------------------------------------------------------------------------!
	217	! 3) Select icy points and fulfill arrays for the vectorial grid.
	218	!------------------------------------------------------------------------------!
[825]	219
	220	DO jl = 1, jpl !loop over ice categories
	221
[921]	222	IF( kt == nit000 .AND. lwp ) THEN
	223	WRITE(numout,*) ' lim_thd : transfer to 1D vectors. Category no : ', jl
	224	WRITE(numout,*) ' ~~~~~~~~'
	225	ENDIF
[825]	226
[4333]	227	zareamin = epsi10
[825]	228	nbpb = 0
	229	DO jj = 1, jpj
	230	DO ji = 1, jpi
	231	IF ( a_i(ji,jj,jl) .gt. zareamin ) THEN
	232	nbpb = nbpb + 1
	233	npb(nbpb) = (jj - 1) * jpi + ji
	234	ENDIF
	235	END DO
	236	END DO
	237
[4333]	238	! debug point to follow
	239	jiindex_1d = 0
	240	IF( ln_nicep ) THEN
	241	DO ji = mi0(jiindx), mi1(jiindx)
	242	DO jj = mj0(jjindx), mj1(jjindx)
	243	jiindex_1d = (jj - 1) * jpi + ji
[4688]	244	WRITE(numout,*) ' lim_thd : Category no : ', jl
[4333]	245	END DO
	246	END DO
	247	ENDIF
	248
[921]	249	!------------------------------------------------------------------------------!
	250	! 4) Thermodynamic computation
	251	!------------------------------------------------------------------------------!
[825]	252
[2715]	253	IF( lk_mpp ) CALL mpp_ini_ice( nbpb , numout )
[869]	254
[1572]	255	IF( nbpb > 0 ) THEN ! If there is no ice, do nothing.
[825]	256
[921]	257	!-------------------------
	258	! 4.1 Move to 1D arrays
	259	!-------------------------
[825]	260
[1572]	261	CALL tab_2d_1d( nbpb, at_i_b (1:nbpb), at_i , jpi, jpj, npb(1:nbpb) )
	262	CALL tab_2d_1d( nbpb, a_i_b (1:nbpb), a_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	263	CALL tab_2d_1d( nbpb, ht_i_b (1:nbpb), ht_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	264	CALL tab_2d_1d( nbpb, ht_s_b (1:nbpb), ht_s(:,:,jl) , jpi, jpj, npb(1:nbpb) )
[825]	265
[1572]	266	CALL tab_2d_1d( nbpb, t_su_b (1:nbpb), t_su(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	267	CALL tab_2d_1d( nbpb, sm_i_b (1:nbpb), sm_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )
[825]	268	DO jk = 1, nlay_s
[1572]	269	CALL tab_2d_1d( nbpb, t_s_b(1:nbpb,jk), t_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
	270	CALL tab_2d_1d( nbpb, q_s_b(1:nbpb,jk), e_s(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
[825]	271	END DO
	272	DO jk = 1, nlay_i
[1572]	273	CALL tab_2d_1d( nbpb, t_i_b(1:nbpb,jk), t_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
	274	CALL tab_2d_1d( nbpb, q_i_b(1:nbpb,jk), e_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
	275	CALL tab_2d_1d( nbpb, s_i_b(1:nbpb,jk), s_i(:,:,jk,jl) , jpi, jpj, npb(1:nbpb) )
[825]	276	END DO
	277
[1572]	278	CALL tab_2d_1d( nbpb, tatm_ice_1d(1:nbpb), tatm_ice(:,:) , jpi, jpj, npb(1:nbpb) )
	279	CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	280	CALL tab_2d_1d( nbpb, fr1_i0_1d (1:nbpb), fr1_i0 , jpi, jpj, npb(1:nbpb) )
	281	CALL tab_2d_1d( nbpb, fr2_i0_1d (1:nbpb), fr2_i0 , jpi, jpj, npb(1:nbpb) )
[4688]	282	CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	283	CALL tab_2d_1d( nbpb, ftr_ice_1d (1:nbpb), ftr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	284	IF( .NOT. lk_cpl ) THEN
	285	CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb), qla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
	286	CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb), dqla_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
	287	ENDIF
[3625]	288	CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
	289	CALL tab_2d_1d( nbpb, t_bo_b (1:nbpb), t_bo , jpi, jpj, npb(1:nbpb) )
	290	CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip , jpi, jpj, npb(1:nbpb) )
[4688]	291	CALL tab_2d_1d( nbpb, fhtur_1d (1:nbpb), fhtur , jpi, jpj, npb(1:nbpb) )
	292	CALL tab_2d_1d( nbpb, qlead_1d (1:nbpb), qlead , jpi, jpj, npb(1:nbpb) )
	293	CALL tab_2d_1d( nbpb, fhld_1d (1:nbpb), fhld , jpi, jpj, npb(1:nbpb) )
[825]	294
[4688]	295	CALL tab_2d_1d( nbpb, wfx_snw_1d (1:nbpb), wfx_snw , jpi, jpj, npb(1:nbpb) )
	296	CALL tab_2d_1d( nbpb, wfx_sub_1d (1:nbpb), wfx_sub , jpi, jpj, npb(1:nbpb) )
	297
	298	CALL tab_2d_1d( nbpb, wfx_bog_1d (1:nbpb), wfx_bog , jpi, jpj, npb(1:nbpb) )
	299	CALL tab_2d_1d( nbpb, wfx_bom_1d (1:nbpb), wfx_bom , jpi, jpj, npb(1:nbpb) )
	300	CALL tab_2d_1d( nbpb, wfx_sum_1d (1:nbpb), wfx_sum , jpi, jpj, npb(1:nbpb) )
	301	CALL tab_2d_1d( nbpb, wfx_sni_1d (1:nbpb), wfx_sni , jpi, jpj, npb(1:nbpb) )
	302	CALL tab_2d_1d( nbpb, wfx_res_1d (1:nbpb), wfx_res , jpi, jpj, npb(1:nbpb) )
	303	CALL tab_2d_1d( nbpb, wfx_spr_1d (1:nbpb), wfx_spr , jpi, jpj, npb(1:nbpb) )
	304
	305	CALL tab_2d_1d( nbpb, sfx_bog_1d (1:nbpb), sfx_bog , jpi, jpj, npb(1:nbpb) )
	306	CALL tab_2d_1d( nbpb, sfx_bom_1d (1:nbpb), sfx_bom , jpi, jpj, npb(1:nbpb) )
	307	CALL tab_2d_1d( nbpb, sfx_sum_1d (1:nbpb), sfx_sum , jpi, jpj, npb(1:nbpb) )
	308	CALL tab_2d_1d( nbpb, sfx_sni_1d (1:nbpb), sfx_sni , jpi, jpj, npb(1:nbpb) )
[3625]	309	CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri , jpi, jpj, npb(1:nbpb) )
[4688]	310	CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res , jpi, jpj, npb(1:nbpb) )
[825]	311
[4688]	312	CALL tab_2d_1d( nbpb, iatte_1d (1:nbpb), iatte , jpi, jpj, npb(1:nbpb) )
	313	CALL tab_2d_1d( nbpb, oatte_1d (1:nbpb), oatte , jpi, jpj, npb(1:nbpb) )
	314
	315	CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd , jpi, jpj, npb(1:nbpb) )
	316	CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr , jpi, jpj, npb(1:nbpb) )
	317	CALL tab_2d_1d( nbpb, hfx_sum_1d (1:nbpb), hfx_sum , jpi, jpj, npb(1:nbpb) )
	318	CALL tab_2d_1d( nbpb, hfx_bom_1d (1:nbpb), hfx_bom , jpi, jpj, npb(1:nbpb) )
	319	CALL tab_2d_1d( nbpb, hfx_bog_1d (1:nbpb), hfx_bog , jpi, jpj, npb(1:nbpb) )
	320	CALL tab_2d_1d( nbpb, hfx_dif_1d (1:nbpb), hfx_dif , jpi, jpj, npb(1:nbpb) )
	321	CALL tab_2d_1d( nbpb, hfx_opw_1d (1:nbpb), hfx_opw , jpi, jpj, npb(1:nbpb) )
	322	CALL tab_2d_1d( nbpb, hfx_snw_1d (1:nbpb), hfx_snw , jpi, jpj, npb(1:nbpb) )
	323	CALL tab_2d_1d( nbpb, hfx_sub_1d (1:nbpb), hfx_sub , jpi, jpj, npb(1:nbpb) )
	324	CALL tab_2d_1d( nbpb, hfx_err_1d (1:nbpb), hfx_err , jpi, jpj, npb(1:nbpb) )
	325	CALL tab_2d_1d( nbpb, hfx_res_1d (1:nbpb), hfx_res , jpi, jpj, npb(1:nbpb) )
	326	CALL tab_2d_1d( nbpb, hfx_err_rem_1d (1:nbpb), hfx_err_rem , jpi, jpj, npb(1:nbpb) )
	327
[921]	328	!--------------------------------
	329	! 4.3) Thermodynamic processes
	330	!--------------------------------
	331
[4688]	332	!---------------------------------!
	333	! Ice/Snow Temperature profile !
	334	!---------------------------------!
	335	CALL lim_thd_dif( 1, nbpb )
[921]	336
[4688]	337	!---------------------------------!
	338	! Ice/Snow thicnkess !
	339	!---------------------------------!
	340	CALL lim_thd_dh( 1, nbpb )
[825]	341
[4688]	342	! --- Ice enthalpy remapping --- !
	343	CALL lim_thd_ent( 1, nbpb, q_i_b(1:nbpb,:) )
	344
	345	!---------------------------------!
	346	! --- Ice salinity --- !
	347	!---------------------------------!
	348	CALL lim_thd_sal( 1, nbpb )
[825]	349
[4688]	350	!---------------------------------!
	351	! --- temperature update --- !
	352	!---------------------------------!
	353	CALL lim_thd_temp( 1, nbpb )
[825]	354
[921]	355	!--------------------------------
	356	! 4.4) Move 1D to 2D vectors
	357	!--------------------------------
[825]	358
[3625]	359	CALL tab_1d_2d( nbpb, at_i , npb, at_i_b (1:nbpb) , jpi, jpj )
	360	CALL tab_1d_2d( nbpb, ht_i(:,:,jl) , npb, ht_i_b (1:nbpb) , jpi, jpj )
	361	CALL tab_1d_2d( nbpb, ht_s(:,:,jl) , npb, ht_s_b (1:nbpb) , jpi, jpj )
	362	CALL tab_1d_2d( nbpb, a_i (:,:,jl) , npb, a_i_b (1:nbpb) , jpi, jpj )
	363	CALL tab_1d_2d( nbpb, t_su(:,:,jl) , npb, t_su_b (1:nbpb) , jpi, jpj )
	364	CALL tab_1d_2d( nbpb, sm_i(:,:,jl) , npb, sm_i_b (1:nbpb) , jpi, jpj )
[825]	365	DO jk = 1, nlay_s
[3625]	366	CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_b (1:nbpb,jk), jpi, jpj)
	367	CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_b (1:nbpb,jk), jpi, jpj)
[825]	368	END DO
	369	DO jk = 1, nlay_i
[3625]	370	CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_b (1:nbpb,jk), jpi, jpj)
	371	CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_b (1:nbpb,jk), jpi, jpj)
	372	CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_b (1:nbpb,jk), jpi, jpj)
[825]	373	END DO
[4688]	374	CALL tab_1d_2d( nbpb, qlead , npb, qlead_1d (1:nbpb) , jpi, jpj )
	375
	376	CALL tab_1d_2d( nbpb, wfx_snw , npb, wfx_snw_1d(1:nbpb) , jpi, jpj )
	377	CALL tab_1d_2d( nbpb, wfx_sub , npb, wfx_sub_1d(1:nbpb) , jpi, jpj )
	378
	379	CALL tab_1d_2d( nbpb, wfx_bog , npb, wfx_bog_1d(1:nbpb) , jpi, jpj )
	380	CALL tab_1d_2d( nbpb, wfx_bom , npb, wfx_bom_1d(1:nbpb) , jpi, jpj )
	381	CALL tab_1d_2d( nbpb, wfx_sum , npb, wfx_sum_1d(1:nbpb) , jpi, jpj )
	382	CALL tab_1d_2d( nbpb, wfx_sni , npb, wfx_sni_1d(1:nbpb) , jpi, jpj )
	383	CALL tab_1d_2d( nbpb, wfx_res , npb, wfx_res_1d(1:nbpb) , jpi, jpj )
	384	CALL tab_1d_2d( nbpb, wfx_spr , npb, wfx_spr_1d(1:nbpb) , jpi, jpj )
	385
	386	CALL tab_1d_2d( nbpb, sfx_bog , npb, sfx_bog_1d(1:nbpb) , jpi, jpj )
	387	CALL tab_1d_2d( nbpb, sfx_bom , npb, sfx_bom_1d(1:nbpb) , jpi, jpj )
	388	CALL tab_1d_2d( nbpb, sfx_sum , npb, sfx_sum_1d(1:nbpb) , jpi, jpj )
	389	CALL tab_1d_2d( nbpb, sfx_sni , npb, sfx_sni_1d(1:nbpb) , jpi, jpj )
	390	CALL tab_1d_2d( nbpb, sfx_res , npb, sfx_res_1d(1:nbpb) , jpi, jpj )
[2528]	391	!
[1572]	392	IF( num_sal == 2 ) THEN
[3625]	393	CALL tab_1d_2d( nbpb, sfx_bri , npb, sfx_bri_1d(1:nbpb) , jpi, jpj )
[825]	394	ENDIF
[4688]	395
	396	CALL tab_1d_2d( nbpb, hfx_thd , npb, hfx_thd_1d(1:nbpb) , jpi, jpj )
	397	CALL tab_1d_2d( nbpb, hfx_spr , npb, hfx_spr_1d(1:nbpb) , jpi, jpj )
	398	CALL tab_1d_2d( nbpb, hfx_sum , npb, hfx_sum_1d(1:nbpb) , jpi, jpj )
	399	CALL tab_1d_2d( nbpb, hfx_bom , npb, hfx_bom_1d(1:nbpb) , jpi, jpj )
	400	CALL tab_1d_2d( nbpb, hfx_bog , npb, hfx_bog_1d(1:nbpb) , jpi, jpj )
	401	CALL tab_1d_2d( nbpb, hfx_dif , npb, hfx_dif_1d(1:nbpb) , jpi, jpj )
	402	CALL tab_1d_2d( nbpb, hfx_opw , npb, hfx_opw_1d(1:nbpb) , jpi, jpj )
	403	CALL tab_1d_2d( nbpb, hfx_snw , npb, hfx_snw_1d(1:nbpb) , jpi, jpj )
	404	CALL tab_1d_2d( nbpb, hfx_sub , npb, hfx_sub_1d(1:nbpb) , jpi, jpj )
	405	CALL tab_1d_2d( nbpb, hfx_err , npb, hfx_err_1d(1:nbpb) , jpi, jpj )
	406	CALL tab_1d_2d( nbpb, hfx_res , npb, hfx_res_1d(1:nbpb) , jpi, jpj )
	407	CALL tab_1d_2d( nbpb, hfx_err_rem , npb, hfx_err_rem_1d(1:nbpb) , jpi, jpj )
[2528]	408	!
[3625]	409	!+++++ temporary stuff for a dummy version
[4688]	410	CALL tab_1d_2d( nbpb, dh_i_surf2D, npb, dh_i_surf(1:nbpb) , jpi, jpj )
	411	CALL tab_1d_2d( nbpb, dh_i_bott2D, npb, dh_i_bott(1:nbpb) , jpi, jpj )
	412	CALL tab_1d_2d( nbpb, s_i_newice , npb, s_i_new (1:nbpb) , jpi, jpj )
	413	CALL tab_1d_2d( nbpb, izero(:,:,jl) , npb, i0 (1:nbpb) , jpi, jpj )
[825]	414	!+++++
[4688]	415	CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj)
	416	CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj )
[2528]	417	!
[1572]	418	IF( lk_mpp ) CALL mpp_comm_free( ncomm_ice ) !RB necessary ??
	419	ENDIF
	420	!
	421	END DO
[825]	422
[921]	423	!------------------------------------------------------------------------------!
	424	! 5) Global variables, diagnostics
	425	!------------------------------------------------------------------------------!
[825]	426
	427	!------------------------
	428	! 5.1) Ice heat content
	429	!------------------------
[4688]	430	! Enthalpies are global variables we have to readjust the units (heat content in Joules)
[825]	431	DO jl = 1, jpl
[921]	432	DO jk = 1, nlay_i
[4688]	433	e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_i(:,:,jl) / ( unit_fac * REAL( nlay_i ) )
[1572]	434	END DO
	435	END DO
[825]	436
	437	!------------------------
	438	! 5.2) Snow heat content
	439	!------------------------
[4688]	440	! Enthalpies are global variables we have to readjust the units (heat content in Joules)
[825]	441	DO jl = 1, jpl
	442	DO jk = 1, nlay_s
[4688]	443	e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_s(:,:,jl) / ( unit_fac * REAL( nlay_s ) )
[1572]	444	END DO
	445	END DO
[825]	446
	447	!----------------------------------
	448	! 5.3) Change thickness to volume
	449	!----------------------------------
	450	CALL lim_var_eqv2glo
	451
	452	!--------------------------------------------
	453	! 5.4) Diagnostic thermodynamic growth rates
	454	!--------------------------------------------
[2528]	455	IF(ln_ctl) THEN ! Control print
[867]	456	CALL prt_ctl_info(' ')
	457	CALL prt_ctl_info(' - Cell values : ')
	458	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
[863]	459	CALL prt_ctl(tab2d_1=area , clinfo1=' lim_thd : cell area :')
	460	CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_thd : at_i :')
	461	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_thd : vt_i :')
	462	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_thd : vt_s :')
	463	DO jl = 1, jpl
[867]	464	CALL prt_ctl_info(' ')
[863]	465	CALL prt_ctl_info(' - Category : ', ivar1=jl)
	466	CALL prt_ctl_info(' ~~~~~~~~~~')
	467	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_thd : a_i : ')
	468	CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_thd : ht_i : ')
	469	CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_thd : ht_s : ')
	470	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_thd : v_i : ')
	471	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_thd : v_s : ')
	472	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_thd : e_s : ')
	473	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_thd : t_su : ')
	474	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_thd : t_snow : ')
	475	CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_thd : sm_i : ')
	476	CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_thd : smv_i : ')
	477	DO jk = 1, nlay_i
[867]	478	CALL prt_ctl_info(' ')
[863]	479	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
	480	CALL prt_ctl_info(' ~~~~~~~')
	481	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_thd : t_i : ')
	482	CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_thd : e_i : ')
	483	END DO
	484	END DO
	485	ENDIF
[2528]	486	!
[4688]	487	! conservation test
	488	IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[4161]	489	!
	490	IF( nn_timing == 1 ) CALL timing_stop('limthd')
[4688]	491	END SUBROUTINE lim_thd
[825]	492
[4688]	493	SUBROUTINE lim_thd_temp( kideb, kiut )
[825]	494	!!-----------------------------------------------------------------------
[4688]	495	!! * ROUTINE lim_thd_temp *
[825]	496	!!
[4688]	497	!! ** Purpose : Computes sea ice temperature (Kelvin) from enthalpy
[825]	498	!!
	499	!! ** Method : Formula (Bitz and Lipscomb, 1999)
	500	!!-------------------------------------------------------------------
[1572]	501	INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop
	502	!!
[2715]	503	INTEGER :: ji, jk ! dummy loop indices
[4688]	504	REAL(wp) :: ztmelts, zswitch, zaaa, zbbb, zccc, zdiscrim ! local scalar
[825]	505	!!-------------------------------------------------------------------
[4688]	506	! Recover ice temperature
	507	DO jk = 1, nlay_i
[825]	508	DO ji = kideb, kiut
[4688]	509	ztmelts = -tmut * s_i_b(ji,jk) + rtt
	510	! Conversion q(S,T) -> T (second order equation)
	511	zaaa = cpic
	512	zbbb = ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_b(ji,jk) / rhoic - lfus
	513	zccc = lfus * ( ztmelts - rtt )
	514	zdiscrim = SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) )
	515	t_i_b(ji,jk) = rtt - ( zbbb + zdiscrim ) / ( 2._wp * zaaa )
	516
	517	! mask temperature
	518	zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_b(ji) ) )
	519	t_i_b(ji,jk) = zswitch * t_i_b(ji,jk) + ( 1._wp - zswitch ) * rtt
	520	END DO
	521	END DO
[825]	522
[4688]	523	END SUBROUTINE lim_thd_temp
[825]	524
	525	SUBROUTINE lim_thd_init
	526	!!-----------------------------------------------------------------------
	527	!! * ROUTINE lim_thd_init *
	528	!!
	529	!! ** Purpose : Physical constants and parameters linked to the ice
[1572]	530	!! thermodynamics
[825]	531	!!
	532	!! ** Method : Read the namicethd namelist and check the ice-thermo
[1572]	533	!! parameter values called at the first timestep (nit000)
[825]	534	!!
	535	!! ** input : Namelist namicether
[2528]	536	!!-------------------------------------------------------------------
[4147]	537	INTEGER :: ios ! Local integer output status for namelist read
[1572]	538	NAMELIST/namicethd/ hmelt , hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb, &
[4688]	539	& hiclim, hnzst, parsub, betas, &
[825]	540	& kappa_i, nconv_i_thd, maxer_i_thd, thcon_i_swi
	541	!!-------------------------------------------------------------------
[2528]	542	!
[1572]	543	IF(lwp) THEN
	544	WRITE(numout,*)
	545	WRITE(numout,*) 'lim_thd : Ice Thermodynamics'
	546	WRITE(numout,*) '~~~~~~~'
	547	ENDIF
[2528]	548	!
[4147]	549	REWIND( numnam_ice_ref ) ! Namelist namicethd in reference namelist : Ice thermodynamics
	550	READ ( numnam_ice_ref, namicethd, IOSTAT = ios, ERR = 901)
	551	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in reference namelist', lwp )
	552
	553	REWIND( numnam_ice_cfg ) ! Namelist namicethd in configuration namelist : Ice thermodynamics
	554	READ ( numnam_ice_cfg, namicethd, IOSTAT = ios, ERR = 902 )
	555	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in configuration namelist', lwp )
[4624]	556	IF(lwm) WRITE ( numoni, namicethd )
[2528]	557	!
[1572]	558	IF(lwp) THEN ! control print
[825]	559	WRITE(numout,*)
[1572]	560	WRITE(numout,*)' Namelist of ice parameters for ice thermodynamic computation '
	561	WRITE(numout,*)' maximum melting at the bottom hmelt = ', hmelt
[4688]	562	WRITE(numout,*)' ice thick. for lateral accretion hiccrit = ', hiccrit
[1572]	563	WRITE(numout,*)' Frazil ice thickness as a function of wind or not fraz_swi = ', fraz_swi
	564	WRITE(numout,*)' Maximum proportion of frazil ice collecting at bottom maxfrazb = ', maxfrazb
	565	WRITE(numout,*)' Thresold relative drift speed for collection of frazil vfrazb = ', vfrazb
	566	WRITE(numout,*)' Squeezing coefficient for collection of frazil Cfrazb = ', Cfrazb
	567	WRITE(numout,*)' minimum ice thickness hiclim = ', hiclim
	568	WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice '
	569	WRITE(numout,*)' thickness of the surf. layer in temp. computation hnzst = ', hnzst
	570	WRITE(numout,*)' switch for snow sublimation (=1) or not (=0) parsub = ', parsub
	571	WRITE(numout,*)' coefficient for ice-lead partition of snowfall betas = ', betas
	572	WRITE(numout,*)' extinction radiation parameter in sea ice (1.0) kappa_i = ', kappa_i
	573	WRITE(numout,*)' maximal n. of iter. for heat diffusion computation nconv_i_thd = ', nconv_i_thd
	574	WRITE(numout,*)' maximal err. on T for heat diffusion computation maxer_i_thd = ', maxer_i_thd
	575	WRITE(numout,*)' switch for comp. of thermal conductivity in the ice thcon_i_swi = ', thcon_i_swi
[4688]	576	WRITE(numout,*)' check heat conservation in the ice/snow con_i = ', con_i
[825]	577	ENDIF
[1572]	578	!
[825]	579	END SUBROUTINE lim_thd_init
	580
	581	#else
[1572]	582	!!----------------------------------------------------------------------
[2528]	583	!! Default option Dummy module NO LIM3 sea-ice model
[1572]	584	!!----------------------------------------------------------------------
[825]	585	#endif
	586
	587	!!======================================================================
	588	END MODULE limthd

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

Download in other formats: