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 @ 3260

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

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

Download in other formats: