Context Navigation

limthd_2.F90 @ 1858

Last change on this file since 1858 was 1858, checked in by gm, 14 years ago
ticket:#665 : step 1 - heat content of freezing-melting ice
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 33.2 KB

Rev	Line
[821]	1	MODULE limthd_2
[3]	2	!!======================================================================
[821]	3	!! * MODULE limthd_2 *
[1857]	4	!! LIM thermo ice model : ice thermodynamic
[3]	5	!!======================================================================
[1858]	6	!! History : 1.0 ! 2000-01 (LIM) Original code
	7	!! 2.0 ! 2002-07 (C. Ethe, G. Madec) F90
	8	!! 2.0 ! 2003-08 (C. Ethe) add lim_thd_init
	9	!! - ! 2008-08 (A. Caubel, G. Madec, E. Maisonnave, S. Masson ) generic coupled interface
	10	!! 2.1 ! 2010-05 (Y. Aksenov G. Madec) add heat associated with mass exchanges
[888]	11	!!---------------------------------------------------------------------
[821]	12	#if defined key_lim2
[3]	13	!!----------------------------------------------------------------------
[821]	14	!! 'key_lim2' : LIM 2.0 sea-ice model
[3]	15	!!----------------------------------------------------------------------
[1857]	16	!! lim_thd_2 : thermodynamic of sea ice
[821]	17	!! lim_thd_init_2 : initialisation of sea-ice thermodynamic
[3]	18	!!----------------------------------------------------------------------
[717]	19	USE phycst ! physical constants
[3]	20	USE dom_oce ! ocean space and time domain variables
[1758]	21	USE domvvl
[717]	22	USE lbclnk
	23	USE in_out_manager ! I/O manager
[1756]	24	USE lib_mpp
[1482]	25	USE iom ! IOM library
[1857]	26	USE ice_2 ! LIM sea-ice variables
[888]	27	USE sbc_oce !
	28	USE sbc_ice !
[1857]	29	USE thd_ice_2 ! LIM thermodynamic sea-ice variables
	30	USE dom_ice_2 ! LIM sea-ice domain
[821]	31	USE limthd_zdf_2
	32	USE limthd_lac_2
	33	USE limtab_2
[258]	34	USE prtctl ! Print control
[1218]	35	USE cpl_oasis3, ONLY : lk_cpl
[1756]	36	USE diaar5, ONLY : lk_diaar5
[3]	37
	38	IMPLICIT NONE
	39	PRIVATE
	40
[888]	41	PUBLIC lim_thd_2 ! called by lim_step
[3]	42
[1218]	43	REAL(wp) :: epsi20 = 1.e-20 ! constant values
	44	REAL(wp) :: epsi16 = 1.e-16 !
	45	REAL(wp) :: epsi04 = 1.e-04 !
	46	REAL(wp) :: rzero = 0.e0 !
	47	REAL(wp) :: rone = 1.e0 !
[70]	48
[3]	49	!! * Substitutions
[70]	50	# include "domzgr_substitute.h90"
[3]	51	# include "vectopt_loop_substitute.h90"
[1857]	52	!!-------- -------------------------------------------------------------
	53	!! NEMO/LIM 3.2, UCL-LOCEAN-IPSL (2009)
[1156]	54	!! $Id$
[888]	55	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[3]	56	!!----------------------------------------------------------------------
	57
	58	CONTAINS
	59
[821]	60	SUBROUTINE lim_thd_2( kt )
[3]	61	!!-------------------------------------------------------------------
[821]	62	!! * ROUTINE lim_thd_2 *
[3]	63	!!
	64	!! ** Purpose : This routine manages the ice thermodynamic.
	65	!!
	66	!! ** Action : - Initialisation of some variables
	67	!! - Some preliminary computation (oceanic heat flux
	68	!! at the ice base, snow acc.,heat budget of the leads)
	69	!! - selection of the icy points and put them in an array
	70	!! - call lim_vert_ther for vert ice thermodynamic
	71	!! - back to the geographic grid
	72	!! - selection of points for lateral accretion
	73	!! - call lim_lat_acc for the ice accretion
	74	!! - back to the geographic grid
	75	!!
[888]	76	!! References : Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90
[3]	77	!!---------------------------------------------------------------------
[508]	78	INTEGER, INTENT(in) :: kt ! number of iteration
[888]	79	!!
[1218]	80	INTEGER :: ji, jj ! dummy loop indices
	81	INTEGER :: nbpb ! nb of icy pts for thermo. cal.
	82	INTEGER :: nbpac ! nb of pts for lateral accretion
[258]	83	CHARACTER (len=22) :: charout
[1218]	84	REAL(wp) :: zfric_umin = 5e-03 ! lower bound for the friction velocity
	85	REAL(wp) :: zfric_umax = 2e-02 ! upper bound for the friction velocity
	86	REAL(wp) :: zinda ! switch for test. the val. of concen.
	87	REAL(wp) :: zindb, zindg ! switches for test. the val of arg
[1559]	88	REAL(wp) :: zfricp ! temporary scalar
[1218]	89	REAL(wp) :: za , zh, zthsnice !
	90	REAL(wp) :: zfric_u ! friction velocity
	91	REAL(wp) :: zfntlat, zpareff ! test. the val. of lead heat budget
[1559]	92	REAL(wp), DIMENSION(jpi,jpj) :: ztmp ! 2D workspace
[1218]	93	REAL(wp), DIMENSION(jpi,jpj) :: zqlbsbq ! link with lead energy budget qldif
[1756]	94	REAL(wp) :: zuice_m, zvice_m ! Sea-ice velocities at U & V-points
	95	REAL(wp) :: zhice_u, zhice_v ! Sea-ice volume at U & V-points
	96	REAL(wp) :: ztr_fram ! Sea-ice transport through Fram strait
	97	REAL(wp) :: zrhoij, zrhoijm1 ! temporary scalars
	98	REAL(wp) :: zztmp ! temporary scalars within a loop
	99	REAL(wp), DIMENSION(jpi,jpj) :: zlicegr ! link with lateral ice growth
[1559]	100	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmsk ! 3D workspace
[1756]	101	!!$ REAL(wp), DIMENSION(jpi,jpj) :: firic !: IR flux over the ice (outputs only)
	102	!!$ REAL(wp), DIMENSION(jpi,jpj) :: fcsic !: Sensible heat flux over the ice (outputs only)
	103	!!$ REAL(wp), DIMENSION(jpi,jpj) :: fleic !: Latent heat flux over the ice (outputs only)
	104	!!$ REAL(wp), DIMENSION(jpi,jpj) :: qlatic !: latent flux (outputs only)
	105	REAL(wp), DIMENSION(jpi,jpj) :: zdvosif !: Variation of volume at surface (outputs only)
	106	REAL(wp), DIMENSION(jpi,jpj) :: zdvobif !: Variation of ice volume at the bottom ice (outputs only)
	107	REAL(wp), DIMENSION(jpi,jpj) :: zdvolif !: Total variation of ice volume (outputs only)
	108	REAL(wp), DIMENSION(jpi,jpj) :: zdvonif !: Surface accretion Snow to Ice transformation (outputs only)
	109	REAL(wp), DIMENSION(jpi,jpj) :: zdvomif !: Bottom variation of ice volume due to melting (outputs only)
	110	REAL(wp), DIMENSION(jpi,jpj) :: zu_imasstr !: Sea-ice transport along i-axis at U-point (outputs only)
	111	REAL(wp), DIMENSION(jpi,jpj) :: zv_imasstr !: Sea-ice transport along j-axis at V-point (outputs only)
[3]	112	!!-------------------------------------------------------------------
	113
[888]	114	IF( kt == nit000 ) CALL lim_thd_init_2 ! Initialization (first time-step only)
[3]	115
	116	!-------------------------------------------!
	117	! Initilization of diagnostic variables !
	118	!-------------------------------------------!
	119
[1218]	120	!!gm needed? yes at least for some of these arrays
[1756]	121	zdvosif(:,:) = 0.e0 ! variation of ice volume at surface
	122	zdvobif(:,:) = 0.e0 ! variation of ice volume at bottom
	123	zdvolif(:,:) = 0.e0 ! total variation of ice volume
	124	zdvonif(:,:) = 0.e0 ! transformation of snow to sea-ice volume
	125	! zdvonif(:,:) = 0.e0 ! lateral variation of ice volume
	126	zlicegr(:,:) = 0.e0 ! lateral variation of ice volume
	127	zdvomif(:,:) = 0.e0 ! variation of ice volume at bottom due to melting only
	128	ztr_fram = 0.e0 ! sea-ice transport through Fram strait
[70]	129	fstric (:,:) = 0.e0 ! part of solar radiation absorbing inside the ice
	130	fscmbq (:,:) = 0.e0 ! linked with fstric
	131	ffltbif(:,:) = 0.e0 ! linked with fstric
	132	qfvbq (:,:) = 0.e0 ! linked with fstric
[1858]	133	rdm_snw(:,:) = 0.e0 ! variation of snow mass over 1 time step
	134	rdq_snw(:,:) = 0.e0 ! heat content associated with rdm_snw
	135	rdm_ice(:,:) = 0.e0 ! variation of ice mass over 1 time step
	136	rdq_ice(:,:) = 0.e0 ! heat content associated with rdm_ice
[258]	137	zmsk (:,:,:) = 0.e0
[3]	138
[1218]	139	! set to zero snow thickness smaller than epsi04
[3]	140	DO jj = 1, jpj
	141	DO ji = 1, jpi
[1218]	142	hsnif(ji,jj) = hsnif(ji,jj) * MAX( rzero, SIGN( rone , hsnif(ji,jj) - epsi04 ) )
[3]	143	END DO
	144	END DO
[1218]	145	!!gm better coded (do not use SIGN...)
	146	! WHERE( hsnif(:,:) < epsi04 ) hsnif(:,:) = 0.e0
	147	!!gm
[258]	148
[1218]	149	IF(ln_ctl) CALL prt_ctl( tab2d_1=hsnif, clinfo1=' lim_thd: hsnif : ' )
[3]	150
	151	!-----------------------------------!
	152	! Treatment of particular cases !
	153	!-----------------------------------!
	154
	155	DO jj = 1, jpj
	156	DO ji = 1, jpi
	157	! snow is transformed into ice if the original ice cover disappears.
[1218]	158	zindg = tms(ji,jj) * MAX( rzero , SIGN( rone , -hicif(ji,jj) ) )
[3]	159	hicif(ji,jj) = hicif(ji,jj) + zindg * rhosn * hsnif(ji,jj) / rau0
[1218]	160	hsnif(ji,jj) = ( rone - zindg ) * hsnif(ji,jj) + zindg * hicif(ji,jj) * ( rau0 - rhoic ) / rhosn
[3]	161	dmgwi(ji,jj) = zindg * (1.0 - frld(ji,jj)) * rhoic * hicif(ji,jj) ! snow/ice mass
	162
	163	! the lead fraction, frld, must be little than or equal to amax (ice ridging).
	164	zthsnice = hsnif(ji,jj) + hicif(ji,jj)
[1218]	165	zindb = tms(ji,jj) * ( 1.0 - MAX( rzero , SIGN( rone , - zthsnice ) ) )
	166	za = zindb * MIN( rone, ( 1.0 - frld(ji,jj) ) * uscomi )
[3]	167	hsnif (ji,jj) = hsnif(ji,jj) * za
	168	hicif (ji,jj) = hicif(ji,jj) * za
	169	qstoif(ji,jj) = qstoif(ji,jj) * za
[1218]	170	frld (ji,jj) = 1.0 - zindb * ( 1.0 - frld(ji,jj) ) / MAX( za, epsi20 )
[3]	171
	172	! the in situ ice thickness, hicif, must be equal to or greater than hiclim.
[1218]	173	zh = MAX( rone , zindb * hiclim / MAX( hicif(ji,jj), epsi20 ) )
[3]	174	hsnif (ji,jj) = hsnif(ji,jj) * zh
	175	hicif (ji,jj) = hicif(ji,jj) * zh
	176	qstoif(ji,jj) = qstoif(ji,jj) * zh
	177	frld (ji,jj) = ( frld(ji,jj) + ( zh - 1.0 ) ) / zh
	178	END DO
	179	END DO
[258]	180
	181	IF(ln_ctl) THEN
[1218]	182	CALL prt_ctl( tab2d_1=hicif , clinfo1=' lim_thd: hicif : ' )
	183	CALL prt_ctl( tab2d_1=hsnif , clinfo1=' lim_thd: hsnif : ' )
	184	CALL prt_ctl( tab2d_1=dmgwi , clinfo1=' lim_thd: dmgwi : ' )
	185	CALL prt_ctl( tab2d_1=qstoif, clinfo1=' lim_thd: qstoif : ' )
	186	CALL prt_ctl( tab2d_1=frld , clinfo1=' lim_thd: frld : ' )
[3]	187	ENDIF
	188
	189
	190	!-------------------------------!
	191	! Thermodynamics of sea ice !
	192	!-------------------------------!
	193
	194	! Partial computation of forcing for the thermodynamic sea ice model.
	195	!--------------------------------------------------------------------------
[1858]	196	!CDIR NOVERRCHK
[3]	197	DO jj = 1, jpj
[1858]	198	!CDIR NOVERRCHK
[3]	199	DO ji = 1, jpi
	200	zthsnice = hsnif(ji,jj) + hicif(ji,jj)
[1218]	201	zindb = tms(ji,jj) * ( 1.0 - MAX( rzero , SIGN( rone , - zthsnice ) ) )
[3]	202	pfrld(ji,jj) = frld(ji,jj)
[1559]	203	zfricp = 1.0 - frld(ji,jj)
	204	zinda = 1.0 - MAX( rzero , SIGN( rone , - zfricp ) )
[3]	205
	206	! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
[70]	207	thcm(ji,jj) = 0.e0
[3]	208
	209	! net downward heat flux from the ice to the ocean, expressed as a function of ocean
	210	! temperature and turbulent mixing (McPhee, 1992)
[1858]	211	zfric_u = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin ) ! friction velocity
	212	fdtcn(ji,jj) = zindb * rau0 * rcp * 0.006 * zfric_u * ( sst_m(ji,jj) + rt0 - tfu(ji,jj) )
[3]	213	qdtcn(ji,jj) = zindb * fdtcn(ji,jj) * frld(ji,jj) * rdt_ice
	214
	215	! partial computation of the lead energy budget (qldif)
[1218]	216	#if defined key_coupled
[1858]	217	qldif(ji,jj) = tms(ji,jj) * rdt_ice &
[1559]	218	& * ( ( qsr_tot(ji,jj) - qsr_ice(ji,jj,1) * zfricp ) * ( 1.0 - thcm(ji,jj) ) &
	219	& + ( qns_tot(ji,jj) - qns_ice(ji,jj,1) * zfricp ) &
[1218]	220	& + frld(ji,jj) * ( fdtcn(ji,jj) + ( 1.0 - zindb ) * fsbbq(ji,jj) ) )
	221	#else
[1858]	222	qldif(ji,jj) = tms(ji,jj) * rdt_ice * frld(ji,jj) &
	223	& * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) ) &
	224	& + qns(ji,jj) + fdtcn(ji,jj) &
	225	& + ( 1.0 - zindb ) * fsbbq(ji,jj) )
[1218]	226	#endif
[3]	227	! parlat : percentage of energy used for lateral ablation (0.0)
[1218]	228	zfntlat = 1.0 - MAX( rzero , SIGN( rone , - qldif(ji,jj) ) )
[3]	229	zpareff = 1.0 + ( parlat - 1.0 ) * zinda * zfntlat
	230	zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / MAX( (1.0 - frld(ji,jj)) * rdt_ice , epsi16 )
	231	qldif (ji,jj) = zpareff * qldif(ji,jj)
	232	qdtcn (ji,jj) = zpareff * qdtcn(ji,jj)
	233
	234	! energy needed to bring ocean surface layer until its freezing
[1858]	235	qcmif (ji,jj) = rau0 * rcp * fse3t_m(ji,jj,1) * ( tfu(ji,jj) - sst_m(ji,jj) - rt0 ) * ( 1 - zinda )
[3]	236
	237	! calculate oceanic heat flux.
	238	fbif (ji,jj) = zindb * ( fsbbq(ji,jj) / MAX( (1.0 - frld(ji,jj)) , epsi20 ) + fdtcn(ji,jj) )
	239
[1482]	240	! computation of the thermodynamic ice production (only needed for output)
	241	hicifp(ji,jj) = hicif(ji,jj) * ( 1.0 - frld(ji,jj) )
[3]	242	END DO
	243	END DO
[1858]	244
[3]	245	! Select icy points and fulfill arrays for the vectorial grid.
	246	!----------------------------------------------------------------------
	247	nbpb = 0
	248	DO jj = 1, jpj
	249	DO ji = 1, jpi
	250	IF ( frld(ji,jj) < 1.0 ) THEN
	251	nbpb = nbpb + 1
	252	npb(nbpb) = (jj - 1) * jpi + ji
	253	ENDIF
	254	END DO
	255	END DO
[258]	256
	257	IF(ln_ctl) THEN
	258	CALL prt_ctl(tab2d_1=pfrld, clinfo1=' lim_thd: pfrld : ', tab2d_2=thcm , clinfo2=' thcm : ')
	259	CALL prt_ctl(tab2d_1=fdtcn, clinfo1=' lim_thd: fdtcn : ', tab2d_2=qdtcn , clinfo2=' qdtcn : ')
	260	CALL prt_ctl(tab2d_1=qldif, clinfo1=' lim_thd: qldif : ', tab2d_2=zqlbsbq, clinfo2=' zqlbsbq : ')
	261	CALL prt_ctl(tab2d_1=qcmif, clinfo1=' lim_thd: qcmif : ', tab2d_2=fbif , clinfo2=' fbif : ')
	262	zmsk(:,:,1) = tms(:,:)
[1482]	263	CALL prt_ctl(tab2d_1=qcmif , clinfo1=' lim_thd: qcmif : ', mask1=zmsk)
	264	CALL prt_ctl(tab2d_1=hicifp, clinfo1=' lim_thd: hicifp : ')
[258]	265	WRITE(charout, FMT="('lim_thd: nbpb = ',I4)") nbpb
	266	CALL prt_ctl_info(charout)
[3]	267	ENDIF
	268
	269
	270	! If there is no ice, do nothing. Otherwise, compute Top and Bottom accretion/ablation
	271	!------------------------------------------------------------------------------------
	272
[1218]	273	IF( nbpb > 0 ) THEN
	274	!
[3]	275	! put the variable in a 1-D array for thermodynamics process
[1463]	276	CALL tab_2d_1d_2( nbpb, frld_1d (1:nbpb) , frld , jpi, jpj, npb(1:nbpb) )
	277	CALL tab_2d_1d_2( nbpb, h_ice_1d (1:nbpb) , hicif , jpi, jpj, npb(1:nbpb) )
	278	CALL tab_2d_1d_2( nbpb, h_snow_1d (1:nbpb) , hsnif , jpi, jpj, npb(1:nbpb) )
	279	CALL tab_2d_1d_2( nbpb, sist_1d (1:nbpb) , sist , jpi, jpj, npb(1:nbpb) )
	280	CALL tab_2d_1d_2( nbpb, tbif_1d (1:nbpb , 1 ), tbif(:,:,1) , jpi, jpj, npb(1:nbpb) )
	281	CALL tab_2d_1d_2( nbpb, tbif_1d (1:nbpb , 2 ), tbif(:,:,2) , jpi, jpj, npb(1:nbpb) )
	282	CALL tab_2d_1d_2( nbpb, tbif_1d (1:nbpb , 3 ), tbif(:,:,3) , jpi, jpj, npb(1:nbpb) )
	283	CALL tab_2d_1d_2( nbpb, qsr_ice_1d (1:nbpb) , qsr_ice(:,:,1) , jpi, jpj, npb(1:nbpb) )
	284	CALL tab_2d_1d_2( nbpb, fr1_i0_1d (1:nbpb) , fr1_i0 , jpi, jpj, npb(1:nbpb) )
	285	CALL tab_2d_1d_2( nbpb, fr2_i0_1d (1:nbpb) , fr2_i0 , jpi, jpj, npb(1:nbpb) )
	286	CALL tab_2d_1d_2( nbpb, qns_ice_1d(1:nbpb) , qns_ice(:,:,1), jpi, jpj, npb(1:nbpb) )
	287	CALL tab_2d_1d_2( nbpb, dqns_ice_1d(1:nbpb) , dqns_ice(:,:,1), jpi, jpj, npb(1:nbpb) )
[1218]	288	IF( .NOT. lk_cpl ) THEN
[1463]	289	CALL tab_2d_1d_2( nbpb, qla_ice_1d (1:nbpb) , qla_ice(:,:,1), jpi, jpj, npb(1:nbpb) )
	290	CALL tab_2d_1d_2( nbpb, dqla_ice_1d(1:nbpb) , dqla_ice(:,:,1), jpi, jpj, npb(1:nbpb) )
[1218]	291	ENDIF
[821]	292	CALL tab_2d_1d_2( nbpb, tfu_1d (1:nbpb) , tfu , jpi, jpj, npb(1:nbpb) )
	293	CALL tab_2d_1d_2( nbpb, sprecip_1d (1:nbpb) , sprecip , jpi, jpj, npb(1:nbpb) )
	294	CALL tab_2d_1d_2( nbpb, fbif_1d (1:nbpb) , fbif , jpi, jpj, npb(1:nbpb) )
	295	CALL tab_2d_1d_2( nbpb, thcm_1d (1:nbpb) , thcm , jpi, jpj, npb(1:nbpb) )
	296	CALL tab_2d_1d_2( nbpb, qldif_1d (1:nbpb) , qldif , jpi, jpj, npb(1:nbpb) )
	297	CALL tab_2d_1d_2( nbpb, qstbif_1d (1:nbpb) , qstoif , jpi, jpj, npb(1:nbpb) )
[1858]	298	!!gm bug ?? rdm_snw, rdq_snw should be transformed into 1D fields ? like rdm_ice and rdq_ice ???
	299	!!gm I think that in fact the transformation is useless for both snow and ice as they are already set to zero
	300	!!gm at the begining of the routine and only the value over icy point is updated in 1D and converted into 2D
	301	!!gm the two lines below are thus useless.
	302	CALL tab_2d_1d_2( nbpb, rdm_ice_1d (1:nbpb) , rdm_ice , jpi, jpj, npb(1:nbpb) )
	303	CALL tab_2d_1d_2( nbpb, rdq_ice_1d (1:nbpb) , rdq_ice , jpi, jpj, npb(1:nbpb) )
	304	!!gm bug end.
[821]	305	CALL tab_2d_1d_2( nbpb, dmgwi_1d (1:nbpb) , dmgwi , jpi, jpj, npb(1:nbpb) )
	306	CALL tab_2d_1d_2( nbpb, qlbbq_1d (1:nbpb) , zqlbsbq , jpi, jpj, npb(1:nbpb) )
[1218]	307	!
[821]	308	CALL lim_thd_zdf_2( 1, nbpb ) ! compute ice growth
[1218]	309	!
[3]	310	! back to the geographic grid.
[821]	311	CALL tab_1d_2d_2( nbpb, frld , npb, frld_1d (1:nbpb) , jpi, jpj )
	312	CALL tab_1d_2d_2( nbpb, hicif , npb, h_ice_1d (1:nbpb) , jpi, jpj )
	313	CALL tab_1d_2d_2( nbpb, hsnif , npb, h_snow_1d (1:nbpb) , jpi, jpj )
	314	CALL tab_1d_2d_2( nbpb, sist , npb, sist_1d (1:nbpb) , jpi, jpj )
	315	CALL tab_1d_2d_2( nbpb, tbif(:,:,1), npb, tbif_1d (1:nbpb , 1 ), jpi, jpj )
	316	CALL tab_1d_2d_2( nbpb, tbif(:,:,2), npb, tbif_1d (1:nbpb , 2 ), jpi, jpj )
	317	CALL tab_1d_2d_2( nbpb, tbif(:,:,3), npb, tbif_1d (1:nbpb , 3 ), jpi, jpj )
	318	CALL tab_1d_2d_2( nbpb, fscmbq , npb, fscbq_1d (1:nbpb) , jpi, jpj )
	319	CALL tab_1d_2d_2( nbpb, ffltbif , npb, fltbif_1d (1:nbpb) , jpi, jpj )
	320	CALL tab_1d_2d_2( nbpb, fstric , npb, fstbif_1d (1:nbpb) , jpi, jpj )
	321	CALL tab_1d_2d_2( nbpb, qldif , npb, qldif_1d (1:nbpb) , jpi, jpj )
	322	CALL tab_1d_2d_2( nbpb, qfvbq , npb, qfvbq_1d (1:nbpb) , jpi, jpj )
	323	CALL tab_1d_2d_2( nbpb, qstoif , npb, qstbif_1d (1:nbpb) , jpi, jpj )
[1858]	324	CALL tab_1d_2d_2( nbpb, rdm_ice , npb, rdm_ice_1d(1:nbpb) , jpi, jpj )
	325	CALL tab_1d_2d_2( nbpb, rdq_ice , npb, rdq_ice_1d(1:nbpb) , jpi, jpj )
	326	CALL tab_1d_2d_2( nbpb, rdm_snw , npb, rdm_snw_1d(1:nbpb) , jpi, jpj )
	327	CALL tab_1d_2d_2( nbpb, rdq_snw , npb, rdq_snw_1d(1:nbpb) , jpi, jpj )
[821]	328	CALL tab_1d_2d_2( nbpb, dmgwi , npb, dmgwi_1d (1:nbpb) , jpi, jpj )
[1756]	329	CALL tab_1d_2d_2( nbpb, zdvosif , npb, dvsbq_1d (1:nbpb) , jpi, jpj )
	330	CALL tab_1d_2d_2( nbpb, zdvobif , npb, dvbbq_1d (1:nbpb) , jpi, jpj )
	331	CALL tab_1d_2d_2( nbpb, zdvomif , npb, rdvomif_1d(1:nbpb) , jpi, jpj )
	332	CALL tab_1d_2d_2( nbpb, zdvolif , npb, dvlbq_1d (1:nbpb) , jpi, jpj )
	333	CALL tab_1d_2d_2( nbpb, zdvonif , npb, dvnbq_1d (1:nbpb) , jpi, jpj )
[1482]	334	CALL tab_1d_2d_2( nbpb, qsr_ice(:,:,1), npb, qsr_ice_1d(1:nbpb) , jpi, jpj )
	335	CALL tab_1d_2d_2( nbpb, qns_ice(:,:,1), npb, qns_ice_1d(1:nbpb) , jpi, jpj )
	336	IF( .NOT. lk_cpl ) CALL tab_1d_2d_2( nbpb, qla_ice(:,:,1), npb, qla_ice_1d(1:nbpb) , jpi, jpj )
[1218]	337	!
[3]	338	ENDIF
	339
[1218]	340	! Up-date sea ice thickness
	341	!--------------------------
[3]	342	DO jj = 1, jpj
	343	DO ji = 1, jpi
	344	phicif(ji,jj) = hicif(ji,jj)
[1218]	345	hicif(ji,jj) = hicif(ji,jj) * ( rone - MAX( rzero, SIGN( rone, - ( 1.0 - frld(ji,jj) ) ) ) )
[3]	346	END DO
	347	END DO
	348
	349
[1218]	350	! Tricky trick : add 2 to frld in the Southern Hemisphere
	351	!--------------------------------------------------------
[192]	352	IF( fcor(1,1) < 0.e0 ) THEN
[421]	353	DO jj = 1, njeqm1
[192]	354	DO ji = 1, jpi
	355	frld(ji,jj) = frld(ji,jj) + 2.0
	356	END DO
[3]	357	END DO
[192]	358	ENDIF
[3]	359
	360
[1218]	361	! Select points for lateral accretion (this occurs when heat exchange
	362	! between ice and ocean is negative; ocean losing heat)
[3]	363	!-----------------------------------------------------------------
	364	nbpac = 0
	365	DO jj = 1, jpj
	366	DO ji = 1, jpi
[70]	367	!i yes! IF ( ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN
	368	IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN
[3]	369	nbpac = nbpac + 1
	370	npac( nbpac ) = (jj - 1) * jpi + ji
	371	ENDIF
	372	END DO
	373	END DO
	374
[258]	375	IF(ln_ctl) THEN
	376	CALL prt_ctl(tab2d_1=phicif, clinfo1=' lim_thd: phicif : ', tab2d_2=hicif, clinfo2=' hicif : ')
	377	WRITE(charout, FMT="('lim_thd: nbpac = ',I4)") nbpac
	378	CALL prt_ctl_info(charout)
[3]	379	ENDIF
	380
[1218]	381
	382	! If ocean gains heat do nothing ; otherwise, one performs lateral accretion
[3]	383	!--------------------------------------------------------------------------------
[70]	384	IF( nbpac > 0 ) THEN
[1218]	385	!
[1858]	386	zlicegr(:,:) = rdm_ice(:,:) ! to output the lateral sea-ice growth
[3]	387	!...Put the variable in a 1-D array for lateral accretion
[821]	388	CALL tab_2d_1d_2( nbpac, frld_1d (1:nbpac) , frld , jpi, jpj, npac(1:nbpac) )
	389	CALL tab_2d_1d_2( nbpac, h_snow_1d (1:nbpac) , hsnif , jpi, jpj, npac(1:nbpac) )
	390	CALL tab_2d_1d_2( nbpac, h_ice_1d (1:nbpac) , hicif , jpi, jpj, npac(1:nbpac) )
	391	CALL tab_2d_1d_2( nbpac, tbif_1d (1:nbpac , 1 ), tbif(:,:,1), jpi, jpj, npac(1:nbpac) )
	392	CALL tab_2d_1d_2( nbpac, tbif_1d (1:nbpac , 2 ), tbif(:,:,2), jpi, jpj, npac(1:nbpac) )
	393	CALL tab_2d_1d_2( nbpac, tbif_1d (1:nbpac , 3 ), tbif(:,:,3), jpi, jpj, npac(1:nbpac) )
	394	CALL tab_2d_1d_2( nbpac, qldif_1d (1:nbpac) , qldif , jpi, jpj, npac(1:nbpac) )
	395	CALL tab_2d_1d_2( nbpac, qcmif_1d (1:nbpac) , qcmif , jpi, jpj, npac(1:nbpac) )
	396	CALL tab_2d_1d_2( nbpac, qstbif_1d (1:nbpac) , qstoif , jpi, jpj, npac(1:nbpac) )
[1858]	397	CALL tab_2d_1d_2( nbpac, rdm_ice_1d(1:nbpac) , rdm_ice , jpi, jpj, npac(1:nbpac) )
	398	CALL tab_2d_1d_2( nbpac, rdq_ice_1d(1:nbpac) , rdq_ice , jpi, jpj, npac(1:nbpac) )
[1756]	399	CALL tab_2d_1d_2( nbpac, dvlbq_1d (1:nbpac) , zdvolif , jpi, jpj, npac(1:nbpac) )
[821]	400	CALL tab_2d_1d_2( nbpac, tfu_1d (1:nbpac) , tfu , jpi, jpj, npac(1:nbpac) )
[1218]	401	!
	402	CALL lim_thd_lac_2( 1 , nbpac ) ! lateral accretion routine.
	403	!
[70]	404	! back to the geographic grid
[821]	405	CALL tab_1d_2d_2( nbpac, frld , npac(1:nbpac), frld_1d (1:nbpac) , jpi, jpj )
	406	CALL tab_1d_2d_2( nbpac, hsnif , npac(1:nbpac), h_snow_1d (1:nbpac) , jpi, jpj )
	407	CALL tab_1d_2d_2( nbpac, hicif , npac(1:nbpac), h_ice_1d (1:nbpac) , jpi, jpj )
	408	CALL tab_1d_2d_2( nbpac, tbif(:,:,1), npac(1:nbpac), tbif_1d (1:nbpac , 1 ), jpi, jpj )
	409	CALL tab_1d_2d_2( nbpac, tbif(:,:,2), npac(1:nbpac), tbif_1d (1:nbpac , 2 ), jpi, jpj )
	410	CALL tab_1d_2d_2( nbpac, tbif(:,:,3), npac(1:nbpac), tbif_1d (1:nbpac , 3 ), jpi, jpj )
	411	CALL tab_1d_2d_2( nbpac, qstoif , npac(1:nbpac), qstbif_1d (1:nbpac) , jpi, jpj )
[1858]	412	CALL tab_1d_2d_2( nbpac, rdm_ice , npac(1:nbpac), rdm_ice_1d(1:nbpac) , jpi, jpj )
	413	CALL tab_1d_2d_2( nbpac, rdq_ice , npac(1:nbpac), rdq_ice_1d(1:nbpac) , jpi, jpj )
[1756]	414	CALL tab_1d_2d_2( nbpac, zdvolif , npac(1:nbpac), dvlbq_1d (1:nbpac) , jpi, jpj )
[1218]	415	!
[70]	416	ENDIF
[3]	417
	418
[1218]	419	! Recover frld values between 0 and 1 in the Southern Hemisphere (tricky trick)
	420	! Update daily thermodynamic ice production.
[70]	421	!------------------------------------------------------------------------------
[3]	422	DO jj = 1, jpj
	423	DO ji = 1, jpi
	424	frld (ji,jj) = MIN( frld(ji,jj), ABS( frld(ji,jj) - 2.0 ) )
[1482]	425	fr_i (ji,jj) = 1.0 - frld(ji,jj)
	426	hicifp(ji,jj) = hicif(ji,jj) * fr_i(ji,jj) - hicifp(ji,jj)
[3]	427	END DO
	428	END DO
	429
[1482]	430	! Outputs
	431	!--------------------------------------------------------------------------------
[1756]	432	ztmp(:,:) = 1. - pfrld(:,:) ! fraction of ice after the dynamic, before the thermodynamic
	433	CALL iom_put( 'ioceflxb', fbif ) ! Oceanic flux at the ice base [W/m2 ???]
	434	CALL iom_put( 'ist_cea', (sist(:,:) - rt0) * ztmp(:,:) ) ! Ice surface temperature [Celius]
	435	CALL iom_put( 'qsr_ai_cea', qsr_ice(:,:,1) * ztmp(:,:) ) ! Solar flux over the ice [W/m2]
	436	CALL iom_put( 'qns_ai_cea', qns_ice(:,:,1) * ztmp(:,:) ) ! Non-solar flux over the ice [W/m2]
[1482]	437	IF( .NOT. lk_cpl ) CALL iom_put( 'qla_ai_cea', qla_ice(:,:,1) * ztmp(:,:) ) ! Latent flux over the ice [W/m2]
	438	!
[1756]	439	CALL iom_put( 'snowthic_cea', hsnif (:,:) * fr_i(:,:) ) ! Snow thickness [m]
	440	CALL iom_put( 'icethic_cea' , hicif (:,:) * fr_i(:,:) ) ! Ice thickness [m]
	441	zztmp = 1.0 / rdt_ice
	442	CALL iom_put( 'iceprod_cea' , hicifp (:,:) * zztmp ) ! Ice produced [m/s]
	443	IF( lk_diaar5 ) THEN
[1858]	444	CALL iom_put( 'snowmel_cea' , rdm_snw(:,:) * zztmp ) ! Snow melt [kg/m2/s]
[1756]	445	zztmp = rhoic / rdt_ice
	446	CALL iom_put( 'sntoice_cea' , zdvonif(:,:) * zztmp ) ! Snow to Ice transformation [kg/m2/s]
	447	CALL iom_put( 'ticemel_cea' , zdvosif(:,:) * zztmp ) ! Melt at Sea Ice top [kg/m2/s]
	448	CALL iom_put( 'bicemel_cea' , zdvomif(:,:) * zztmp ) ! Melt at Sea Ice bottom [kg/m2/s]
[1858]	449	zlicegr(:,:) = MAX( 0.e0, rdm_ice(:,:)-zlicegr(:,:) )
[1756]	450	CALL iom_put( 'licepro_cea' , zlicegr(:,:) * zztmp ) ! Latereal sea ice growth [kg/m2/s]
	451	ENDIF
[1482]	452	!
[1756]	453	! Compute the Eastward & Northward sea-ice transport
	454	zztmp = 0.25 * rhoic
	455	DO jj = 1, jpjm1
	456	DO ji = 1, jpim1 ! NO vector opt.
	457	! Ice velocities, volume & transport at U & V-points
	458	zuice_m = u_ice(ji+1,jj+1) + u_ice(ji+1,jj )
	459	zvice_m = v_ice(ji+1,jj+1) + v_ice(ji ,jj+1)
	460	zhice_u = hicif(ji,jj)e2t(ji,jj)fr_i(ji,jj) + hicif(ji+1,jj )e2t(ji+1,jj )fr_i(ji+1,jj )
	461	zhice_v = hicif(ji,jj)e1t(ji,jj)fr_i(ji,jj) + hicif(ji ,jj+1)e1t(ji ,jj+1)fr_i(ji ,jj+1)
	462	zu_imasstr(ji,jj) = zztmp * zhice_u * zuice_m
	463	zv_imasstr(ji,jj) = zztmp * zhice_v * zvice_m
	464	END DO
	465	END DO
	466	CALL lbc_lnk( zu_imasstr, 'U', -1. ) ; CALL lbc_lnk( zv_imasstr, 'V', -1. )
	467	CALL iom_put( 'u_imasstr', zu_imasstr(:,:) ) ! Ice transport along i-axis at U-point [kg/s]
	468	CALL iom_put( 'v_imasstr', zv_imasstr(:,:) ) ! Ice transport along j-axis at V-point [kg/s]
[1482]	469
[1756]	470	!! Fram Strait sea-ice transport (sea-ice + snow) (in ORCA2 = 5 points)
	471	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
	472	DO jj = mj0(137), mj1(137) ! B grid
	473	IF( mj0(jj-1) >= nldj ) THEN
	474	DO ji = MAX(mi0(134),nldi), MIN(mi1(138),nlei)
	475	zrhoij = e1t(ji,jj ) * fr_i(ji,jj ) * ( rhoichicif(ji,jj ) + rhosnhsnif(ji,jj ) )
	476	zrhoijm1 = e1t(ji,jj-1) * fr_i(ji,jj-1) * ( rhoichicif(ji,jj-1) + rhosnhsnif(ji,jj-1) )
	477	ztr_fram = ztr_fram - 0.25 * ( v_ice(ji,jj)+ v_ice(ji+1,jj) ) * ( zrhoij + zrhoijm1 )
	478	END DO
	479	ENDIF
	480	END DO
	481	IF( lk_mpp ) CALL mpp_sum( ztr_fram )
	482	CALL iom_put( 'fram_trans', ztr_fram ) ! Ice transport through Fram strait [kg/s]
	483	ENDIF
	484
	485	ztmp(:,:) = 1. - AINT( frld(:,:), wp ) ! return 1 as soon as there is ice
	486	CALL iom_put( 'ice_pres' , ztmp ) ! Ice presence [-]
	487	CALL iom_put( 'ist_ipa' , ( sist(:,:) - rt0 ) * ztmp(:,:) ) ! Ice surface temperature [Celius]
	488	CALL iom_put( 'uice_ipa' , u_ice(:,:) * ztmp(:,:) ) ! Ice velocity along i-axis at I-point [m/s]
	489	CALL iom_put( 'vice_ipa' , v_ice(:,:) * ztmp(:,:) ) ! Ice velocity along j-axis at I-point [m/s]
	490
[258]	491	IF(ln_ctl) THEN
	492	CALL prt_ctl_info(' lim_thd end ')
[1218]	493	CALL prt_ctl( tab2d_1=hicif , clinfo1=' lim_thd: hicif : ', tab2d_2=hsnif , clinfo2=' hsnif : ' )
	494	CALL prt_ctl( tab2d_1=frld , clinfo1=' lim_thd: frld : ', tab2d_2=hicifp, clinfo2=' hicifp : ' )
	495	CALL prt_ctl( tab2d_1=phicif , clinfo1=' lim_thd: phicif : ', tab2d_2=pfrld , clinfo2=' pfrld : ' )
	496	CALL prt_ctl( tab2d_1=sist , clinfo1=' lim_thd: sist : ' )
	497	CALL prt_ctl( tab2d_1=tbif(:,:,1), clinfo1=' lim_thd: tbif 1 : ' )
	498	CALL prt_ctl( tab2d_1=tbif(:,:,2), clinfo1=' lim_thd: tbif 2 : ' )
	499	CALL prt_ctl( tab2d_1=tbif(:,:,3), clinfo1=' lim_thd: tbif 3 : ' )
	500	CALL prt_ctl( tab2d_1=fdtcn , clinfo1=' lim_thd: fdtcn : ', tab2d_2=qdtcn , clinfo2=' qdtcn : ' )
	501	CALL prt_ctl( tab2d_1=qstoif , clinfo1=' lim_thd: qstoif : ', tab2d_2=fsbbq , clinfo2=' fsbbq : ' )
[3]	502	ENDIF
[888]	503	!
[821]	504	END SUBROUTINE lim_thd_2
[3]	505
[719]	506
[821]	507	SUBROUTINE lim_thd_init_2
[3]	508	!!-------------------------------------------------------------------
[821]	509	!! * ROUTINE lim_thd_init_2 *
[3]	510	!!
	511	!! ** Purpose : Physical constants and parameters linked to the ice
	512	!! thermodynamics
	513	!!
	514	!! ** Method : Read the namicethd namelist and check the ice-thermo
	515	!! parameter values called at the first timestep (nit000)
	516	!!
	517	!! ** input : Namelist namicether
	518	!!-------------------------------------------------------------------
	519	NAMELIST/namicethd/ hmelt , hiccrit, hicmin, hiclim, amax , &
	520	& swiqst, sbeta , parlat, hakspl, hibspl, exld, &
	521	& hakdif, hnzst , thth , parsub, alphs
	522	!!-------------------------------------------------------------------
[1218]	523	!
	524	REWIND( numnam_ice ) ! read namelist
[3]	525	READ ( numnam_ice , namicethd )
[1218]	526	IF( lk_cpl .AND. parsub /= 0.0 ) CALL ctl_stop( 'In coupled mode, use parsub = 0. or send dqla' )
	527	!
	528	IF(lwp) THEN ! control print
[3]	529	WRITE(numout,*)
[821]	530	WRITE(numout,*)'lim_thd_init_2: ice parameters for ice thermodynamic computation '
	531	WRITE(numout,*)'~~~~~~~~~~~~~~'
[70]	532	WRITE(numout,*)' maximum melting at the bottom hmelt = ', hmelt
	533	WRITE(numout,*)' ice thick. for lateral accretion in NH (SH) hiccrit(1/2) = ', hiccrit
	534	WRITE(numout,*)' ice thick. corr. to max. energy stored in brine pocket hicmin = ', hicmin
	535	WRITE(numout,*)' minimum ice thickness hiclim = ', hiclim
	536	WRITE(numout,*)' maximum lead fraction amax = ', amax
[1218]	537	WRITE(numout,*)' energy stored in brine pocket (=1) or not (=0) swiqst = ', swiqst
[70]	538	WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice '
	539	WRITE(numout,*)' Cranck-Nicholson (=0.5), implicit (=1), explicit (=0) sbeta = ', sbeta
	540	WRITE(numout,*)' percentage of energy used for lateral ablation parlat = ', parlat
	541	WRITE(numout,*)' slope of distr. for Hakkinen-Mellor lateral melting hakspl = ', hakspl
	542	WRITE(numout,*)' slope of distribution for Hibler lateral melting hibspl = ', hibspl
	543	WRITE(numout,*)' exponent for leads-closure rate exld = ', exld
	544	WRITE(numout,*)' coefficient for diffusions of ice and snow hakdif = ', hakdif
	545	WRITE(numout,*)' threshold thick. for comp. of eq. thermal conductivity zhth = ', thth
	546	WRITE(numout,*)' thickness of the surf. layer in temp. computation hnzst = ', hnzst
	547	WRITE(numout,*)' switch for snow sublimation (=1) or not (=0) parsub = ', parsub
	548	WRITE(numout,*)' coefficient for snow density when snow ice formation alphs = ', alphs
[3]	549	ENDIF
[1218]	550	!
[3]	551	uscomi = 1.0 / ( 1.0 - amax ) ! inverse of minimum lead fraction
	552	rcdsn = hakdif * rcdsn
	553	rcdic = hakdif * rcdic
[1218]	554	!
	555	IF( hsndif > 100.e0 .OR. hicdif > 100.e0 ) THEN
[3]	556	cnscg = 0.e0
	557	ELSE
	558	cnscg = rcpsn / rcpic ! ratio rcpsn/rcpic
	559	ENDIF
[1218]	560	!
[821]	561	END SUBROUTINE lim_thd_init_2
[3]	562
	563	#else
[70]	564	!!----------------------------------------------------------------------
[821]	565	!! Default option Dummy module NO LIM 2.0 sea-ice model
[70]	566	!!----------------------------------------------------------------------
[1857]	567	CONTAINS
	568	SUBROUTINE lim_thd_2 ! Dummy routine
	569	END SUBROUTINE lim_thd_2
[3]	570	#endif
	571
	572	!!======================================================================
[821]	573	END MODULE limthd_2

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

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

Context Navigation

source: branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_2.F90 @ 1858

Download in other formats: