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_2.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90 @ 3584

Last change on this file since 3584 was 3558, checked in by rblod, 12 years ago
Fix issues when using key_nosignedzeo, see ticket #996
Property svn:keywords set to `Id`
File size: 33.7 KB

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

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

Download in other formats: