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.

iceistate.F90 in NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/iceistate.F90 @ 11413

Last change on this file since 11413 was 11413, checked in by gsamson, 5 years ago

dev_r11265_ABL : see #2131

merge src and cfgs from HPC-13_IRRMANN_BDY_optimization branch @ r11402 with dev_r11265_ABL branch @ r11363
change ORCA2 results due to ice rheology "cleaning" (see commit r11377)

Property svn:keywords set to Id

File size: 25.8 KB

Rev	Line
[8586]	1	MODULE iceistate
	2	!!======================================================================
	3	!! * MODULE iceistate *
	4	!! sea-ice : Initialization of ice variables
	5	!!======================================================================
[9604]	6	!! History : 2.0 ! 2004-01 (C. Ethe, G. Madec) Original code
	7	!! 3.0 ! 2007 (M. Vancoppenolle) Rewrite for ice cats
	8	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
[8586]	9	!!----------------------------------------------------------------------
[9570]	10	#if defined key_si3
[8586]	11	!!----------------------------------------------------------------------
[9570]	12	!! 'key_si3' SI3 sea-ice model
[8586]	13	!!----------------------------------------------------------------------
	14	!! ice_istate : initialization of diagnostics ice variables
	15	!! ice_istate_init : initialization of ice state and namelist read
	16	!!----------------------------------------------------------------------
	17	USE phycst ! physical constant
	18	USE oce ! dynamics and tracers variables
	19	USE dom_oce ! ocean domain
	20	USE sbc_oce , ONLY : sst_m, sss_m, ln_ice_embd
	21	USE sbc_ice , ONLY : tn_ice, snwice_mass, snwice_mass_b
	22	USE eosbn2 ! equation of state
	23	USE domvvl ! Variable volume
[11413]	24	USE ice ! sea-ice: variables
	25	USE ice1D ! sea-ice: thermodynamics variables
	26	USE icetab ! sea-ice: 1D <==> 2D transformation
	27	USE icevar ! sea-ice: operations
[8586]	28	!
	29	USE in_out_manager ! I/O manager
	30	USE iom ! I/O manager library
	31	USE lib_mpp ! MPP library
	32	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
	33	USE fldread ! read input fields
	34
	35	IMPLICIT NONE
	36	PRIVATE
	37
	38	PUBLIC ice_istate ! called by icestp.F90
	39	PUBLIC ice_istate_init ! called by icestp.F90
	40	!
[8813]	41	! !! namelist (namini)
[11413]	42	LOGICAL, PUBLIC :: ln_iceini !: Ice initialization or not
	43	LOGICAL, PUBLIC :: ln_iceini_file !: Ice initialization from 2D netcdf file
	44	REAL(wp) :: rn_thres_sst
	45	REAL(wp) :: rn_hti_ini_n, rn_hts_ini_n, rn_ati_ini_n, rn_smi_ini_n, rn_tmi_ini_n, rn_tsu_ini_n, rn_tms_ini_n
	46	REAL(wp) :: rn_hti_ini_s, rn_hts_ini_s, rn_ati_ini_s, rn_smi_ini_s, rn_tmi_ini_s, rn_tsu_ini_s, rn_tms_ini_s
	47	REAL(wp) :: rn_apd_ini_n, rn_hpd_ini_n
	48	REAL(wp) :: rn_apd_ini_s, rn_hpd_ini_s
	49	!
	50	! ! if ln_iceini_file = T
	51	INTEGER , PARAMETER :: jpfldi = 9 ! maximum number of files to read
	52	INTEGER , PARAMETER :: jp_hti = 1 ! index of ice thickness (m)
	53	INTEGER , PARAMETER :: jp_hts = 2 ! index of snw thickness (m)
	54	INTEGER , PARAMETER :: jp_ati = 3 ! index of ice fraction (-)
	55	INTEGER , PARAMETER :: jp_smi = 4 ! index of ice salinity (g/kg)
	56	INTEGER , PARAMETER :: jp_tmi = 5 ! index of ice temperature (K)
	57	INTEGER , PARAMETER :: jp_tsu = 6 ! index of ice surface temp (K)
	58	INTEGER , PARAMETER :: jp_tms = 7 ! index of snw temperature (K)
	59	INTEGER , PARAMETER :: jp_apd = 8 ! index of pnd fraction (-)
	60	INTEGER , PARAMETER :: jp_hpd = 9 ! index of pnd depth (m)
	61	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: si ! structure of input fields (file informations, fields read)
	62	!
[8586]	63	!!----------------------------------------------------------------------
[9598]	64	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
[10069]	65	!! $Id$
[10413]	66	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[8586]	67	!!----------------------------------------------------------------------
	68	CONTAINS
	69
[11413]	70	SUBROUTINE ice_istate( kt )
[8586]	71	!!-------------------------------------------------------------------
	72	!! * ROUTINE ice_istate *
	73	!!
	74	!! ** Purpose : defined the sea-ice initial state
	75	!!
	76	!! ** Method : This routine will put some ice where ocean
	77	!! is at the freezing point, then fill in ice
	78	!! state variables using prescribed initial
	79	!! values in the namelist
	80	!!
	81	!! ** Steps : 1) Set initial surface and basal temperatures
	82	!! 2) Recompute or read sea ice state variables
	83	!! 3) Fill in the ice thickness distribution using gaussian
	84	!! 4) Fill in space-dependent arrays for state variables
	85	!! 5) snow-ice mass computation
	86	!! 6) store before fields
	87	!!
	88	!! ** Notes : o_i, t_su, t_s, t_i, sz_i must be filled everywhere, even
[11413]	89	!! where there is no ice
[8586]	90	!!--------------------------------------------------------------------
[11413]	91	INTEGER, INTENT(in) :: kt ! time step
	92	!!
[8813]	93	INTEGER :: ji, jj, jk, jl ! dummy loop indices
[11413]	94	REAL(wp) :: ztmelts
[8586]	95	INTEGER , DIMENSION(4) :: itest
	96	REAL(wp), DIMENSION(jpi,jpj) :: z2d
	97	REAL(wp), DIMENSION(jpi,jpj) :: zswitch ! ice indicator
[11413]	98	REAL(wp), DIMENSION(jpi,jpj) :: zht_i_ini, zat_i_ini, ztm_s_ini !data from namelist or nc file
	99	REAL(wp), DIMENSION(jpi,jpj) :: zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
	100	REAL(wp), DIMENSION(jpi,jpj) :: zapnd_ini, zhpnd_ini !data from namelist or nc file
	101	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zti_3d , zts_3d !temporary arrays
	102	!!
	103	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d
[8586]	104	!--------------------------------------------------------------------
	105
	106	IF(lwp) WRITE(numout,*)
	107	IF(lwp) WRITE(numout,*) 'ice_istate: sea-ice initialization '
	108	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
	109
[11413]	110	!---------------------------
	111	! 1) 1st init. of the fields
	112	!---------------------------
[8586]	113	!
[11413]	114	! basal temperature (considered at freezing point) [Kelvin]
	115	CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
	116	t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1)
	117	!
	118	! surface temperature and conductivity
[8586]	119	DO jl = 1, jpl
[9416]	120	t_su (:,:,jl) = rt0 * tmask(:,:,1) ! temp at the surface
[10534]	121	cnd_ice(:,:,jl) = 0._wp ! initialisation of the effective conductivity at the top of ice/snow (ln_cndflx=T)
[8586]	122	END DO
[8813]	123	!
[11413]	124	! ice and snw temperatures
	125	DO jl = 1, jpl
	126	DO jk = 1, nlay_i
	127	t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
	128	END DO
	129	DO jk = 1, nlay_s
	130	t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
	131	END DO
	132	END DO
	133	!
	134	! specific temperatures for coupled runs
	135	tn_ice (:,:,:) = t_i (:,:,1,:)
	136	t1_ice (:,:,:) = t_i (:,:,1,:)
[8586]	137
[11413]	138	! heat contents
	139	e_i (:,:,:,:) = 0._wp
	140	e_s (:,:,:,:) = 0._wp
	141
	142	! general fields
	143	a_i (:,:,:) = 0._wp
	144	v_i (:,:,:) = 0._wp
	145	v_s (:,:,:) = 0._wp
	146	sv_i(:,:,:) = 0._wp
	147	oa_i(:,:,:) = 0._wp
	148	!
	149	h_i (:,:,:) = 0._wp
	150	h_s (:,:,:) = 0._wp
	151	s_i (:,:,:) = 0._wp
	152	o_i (:,:,:) = 0._wp
	153	!
	154	! melt ponds
	155	a_ip (:,:,:) = 0._wp
	156	v_ip (:,:,:) = 0._wp
	157	a_ip_frac(:,:,:) = 0._wp
	158	h_ip (:,:,:) = 0._wp
	159	!
	160	! ice velocities
	161	u_ice (:,:) = 0._wp
	162	v_ice (:,:) = 0._wp
	163	!
	164	!------------------------------------------------------------------------
	165	! 2) overwrite some of the fields with namelist parameters or netcdf file
	166	!------------------------------------------------------------------------
[8586]	167	IF( ln_iceini ) THEN
	168	! !---------------!
	169	IF( ln_iceini_file )THEN ! Read a file !
	170	! !---------------!
[11413]	171	WHERE( ff_t(:,:) >= 0._wp ) ; zswitch(:,:) = 1._wp
	172	ELSEWHERE ; zswitch(:,:) = 0._wp
	173	END WHERE
[8586]	174	!
[11413]	175	CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step
[8586]	176	!
[11413]	177	! -- mandatory fields -- !
	178	zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1)
	179	zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1)
	180	zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1)
	181
	182	! -- optional fields -- !
	183	! if fields do not exist then set them to the values present in the namelist (except for snow and surface temperature)
	184	!
	185	! ice salinity
	186	IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) &
	187	& si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
	188	zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1)
	189	!
	190	! ice temperature
	191	IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) &
	192	& si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
	193	ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1)
	194	!
	195	! surface temperature => set to ice temperature if it exists
	196	IF ( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) THEN
	197	si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
	198	ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN
	199	si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
	200	ENDIF
	201	zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1)
	202	!
	203	! snow temperature => set to ice temperature if it exists
	204	IF ( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) THEN
	205	si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
	206	ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN
	207	si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
	208	ENDIF
	209	ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1)
	210	!
	211	! ponds
	212	IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) &
	213	& si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
	214	zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1)
	215	IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) &
	216	& si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
	217	zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1)
	218	!
	219	! change the switch for the following
	220	WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1)
	221	ELSEWHERE ; zswitch(:,:) = 0._wp
[8586]	222	END WHERE
[8813]	223	! !---------------!
[8586]	224	ELSE ! Read namelist !
	225	! !---------------!
[11413]	226	! no ice if (sst - Tfreez) >= thresold
[8813]	227	WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp
	228	ELSEWHERE ; zswitch(:,:) = tmask(:,:,1)
[8586]	229	END WHERE
[8813]	230	!
[8586]	231	! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
	232	WHERE( ff_t(:,:) >= 0._wp )
	233	zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:)
	234	zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:)
	235	zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:)
	236	zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:)
	237	ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
[11413]	238	zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:)
	239	ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:)
	240	zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:)
	241	zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:)
[8586]	242	ELSEWHERE
	243	zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:)
	244	zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:)
	245	zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:)
	246	zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:)
	247	ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
[11413]	248	zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:)
	249	ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:)
	250	zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:)
	251	zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:)
[8586]	252	END WHERE
	253	!
	254	ENDIF
[11413]	255	!-------------!
	256	! fill fields !
	257	!-------------!
	258	! select ice covered grid points
	259	npti = 0 ; nptidx(:) = 0
	260	DO jj = 1, jpj
	261	DO ji = 1, jpi
	262	IF ( zht_i_ini(ji,jj) > 0._wp ) THEN
	263	npti = npti + 1
	264	nptidx(npti) = (jj - 1) * jpi + ji
	265	ENDIF
	266	END DO
	267	END DO
	268
	269	! move to 1D arrays: (jpi,jpj) -> (jpi*jpj)
	270	CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti) , zht_i_ini )
	271	CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti) , zht_s_ini )
	272	CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti) , zat_i_ini )
	273	CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini )
	274	CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini )
	275	CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti) , zt_su_ini )
	276	CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti) , zsm_i_ini )
	277
	278	! allocate temporary arrays
	279	ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), &
	280	& zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl) )
[8586]	281
[11413]	282	! distribute 1-cat into jpl-cat: (jpijpj) -> (jpijpj,jpl)
	283	CALL ice_var_itd( h_i_1d(1:npti) , h_s_1d(1:npti) , at_i_1d(1:npti), zhi_2d, zhs_2d, zai_2d , &
	284	& t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), s_i_1d(1:npti), zti_2d, zts_2d, ztsu_2d, zsi_2d )
[8586]	285
[11413]	286	! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl)
	287	DO jl = 1, jpl
	288	zti_3d(:,:,jl) = rt0 * tmask(:,:,1)
	289	zts_3d(:,:,jl) = rt0 * tmask(:,:,1)
	290	END DO
	291	CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d , h_i )
	292	CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d , h_s )
	293	CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d , a_i )
	294	CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d , zti_3d )
	295	CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d , zts_3d )
	296	CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d , t_su )
	297	CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d , s_i )
[10413]	298
[11413]	299	! deallocate temporary arrays
	300	DEALLOCATE( zhi_2d, zhs_2d, zai_2d , &
	301	& zti_2d, zts_2d, ztsu_2d, zsi_2d )
[8586]	302
[11413]	303	! Melt ponds: distribute uniformely over the categories
	304	IF ( ln_pnd_CST .OR. ln_pnd_H12 ) THEN
	305	DO jl = 1, jpl
	306	a_ip_frac(:,:,jl) = zapnd_ini(:,:)
	307	h_ip (:,:,jl) = zhpnd_ini(:,:)
	308	a_ip (:,:,jl) = a_ip_frac(:,:,jl) * a_i (:,:,jl)
	309	v_ip (:,:,jl) = h_ip (:,:,jl) * a_ip(:,:,jl)
[10413]	310	END DO
	311	ENDIF
[11413]	312
	313	! calculate extensive and intensive variables
	314	CALL ice_var_salprof ! for sz_i
	315	DO jl = 1, jpl
[8586]	316	DO jj = 1, jpj
	317	DO ji = 1, jpi
[11413]	318	v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)
	319	v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)
	320	sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl)
[8586]	321	END DO
	322	END DO
	323	END DO
[8813]	324	!
[11413]	325	DO jl = 1, jpl
	326	DO jk = 1, nlay_s
[8586]	327	DO jj = 1, jpj
	328	DO ji = 1, jpi
[11413]	329	t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl)
	330	e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * &
	331	& rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus )
[8586]	332	END DO
	333	END DO
	334	END DO
	335	END DO
[8813]	336	!
[11413]	337	DO jl = 1, jpl
	338	DO jk = 1, nlay_i
[8586]	339	DO jj = 1, jpj
	340	DO ji = 1, jpi
[11413]	341	t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl)
	342	ztmelts = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K
	343	e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * &
	344	& rhoi * ( rcpi * ( ztmelts - t_i(ji,jj,jk,jl) ) + &
	345	& rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) &
	346	& - rcp * ( ztmelts - rt0 ) )
[8586]	347	END DO
	348	END DO
	349	END DO
	350	END DO
[8933]	351
[11413]	352	! specific temperatures for coupled runs
	353	tn_ice(:,:,:) = t_su(:,:,:)
	354	t1_ice(:,:,:) = t_i (:,:,1,:)
[8813]	355	!
[8586]	356	ENDIF ! ln_iceini
[8813]	357	!
[11413]	358	at_i(:,:) = SUM( a_i, dim=3 )
[8586]	359	!
	360	!----------------------------------------------
[11413]	361	! 3) Snow-ice mass (case ice is fully embedded)
[8586]	362	!----------------------------------------------
[9935]	363	snwice_mass (:,:) = tmask(:,:,1) * SUM( rhos * v_s(:,:,:) + rhoi * v_i(:,:,:), dim=3 ) ! snow+ice mass
[8586]	364	snwice_mass_b(:,:) = snwice_mass(:,:)
	365	!
	366	IF( ln_ice_embd ) THEN ! embedded sea-ice: deplete the initial ssh below sea-ice area
[8813]	367	!
[8586]	368	sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
	369	sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
[8813]	370	!
[8586]	371	IF( .NOT.ln_linssh ) THEN
[8813]	372	!
[8586]	373	WHERE( ht_0(:,:) > 0 ) ; z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:)
	374	ELSEWHERE ; z2d(:,:) = 1._wp ; END WHERE
[8813]	375	!
[8586]	376	DO jk = 1,jpkm1 ! adjust initial vertical scale factors
	377	e3t_n(:,:,jk) = e3t_0(:,:,jk) * z2d(:,:)
	378	e3t_b(:,:,jk) = e3t_n(:,:,jk)
	379	e3t_a(:,:,jk) = e3t_n(:,:,jk)
	380	END DO
[8813]	381	!
[8586]	382	! Reconstruction of all vertical scale factors at now and before time-steps
	383	! =========================================================================
	384	! Horizontal scale factor interpolations
	385	! --------------------------------------
	386	CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
	387	CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
	388	CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
	389	CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
	390	CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
	391	! Vertical scale factor interpolations
	392	! ------------------------------------
	393	CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' )
	394	CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
	395	CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
	396	CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
	397	CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
	398	! t- and w- points depth
	399	! ----------------------
	400	!!gm not sure of that....
	401	gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
	402	gdepw_n(:,:,1) = 0.0_wp
	403	gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
	404	DO jk = 2, jpk
	405	gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk )
	406	gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
	407	gde3w_n(:,:,jk) = gdept_n(:,:,jk ) - sshn (:,:)
	408	END DO
	409	ENDIF
	410	ENDIF
	411
	412	!------------------------------------
[11413]	413	! 4) store fields at before time-step
[8586]	414	!------------------------------------
	415	! it is only necessary for the 1st interpolation by Agrif
	416	a_i_b (:,:,:) = a_i (:,:,:)
	417	e_i_b (:,:,:,:) = e_i (:,:,:,:)
	418	v_i_b (:,:,:) = v_i (:,:,:)
	419	v_s_b (:,:,:) = v_s (:,:,:)
	420	e_s_b (:,:,:,:) = e_s (:,:,:,:)
	421	sv_i_b (:,:,:) = sv_i (:,:,:)
	422	oa_i_b (:,:,:) = oa_i (:,:,:)
	423	u_ice_b(:,:) = u_ice(:,:)
	424	v_ice_b(:,:) = v_ice(:,:)
[10527]	425	! total concentration is needed for Lupkes parameterizations
	426	at_i_b (:,:) = at_i (:,:)
[8586]	427
[8885]	428	!!clem: output of initial state should be written here but it is impossible because
	429	!! the ocean and ice are in the same file
[10425]	430	!! CALL dia_wri_state( 'output.init' )
[8813]	431	!
[8586]	432	END SUBROUTINE ice_istate
	433
[8813]	434
[8586]	435	SUBROUTINE ice_istate_init
	436	!!-------------------------------------------------------------------
	437	!! * ROUTINE ice_istate_init *
	438	!!
[8813]	439	!! ** Purpose : Definition of initial state of the ice
[8586]	440	!!
[8813]	441	!! ** Method : Read the namini namelist and check the parameter
	442	!! values called at the first timestep (nit000)
[8586]	443	!!
[8813]	444	!! ** input : Namelist namini
[8586]	445	!!
	446	!!-----------------------------------------------------------------------------
[11263]	447	INTEGER :: ios ! Local integer output status for namelist read
[8813]	448	INTEGER :: ifpr, ierror
[8586]	449	!
	450	CHARACTER(len=256) :: cn_dir ! Root directory for location of ice files
[11413]	451	TYPE(FLD_N) :: sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd
[8586]	452	TYPE(FLD_N), DIMENSION(jpfldi) :: slf_i ! array of namelist informations on the fields to read
	453	!
[11413]	454	NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, &
	455	& rn_hti_ini_n, rn_hti_ini_s, rn_hts_ini_n, rn_hts_ini_s, &
	456	& rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, &
	457	& rn_tmi_ini_n, rn_tmi_ini_s, rn_tsu_ini_n, rn_tsu_ini_s, rn_tms_ini_n, rn_tms_ini_s, &
	458	& rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, &
	459	& sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, cn_dir
[8586]	460	!!-----------------------------------------------------------------------------
	461	!
	462	REWIND( numnam_ice_ref ) ! Namelist namini in reference namelist : Ice initial state
	463	READ ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901)
[11348]	464	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in reference namelist' )
[8586]	465	REWIND( numnam_ice_cfg ) ! Namelist namini in configuration namelist : Ice initial state
	466	READ ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 )
[11348]	467	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namini in configuration namelist' )
[8586]	468	IF(lwm) WRITE ( numoni, namini )
[8813]	469	!
[8586]	470	slf_i(jp_hti) = sn_hti ; slf_i(jp_hts) = sn_hts
[11413]	471	slf_i(jp_ati) = sn_ati ; slf_i(jp_smi) = sn_smi
	472	slf_i(jp_tmi) = sn_tmi ; slf_i(jp_tsu) = sn_tsu ; slf_i(jp_tms) = sn_tms
	473	slf_i(jp_apd) = sn_apd ; slf_i(jp_hpd) = sn_hpd
[8586]	474	!
	475	IF(lwp) THEN ! control print
	476	WRITE(numout,*)
	477	WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation '
	478	WRITE(numout,*) '~~~~~~~~~~~~~~~'
	479	WRITE(numout,*) ' Namelist namini:'
[11413]	480	WRITE(numout,*) ' ice initialization (T) or not (F) ln_iceini = ', ln_iceini
	481	WRITE(numout,*) ' ice initialization from a netcdf file ln_iceini_file = ', ln_iceini_file
	482	WRITE(numout,*) ' max ocean temp. above Tfreeze with initial ice rn_thres_sst = ', rn_thres_sst
	483	IF( ln_iceini .AND. .NOT.ln_iceini_file ) THEN
	484	WRITE(numout,*) ' initial snw thickness in the north-south rn_hts_ini = ', rn_hts_ini_n,rn_hts_ini_s
	485	WRITE(numout,*) ' initial ice thickness in the north-south rn_hti_ini = ', rn_hti_ini_n,rn_hti_ini_s
	486	WRITE(numout,*) ' initial ice concentr in the north-south rn_ati_ini = ', rn_ati_ini_n,rn_ati_ini_s
	487	WRITE(numout,*) ' initial ice salinity in the north-south rn_smi_ini = ', rn_smi_ini_n,rn_smi_ini_s
	488	WRITE(numout,*) ' initial surf temperat in the north-south rn_tsu_ini = ', rn_tsu_ini_n,rn_tsu_ini_s
	489	WRITE(numout,*) ' initial ice temperat in the north-south rn_tmi_ini = ', rn_tmi_ini_n,rn_tmi_ini_s
	490	WRITE(numout,*) ' initial snw temperat in the north-south rn_tms_ini = ', rn_tms_ini_n,rn_tms_ini_s
	491	WRITE(numout,*) ' initial pnd fraction in the north-south rn_apd_ini = ', rn_apd_ini_n,rn_apd_ini_s
	492	WRITE(numout,*) ' initial pnd depth in the north-south rn_hpd_ini = ', rn_hpd_ini_n,rn_hpd_ini_s
	493	ENDIF
[8586]	494	ENDIF
[8813]	495	!
[8586]	496	IF( ln_iceini_file ) THEN ! Ice initialization using input file
	497	!
	498	! set si structure
	499	ALLOCATE( si(jpfldi), STAT=ierror )
	500	IF( ierror > 0 ) THEN
[11413]	501	CALL ctl_stop( 'ice_istate_ini in iceistate: unable to allocate si structure' ) ; RETURN
[8586]	502	ENDIF
[8813]	503	!
[8586]	504	DO ifpr = 1, jpfldi
	505	ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
[11413]	506	IF( slf_i(ifpr)%ln_tint ) ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
[8586]	507	END DO
[8813]	508	!
[8586]	509	! fill si with slf_i and control print
[11413]	510	CALL fld_fill( si, slf_i, cn_dir, 'ice_istate_ini', 'initialization of sea ice fields', 'numnam_ice' )
[8813]	511	!
[8586]	512	ENDIF
[8813]	513	!
[8586]	514	END SUBROUTINE ice_istate_init
	515
	516	#else
	517	!!----------------------------------------------------------------------
[9570]	518	!! Default option : Empty module NO SI3 sea-ice model
[8586]	519	!!----------------------------------------------------------------------
	520	#endif
	521
	522	!!======================================================================
	523	END MODULE iceistate

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

Download in other formats: