Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

iceistate.F90 @ 15795

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

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

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

Context Navigation

source: NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/ICE/iceistate.F90 @ 15795

Download in other formats: