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/trunk/src/ICE – NEMO

Context Navigation

source: NEMO/trunk/src/ICE/iceistate.F90 @ 14137

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

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

Download in other formats: