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.

sbcice_lim.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90 @ 8291

Last change on this file since 8291 was 8291, checked in by clem, 7 years ago
correct last commit changes. This is the start version for the modifications of the ice model. All sette tests passed except restartability for agrif and sas (though sas_biper is restartable). Configurations CREG, SPITZ, FER & SASBIPER are running.
Property svn:keywords set to `Id`
File size: 34.9 KB

Rev	Line
[888]	1	MODULE sbcice_lim
	2	!!======================================================================
	3	!! * MODULE sbcice_lim *
[921]	4	!! Surface module : update the ocean surface boundary condition over ice
	5	!! & covered area using LIM sea-ice model
[2715]	6	!! Sea-Ice model : LIM-3 Sea ice model time-stepping
[2528]	7	!!=====================================================================
	8	!! History : 2.0 ! 2006-12 (M. Vancoppenolle) Original code
	9	!! 3.0 ! 2008-02 (C. Talandier) Surface module from icestp.F90
	10	!! - ! 2008-04 (G. Madec) sltyle and lim_ctl routine
	11	!! 3.3 ! 2010-11 (G. Madec) ice-ocean stress always computed at each ocean time-step
[3625]	12	!! 3.4 ! 2011-01 (A Porter) dynamical allocation
[4161]	13	!! - ! 2012-10 (C. Rousset) add lim_diahsb
[4990]	14	!! 3.6 ! 2014-07 (M. Vancoppenolle, G. Madec, O. Marti) revise coupled interface
[7646]	15	!! 4.0 ! 2016-06 (L. Brodeau) new unified bulk routine (based on AeroBulk)
[888]	16	!!----------------------------------------------------------------------
	17	#if defined key_lim3
	18	!!----------------------------------------------------------------------
	19	!! 'key_lim3' : LIM 3.0 sea-ice model
	20	!!----------------------------------------------------------------------
[921]	21	!! sbc_ice_lim : sea-ice model time-stepping and update ocean sbc over ice-covered area
[888]	22	!!----------------------------------------------------------------------
	23	USE oce ! ocean dynamics and tracers
	24	USE dom_oce ! ocean space and time domain
[2528]	25	USE ice ! LIM-3: ice variables
[5123]	26	USE thd_ice ! LIM-3: thermodynamical variables
[6140]	27	!
[888]	28	USE sbc_oce ! Surface boundary condition: ocean fields
	29	USE sbc_ice ! Surface boundary condition: ice fields
[7646]	30	USE usrdef_sbc ! user defined: surface boundary condition
	31	USE sbcblk ! Surface boundary condition: bulk
[4161]	32	USE sbccpl ! Surface boundary condition: coupled interface
[2528]	33	USE albedo ! ocean & ice albedo
[6140]	34	!
[888]	35	USE phycst ! Define parameters for the routines
	36	USE eosbn2 ! equation of state
	37	USE limdyn ! Ice dynamics
	38	USE limtrp ! Ice transport
[5429]	39	USE limhdf ! Ice horizontal diffusion
[888]	40	USE limthd ! Ice thermodynamics
	41	USE limitd_me ! Mechanics on ice thickness distribution
	42	USE limsbc ! sea surface boundary condition
[4161]	43	USE limdiahsb ! Ice budget diagnostics
[888]	44	USE limwri ! Ice outputs
	45	USE limrst ! Ice restarts
[5123]	46	USE limupdate1 ! update of global variables
	47	USE limupdate2 ! update of global variables
[888]	48	USE limvar ! Ice variables switch
[7646]	49	USE limctl !
[8233]	50	! MV MP 2016
	51	USE limmp
	52	! END MV MP 2016
	53
[5123]	54	USE limistate ! LIM initial state
	55	USE limthd_sal ! LIM ice thermodynamics: salinity
[6140]	56	!
[2528]	57	USE c1d ! 1D vertical configuration
[6140]	58	USE in_out_manager ! I/O manager
	59	USE iom ! I/O manager library
	60	USE prtctl ! Print control
[7646]	61	USE lib_fortran !
[2528]	62	USE lbclnk ! lateral boundary condition - MPP link
[2715]	63	USE lib_mpp ! MPP library
[3294]	64	USE wrk_nemo ! work arrays
[4161]	65	USE timing ! Timing
[888]	66
[7646]	67	USE bdy_oce , ONLY: ln_bdy
	68	USE bdyice_lim ! unstructured open boundary data (bdy_ice_lim routine)
	69	# if defined key_agrif
	70	USE agrif_ice
	71	USE agrif_lim3_update
	72	USE agrif_lim3_interp
	73	# endif
[4161]	74
[888]	75	IMPLICIT NONE
	76	PRIVATE
	77
	78	PUBLIC sbc_ice_lim ! routine called by sbcmod.F90
[5123]	79	PUBLIC sbc_lim_init ! routine called by sbcmod.F90
[7646]	80
[888]	81	!! * Substitutions
	82	# include "vectopt_loop_substitute.h90"
	83	!!----------------------------------------------------------------------
[2715]	84	!! NEMO/OPA 4.0 , UCL NEMO Consortium (2011)
[1146]	85	!! $Id$
[2528]	86	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[888]	87	!!----------------------------------------------------------------------
	88	CONTAINS
	89
[7646]	90	SUBROUTINE sbc_ice_lim( kt, ksbc )
[888]	91	!!---------------------------------------------------------------------
	92	!! * ROUTINE sbc_ice_lim *
	93	!!
[7646]	94	!! ** Purpose : update the ocean surface boundary condition via the
	95	!! Louvain la Neuve Sea Ice Model time stepping
	96	!!
[888]	97	!! ** Method : ice model time stepping
[7646]	98	!! - call the ice dynamics routine
	99	!! - call the ice advection/diffusion routine
	100	!! - call the ice thermodynamics routine
	101	!! - call the routine that computes mass and
[888]	102	!! heat fluxes at the ice/ocean interface
[7646]	103	!! - save the outputs
[888]	104	!! - save the outputs for restart when necessary
	105	!!
	106	!! ** Action : - time evolution of the LIM sea-ice model
	107	!! - update all sbc variables below sea-ice:
[7646]	108	!! utau, vtau, taum, wndm, qns , qsr, emp , sfx
[888]	109	!!---------------------------------------------------------------------
	110	INTEGER, INTENT(in) :: kt ! ocean time step
[7646]	111	INTEGER, INTENT(in) :: ksbc ! type of sbc flux ( 1 = user defined formulation,
	112	! 3 = bulk formulation,
	113	! 4 = Pure Coupled formulation)
[888]	114	!!
[7753]	115	INTEGER :: jl ! dummy loop index
[4990]	116	REAL(wp), POINTER, DIMENSION(:,:,:) :: zalb_os, zalb_cs ! ice albedo under overcast/clear sky
[5407]	117	REAL(wp), POINTER, DIMENSION(:,: ) :: zutau_ice, zvtau_ice
[888]	118	!!----------------------------------------------------------------------
	119
[7646]	120	IF( nn_timing == 1 ) CALL timing_start('sbc_ice_lim')
[4161]	121
[5410]	122	!-----------------------!
	123	! --- Ice time step --- !
	124	!-----------------------!
	125	IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
[5407]	126
[7646]	127	# if defined key_agrif
	128	IF( .NOT. Agrif_Root() ) lim_nbstep = MOD( lim_nbstep, Agrif_irhot() * Agrif_Parent(nn_fsbc) / nn_fsbc ) + 1
	129	# endif
	130
[5410]	131	! mean surface ocean current at ice velocity point (C-grid dynamics : U- & V-points as the ocean)
[7753]	132	u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)
	133	v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)
[7646]	134
[5123]	135	! masked sea surface freezing temperature [Kelvin] (set to rt0 over land)
[5541]	136	CALL eos_fzp( sss_m(:,:) , t_bo(:,:) )
[7753]	137	t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
[7646]	138
[5123]	139	! Mask sea ice surface temperature (set to rt0 over land)
[1109]	140	DO jl = 1, jpl
[7753]	141	t_su(:,:,jl) = t_su(:,:,jl) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
[7646]	142	END DO
[888]	143	!
[5410]	144	!------------------------------------------------!
[7646]	145	! --- Dynamical coupling with the atmosphere --- !
	146	!------------------------------------------------!
[5410]	147	! It provides the following fields:
	148	! utau_ice, vtau_ice : surface ice stress (U- & V-points) [N/m2]
	149	!-----------------------------------------------------------------
[7646]	150	CALL sbc_lim_bef ! Store previous ice values
	151	SELECT CASE( ksbc )
	152	CASE( jp_usr ) ; CALL usrdef_sbc_ice_tau( kt ) ! user defined formulation
	153	CASE( jp_blk ) ; CALL blk_ice_tau ! Bulk formulation
	154	CASE( jp_purecpl ) ; CALL sbc_cpl_ice_tau( utau_ice , vtau_ice ) ! Coupled formulation
[888]	155	END SELECT
[7646]	156
	157	IF( ln_mixcpl) THEN ! Case of a mixed Bulk/Coupled formulation
[5407]	158	CALL wrk_alloc( jpi,jpj , zutau_ice, zvtau_ice)
[7646]	159	CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice )
[7753]	160	utau_ice(:,:) = utau_ice(:,:) * xcplmask(:,:,0) + zutau_ice(:,:) * ( 1. - xcplmask(:,:,0) )
	161	vtau_ice(:,:) = vtau_ice(:,:) * xcplmask(:,:,0) + zvtau_ice(:,:) * ( 1. - xcplmask(:,:,0) )
[5407]	162	CALL wrk_dealloc( jpi,jpj , zutau_ice, zvtau_ice)
	163	ENDIF
	164
[5410]	165	!-------------------------------------------------------!
	166	! --- ice dynamics and transport (except in 1D case) ---!
	167	!-------------------------------------------------------!
	168	numit = numit + nn_fsbc ! Ice model time step
[921]	169	!
[7646]	170	CALL sbc_lim_diag0 ! set diag of mass, heat and salt fluxes to 0
	171	CALL lim_rst_opn( kt ) ! Open Ice restart file
	172	!
	173	! --- zap this if no ice dynamics --- !
	174	IF( .NOT. lk_c1d .AND. ln_limdyn ) THEN
[5410]	175	!
[7646]	176	IF( nn_limdyn /= 0 ) THEN ! -- Ice dynamics
	177	CALL lim_dyn( kt ) ! rheology
	178	ELSE
[7753]	179	u_ice(:,:) = rn_uice * umask(:,:,1) ! or prescribed velocity
	180	v_ice(:,:) = rn_vice * vmask(:,:,1)
[8291]	181	!!CALL RANDOM_NUMBER(u_ice(:,:))
	182	!!CALL RANDOM_NUMBER(v_ice(:,:))
[7646]	183	ENDIF
	184	CALL lim_trp( kt ) ! -- Ice transport (Advection/diffusion)
	185	IF( nn_limdyn == 2 .AND. nn_monocat /= 2 ) & ! -- Mechanical redistribution (ridging/rafting)
	186	& CALL lim_itd_me
	187	IF( nn_limdyn == 2 ) CALL lim_update1( kt ) ! -- Corrections
[5410]	188	!
[7646]	189	ENDIF
	190
	191	! ---
	192	#if defined key_agrif
	193	IF( .NOT. Agrif_Root() ) CALL agrif_interp_lim3('T')
[4333]	194	#endif
[7761]	195	IF( ln_limthd .AND. ln_bdy ) CALL bdy_ice_lim( kt ) ! -- bdy ice thermo
[5123]	196	! previous lead fraction and ice volume for flux calculations
[7646]	197	CALL sbc_lim_bef
	198	CALL lim_var_glo2eqv ! ht_i and ht_s for ice albedo calculation
	199	CALL lim_var_agg(1) ! at_i for coupling (via pfrld)
	200	!
[7753]	201	pfrld(:,:) = 1._wp - at_i(:,:)
	202	phicif(:,:) = vt_i(:,:)
[7646]	203
[5410]	204	!------------------------------------------------------!
[7646]	205	! --- Thermodynamical coupling with the atmosphere --- !
	206	!------------------------------------------------------!
[5410]	207	! It provides the following fields:
	208	! qsr_ice , qns_ice : solar & non solar heat flux over ice (T-point) [W/m2]
	209	! qla_ice : latent heat flux over ice (T-point) [W/m2]
	210	! dqns_ice, dqla_ice : non solar & latent heat sensistivity (T-point) [W/m2]
	211	! tprecip , sprecip : total & solid precipitation (T-point) [Kg/m2/s]
	212	! fr1_i0 , fr2_i0 : 1sr & 2nd fraction of qsr penetration in ice [%]
	213	!----------------------------------------------------------------------------------------
[6416]	214	CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs )
[7646]	215
[8233]	216	CALL albedo_ice( t_su, ht_i, ht_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos MV MP 2016
	217
[7646]	218	SELECT CASE( ksbc )
	219	CASE( jp_usr ) ; CALL usrdef_sbc_ice_flx( kt ) ! user defined formulation
	220	CASE( jp_blk ) ! bulk formulation
	221	! albedo depends on cloud fraction because of non-linear spectral effects
[7753]	222	alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
[7646]	223	CALL blk_ice_flx( t_su, alb_ice )
	224	IF( ln_mixcpl ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su )
	225	IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
	226	CASE ( jp_purecpl )
	227	! albedo depends on cloud fraction because of non-linear spectral effects
[7753]	228	alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
[7646]	229	CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su )
	230	IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
	231	END SELECT
[5410]	232
[6416]	233	CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs )
[5407]	234
[5410]	235	!----------------------------!
	236	! --- ice thermodynamics --- !
	237	!----------------------------!
[7646]	238	! --- zap this if no ice thermo --- !
	239	IF( ln_limthd ) CALL lim_thd( kt ) ! -- Ice thermodynamics
[8233]	240
	241	! MV MP 2016
	242	IF ( ln_pnd ) CALL lim_mp( kt ) ! -- Melt ponds
	243	! END MV MP 2016
	244
[7646]	245	IF( ln_limthd ) CALL lim_update2( kt ) ! -- Corrections
	246	! ---
	247	# if defined key_agrif
	248	IF( .NOT. Agrif_Root() ) CALL agrif_update_lim3( kt )
	249	# endif
	250	CALL lim_var_glo2eqv ! necessary calls (at least for coupling)
	251	CALL lim_var_agg( 2 ) ! necessary calls (at least for coupling)
	252	!
[7761]	253	!! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work
	254	!! moreover it should only be called at the update frequency (as in agrif_lim3_update.F90)
	255	!!# if defined key_agrif
	256	!! IF( .NOT. Agrif_Root() ) CALL Agrif_ChildGrid_To_ParentGrid()
	257	!!# endif
[7646]	258	CALL lim_sbc_flx( kt ) ! -- Update surface ocean mass, heat and salt fluxes
[7761]	259	!!# if defined key_agrif
	260	!! IF( .NOT. Agrif_Root() ) CALL Agrif_ParentGrid_To_ChildGrid()
	261	!!# endif
[7646]	262	IF( ln_limdiahsb ) CALL lim_diahsb( kt ) ! -- Diagnostics and outputs
[921]	263	!
[7646]	264	CALL lim_wri( 1 ) ! -- Ice outputs
[888]	265	!
[4205]	266	IF( kt == nit000 .AND. ln_rstart ) &
[7646]	267	& CALL iom_close( numrir ) ! close input ice restart file
[4205]	268	!
[7646]	269	IF( lrst_ice ) CALL lim_rst_write( kt ) ! -- Ice restart file
[921]	270	!
[7646]	271	IF( ln_limctl ) CALL lim_ctl( kt ) ! alerts in case of model crash
[921]	272	!
[5123]	273	ENDIF ! End sea-ice time step only
[888]	274
[5410]	275	!-------------------------!
	276	! --- Ocean time step --- !
	277	!-------------------------!
	278	! Update surface ocean stresses (only in ice-dynamic case) otherwise the atm.-ocean stresses are used everywhere
[7646]	279	! using before instantaneous surf. currents
	280	IF( ln_limdyn ) CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) )
[2528]	281	!!gm remark, the ocean-ice stress is not saved in ice diag call above ..... find a solution!!!
[2715]	282	!
[5410]	283	IF( nn_timing == 1 ) CALL timing_stop('sbc_ice_lim')
[4161]	284	!
[921]	285	END SUBROUTINE sbc_ice_lim
[5123]	286
[7646]	287
[5123]	288	SUBROUTINE sbc_lim_init
	289	!!----------------------------------------------------------------------
	290	!! * ROUTINE sbc_lim_init *
	291	!!
	292	!! ** purpose : Allocate all the dynamic arrays of the LIM-3 modules
	293	!!----------------------------------------------------------------------
[7753]	294	INTEGER :: ji, jj, ierr
[5123]	295	!!----------------------------------------------------------------------
	296	IF(lwp) WRITE(numout,*)
[7646]	297	IF(lwp) WRITE(numout,*) 'sbc_lim_init : update ocean surface boudary condition'
[5123]	298	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ via Louvain la Neuve Ice Model (LIM-3) time stepping'
	299	!
[7646]	300	! ! Open the reference and configuration namelist files and namelist output file
	301	CALL ctl_opn( numnam_ice_ref, 'namelist_ice_ref', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
[5123]	302	CALL ctl_opn( numnam_ice_cfg, 'namelist_ice_cfg', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
	303	IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 )
[6140]	304	!
[7646]	305	CALL lim_run_init ! set some ice run parameters
[5123]	306	!
	307	! ! Allocate the ice arrays
	308	ierr = ice_alloc () ! ice variables
	309	ierr = ierr + sbc_ice_alloc () ! surface forcing
	310	ierr = ierr + thd_ice_alloc () ! thermodynamics
[7646]	311	IF( ln_limdyn ) ierr = ierr + lim_itd_me_alloc () ! ice thickness distribution - mechanics
[5123]	312	!
	313	IF( lk_mpp ) CALL mpp_sum( ierr )
	314	IF( ierr /= 0 ) CALL ctl_stop('STOP', 'sbc_lim_init : unable to allocate ice arrays')
	315	!
[7646]	316	CALL lim_dyn_init ! set ice dynamics parameters
[5123]	317	!
	318	CALL lim_itd_init ! ice thickness distribution initialization
	319	!
[5429]	320	CALL lim_hdf_init ! set ice horizontal diffusion computation parameters
	321	!
[5123]	322	CALL lim_thd_init ! set ice thermodynics parameters
	323	!
	324	CALL lim_thd_sal_init ! set ice salinity parameters
[8233]	325
	326	! MV MP 2016
	327	CALL lim_mp_init ! set melt ponds parameters
	328	! END MV MP 2016
	329
[7646]	330	IF( ln_limdyn ) CALL lim_itd_me_init ! ice thickness distribution initialization for mecanical deformation
[5123]	331	! ! Initial sea-ice state
	332	IF( .NOT. ln_rstart ) THEN ! start from rest: sea-ice deduced from sst
	333	numit = 0
	334	numit = nit000 - 1
	335	CALL lim_istate
	336	ELSE ! start from a restart file
	337	CALL lim_rst_read
	338	numit = nit000 - 1
	339	ENDIF
[7646]	340	CALL lim_var_agg(2)
[5123]	341	CALL lim_var_glo2eqv
	342	!
[7646]	343	CALL lim_sbc_init ! ice surface boundary condition
[5123]	344	!
[7646]	345	IF( ln_limdiahsb) CALL lim_diahsb_init ! initialization for diags
	346	!
[7753]	347	fr_i(:,:) = at_i(:,:) ! initialisation of sea-ice fraction
	348	tn_ice(:,:,:) = t_su(:,:,:) ! initialisation of surface temp for coupled simu
[5123]	349	!
[6403]	350	DO jj = 1, jpj
	351	DO ji = 1, jpi
	352	IF( gphit(ji,jj) > 0._wp ) THEN ; rn_amax_2d(ji,jj) = rn_amax_n ! NH
	353	ELSE ; rn_amax_2d(ji,jj) = rn_amax_s ! SH
	354	ENDIF
[7646]	355	END DO
	356	END DO
[6403]	357	!
[7646]	358	nstart = numit + nn_fsbc
	359	nitrun = nitend - nit000 + 1
	360	nlast = numit + nitrun
[5123]	361	!
	362	IF( nstock == 0 ) nstock = nlast + 1
	363	!
	364	END SUBROUTINE sbc_lim_init
	365
	366
[7646]	367	SUBROUTINE lim_run_init
[5123]	368	!!-------------------------------------------------------------------
[7646]	369	!! * ROUTINE lim_run_init *
	370	!!
[5123]	371	!! ** Purpose : Definition some run parameter for ice model
	372	!!
[7646]	373	!! ** Method : Read the namicerun namelist and check the parameter
[5123]	374	!! values called at the first timestep (nit000)
	375	!!
	376	!! ** input : Namelist namicerun
	377	!!-------------------------------------------------------------------
	378	INTEGER :: ios ! Local integer output status for namelist read
[7646]	379	NAMELIST/namicerun/ jpl, nlay_i, nlay_s, rn_amax_n, rn_amax_s, cn_icerst_in, cn_icerst_indir, &
	380	& cn_icerst_out, cn_icerst_outdir, ln_limthd, ln_limdyn, nn_limdyn, rn_uice, rn_vice
	381	NAMELIST/namicediag/ ln_limdiachk, ln_limdiahsb, ln_limctl, iiceprt, jiceprt
[5123]	382	!!-------------------------------------------------------------------
[7646]	383	!
[5123]	384	REWIND( numnam_ice_ref ) ! Namelist namicerun in reference namelist : Parameters for ice
	385	READ ( numnam_ice_ref, namicerun, IOSTAT = ios, ERR = 901)
[7646]	386	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp )
	387
[5123]	388	REWIND( numnam_ice_cfg ) ! Namelist namicerun in configuration namelist : Parameters for ice
	389	READ ( numnam_ice_cfg, namicerun, IOSTAT = ios, ERR = 902 )
[7646]	390	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp )
	391	IF(lwm) WRITE ( numoni, namicerun )
[5123]	392	!
[7646]	393	REWIND( numnam_ice_ref ) ! Namelist namicediag in reference namelist : Parameters for ice
	394	READ ( numnam_ice_ref, namicediag, IOSTAT = ios, ERR = 903)
	395	903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in reference namelist', lwp )
	396
	397	REWIND( numnam_ice_cfg ) ! Namelist namicediag in configuration namelist : Parameters for ice
	398	READ ( numnam_ice_cfg, namicediag, IOSTAT = ios, ERR = 904 )
	399	904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in configuration namelist', lwp )
	400	IF(lwm) WRITE ( numoni, namicediag )
	401	!
[5123]	402	IF(lwp) THEN ! control print
	403	WRITE(numout,*)
[7646]	404	WRITE(numout,*) 'lim_run_init : ice share parameters for dynamics/advection/thermo of sea-ice'
[5123]	405	WRITE(numout,*) ' ~~~~~~'
	406	WRITE(numout,*) ' number of ice categories = ', jpl
	407	WRITE(numout,*) ' number of ice layers = ', nlay_i
	408	WRITE(numout,*) ' number of snow layers = ', nlay_s
[6403]	409	WRITE(numout,*) ' maximum ice concentration for NH = ', rn_amax_n
	410	WRITE(numout,*) ' maximum ice concentration for SH = ', rn_amax_s
[7646]	411	WRITE(numout,*) ' Ice thermodynamics (T) or not (F) ln_limthd = ', ln_limthd
	412	WRITE(numout,*) ' Ice dynamics (T) or not (F) ln_limdyn = ', ln_limdyn
	413	WRITE(numout,*) ' (ln_limdyn=T) Ice dynamics switch nn_limdyn = ', nn_limdyn
	414	WRITE(numout,*) ' 2: total'
	415	WRITE(numout,*) ' 1: advection only (no diffusion, no ridging/rafting)'
	416	WRITE(numout,*) ' 0: advection only (as 1 + prescribed velocity, bypass rheology)'
	417	WRITE(numout,*) ' (ln_limdyn=T) prescribed u-vel (case nn_limdyn=0) = ', rn_uice
	418	WRITE(numout,*) ' (ln_limdyn=T) prescribed v-vel (case nn_limdyn=0) = ', rn_vice
	419	WRITE(numout,*)
	420	WRITE(numout,*) '...and ice diagnostics'
	421	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~'
	422	WRITE(numout,*) ' Diagnose online heat/mass/salt budget ln_limdiachk = ', ln_limdiachk
	423	WRITE(numout,*) ' Output heat/mass/salt budget ln_limdiahsb = ', ln_limdiahsb
	424	WRITE(numout,*) ' control prints in ocean.out for (i,j)=(iiceprt,jiceprt) = ', ln_limctl
	425	WRITE(numout,*) ' i-index for control prints (ln_limctl=true) = ', iiceprt
	426	WRITE(numout,*) ' j-index for control prints (ln_limctl=true) = ', jiceprt
[5123]	427	ENDIF
	428	!
	429	! sea-ice timestep and inverse
[7646]	430	rdt_ice = REAL(nn_fsbc) * rdt
[5123]	431	r1_rdtice = 1._wp / rdt_ice
	432
	433	! inverse of nlay_i and nlay_s
	434	r1_nlay_i = 1._wp / REAL( nlay_i, wp )
	435	r1_nlay_s = 1._wp / REAL( nlay_s, wp )
	436	!
[7646]	437	IF( lwp .AND. ln_bdy .AND. ln_limdiachk ) &
[7761]	438	& CALL ctl_warn('online conservation check activated but it does not work with BDY')
[5167]	439	!
[7646]	440	IF( lwp ) WRITE(numout,*) ' ice timestep rdt_ice = ', rdt_ice
	441	!
	442	END SUBROUTINE lim_run_init
[5123]	443
	444
	445	SUBROUTINE lim_itd_init
	446	!!------------------------------------------------------------------
	447	!! * ROUTINE lim_itd_init *
	448	!!
	449	!! ** Purpose : Initializes the ice thickness distribution
	450	!! ** Method : ...
	451	!! ** input : Namelist namiceitd
	452	!!-------------------------------------------------------------------
	453	INTEGER :: ios ! Local integer output status for namelist read
	454	NAMELIST/namiceitd/ nn_catbnd, rn_himean
	455	!
	456	INTEGER :: jl ! dummy loop index
	457	REAL(wp) :: zc1, zc2, zc3, zx1 ! local scalars
	458	REAL(wp) :: zhmax, znum, zden, zalpha !
	459	!!------------------------------------------------------------------
	460	!
	461	REWIND( numnam_ice_ref ) ! Namelist namiceitd in reference namelist : Parameters for ice
[7646]	462	READ ( numnam_ice_ref, namiceitd, IOSTAT = ios, ERR = 905)
	463	905 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in reference namelist', lwp )
	464
[5123]	465	REWIND( numnam_ice_cfg ) ! Namelist namiceitd in configuration namelist : Parameters for ice
[7646]	466	READ ( numnam_ice_cfg, namiceitd, IOSTAT = ios, ERR = 906 )
	467	906 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in configuration namelist', lwp )
	468	IF(lwm) WRITE ( numoni, namiceitd )
[5123]	469	!
	470	IF(lwp) THEN ! control print
	471	WRITE(numout,*)
[7646]	472	WRITE(numout,*) 'lim_itd_init : Initialization of ice cat distribution '
	473	WRITE(numout,*) '~~~~~~~~~~~~'
	474	WRITE(numout,*) ' shape of ice categories distribution nn_catbnd = ', nn_catbnd
	475	WRITE(numout,*) ' mean ice thickness in the domain (only active if nn_catbnd=2) rn_himean = ', rn_himean
[5123]	476	ENDIF
[6140]	477	!
[5123]	478	!----------------------------------
[7646]	479	!- Thickness categories boundaries
[5123]	480	!----------------------------------
[6140]	481	!
[5123]	482	hi_max(:) = 0._wp
[6140]	483	!
	484	SELECT CASE ( nn_catbnd ) ! type of ice categories distribution
	485	!
	486	CASE (1) !== tanh function (CICE) ==!
[5123]	487	zc1 = 3._wp / REAL( jpl, wp )
	488	zc2 = 10._wp * zc1
	489	zc3 = 3._wp
	490	DO jl = 1, jpl
	491	zx1 = REAL( jl-1, wp ) / REAL( jpl, wp )
	492	hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1._wp + TANH( zc3 * (zx1 - 1._wp ) ) )
	493	END DO
[6140]	494	!
	495	CASE (2) !== h^(-alpha) function ==!
	496	zalpha = 0.05_wp
	497	zhmax = 3._wp * rn_himean
[7646]	498	DO jl = 1, jpl
[5123]	499	znum = jpl * ( zhmax+1 )**zalpha
[6140]	500	zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
[5123]	501	hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
	502	END DO
[6140]	503	!
[5123]	504	END SELECT
[6140]	505	!
	506	DO jl = 1, jpl ! mean thickness by category
[5123]	507	hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
	508	END DO
[6140]	509	!
	510	hi_max(jpl) = 99._wp ! set to a big value to ensure that all ice is thinner than hi_max(jpl)
	511	!
[7646]	512	IF(lwp) WRITE(numout,*)
	513	IF(lwp) WRITE(numout,*) ' Thickness category boundaries '
	514	IF(lwp) WRITE(numout,*) ' hi_max ', hi_max(0:jpl)
[5123]	515	!
	516	END SUBROUTINE lim_itd_init
	517
[7646]	518
	519	SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_limflx )
[4990]	520	!!---------------------------------------------------------------------
[5123]	521	!! * ROUTINE ice_lim_flx *
[7646]	522	!!
[4990]	523	!! ** Purpose : update the ice surface boundary condition by averaging and / or
[7646]	524	!! redistributing fluxes on ice categories
[4990]	525	!!
[7646]	526	!! ** Method : average then redistribute
[4990]	527	!!
[7646]	528	!! ** Action :
[4990]	529	!!---------------------------------------------------------------------
[7646]	530	INTEGER , INTENT(in ) :: k_limflx ! =-1 do nothing; =0 average ;
	531	! ! = 1 average and redistribute ; =2 redistribute
	532	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptn_ice ! ice surface temperature
[4990]	533	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: palb_ice ! ice albedo
	534	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pqns_ice ! non solar flux
	535	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pqsr_ice ! net solar flux
	536	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdqn_ice ! non solar flux sensitivity
[5407]	537	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pevap_ice ! sublimation
	538	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdevap_ice ! sublimation sensitivity
[4990]	539	!
[7753]	540	INTEGER :: jl ! dummy loop index
[4990]	541	!
	542	REAL(wp), POINTER, DIMENSION(:,:) :: zalb_m ! Mean albedo over all categories
	543	REAL(wp), POINTER, DIMENSION(:,:) :: ztem_m ! Mean temperature over all categories
	544	!
	545	REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_m ! Mean solar heat flux over all categories
	546	REAL(wp), POINTER, DIMENSION(:,:) :: z_qns_m ! Mean non solar heat flux over all categories
[5407]	547	REAL(wp), POINTER, DIMENSION(:,:) :: z_evap_m ! Mean sublimation over all categories
[4990]	548	REAL(wp), POINTER, DIMENSION(:,:) :: z_dqn_m ! Mean d(qns)/dT over all categories
[5407]	549	REAL(wp), POINTER, DIMENSION(:,:) :: z_devap_m ! Mean d(evap)/dT over all categories
[4990]	550	!!----------------------------------------------------------------------
[6140]	551	!
[4990]	552	IF( nn_timing == 1 ) CALL timing_start('ice_lim_flx')
	553	!
	554	SELECT CASE( k_limflx ) !== averaged on all ice categories ==!
	555	CASE( 0 , 1 )
[5407]	556	CALL wrk_alloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m)
[4990]	557	!
[5410]	558	z_qns_m (:,:) = fice_ice_ave ( pqns_ice (:,:,:) )
	559	z_qsr_m (:,:) = fice_ice_ave ( pqsr_ice (:,:,:) )
	560	z_dqn_m (:,:) = fice_ice_ave ( pdqn_ice (:,:,:) )
	561	z_evap_m (:,:) = fice_ice_ave ( pevap_ice (:,:,:) )
[5407]	562	z_devap_m(:,:) = fice_ice_ave ( pdevap_ice (:,:,:) )
[4990]	563	DO jl = 1, jpl
[7753]	564	pdqn_ice (:,:,jl) = z_dqn_m(:,:)
	565	pdevap_ice(:,:,jl) = z_devap_m(:,:)
[4990]	566	END DO
	567	!
	568	DO jl = 1, jpl
[7753]	569	pqns_ice (:,:,jl) = z_qns_m(:,:)
	570	pqsr_ice (:,:,jl) = z_qsr_m(:,:)
	571	pevap_ice(:,:,jl) = z_evap_m(:,:)
[4990]	572	END DO
	573	!
[5407]	574	CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m)
[4990]	575	END SELECT
[6140]	576	!
[4990]	577	SELECT CASE( k_limflx ) !== redistribution on all ice categories ==!
	578	CASE( 1 , 2 )
	579	CALL wrk_alloc( jpi,jpj, zalb_m, ztem_m )
	580	!
[7646]	581	zalb_m(:,:) = fice_ice_ave ( palb_ice (:,:,:) )
	582	ztem_m(:,:) = fice_ice_ave ( ptn_ice (:,:,:) )
[4990]	583	DO jl = 1, jpl
[7753]	584	pqns_ice (:,:,jl) = pqns_ice (:,:,jl) + pdqn_ice (:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
	585	pevap_ice(:,:,jl) = pevap_ice(:,:,jl) + pdevap_ice(:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
	586	pqsr_ice (:,:,jl) = pqsr_ice (:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) )
[4990]	587	END DO
	588	!
	589	CALL wrk_dealloc( jpi,jpj, zalb_m, ztem_m )
	590	END SELECT
	591	!
	592	IF( nn_timing == 1 ) CALL timing_stop('ice_lim_flx')
	593	!
	594	END SUBROUTINE ice_lim_flx
[888]	595
[6140]	596
[5167]	597	SUBROUTINE sbc_lim_bef
[5123]	598	!!----------------------------------------------------------------------
[5167]	599	!! * ROUTINE sbc_lim_bef *
[5123]	600	!!
[7646]	601	!! ** purpose : store ice variables at "before" time step
[5123]	602	!!----------------------------------------------------------------------
[7753]	603	a_i_b (:,:,:) = a_i (:,:,:) ! ice area
	604	e_i_b (:,:,:,:) = e_i (:,:,:,:) ! ice thermal energy
	605	v_i_b (:,:,:) = v_i (:,:,:) ! ice volume
	606	v_s_b (:,:,:) = v_s (:,:,:) ! snow volume
	607	e_s_b (:,:,:,:) = e_s (:,:,:,:) ! snow thermal energy
	608	smv_i_b(:,:,:) = smv_i(:,:,:) ! salt content
	609	oa_i_b (:,:,:) = oa_i (:,:,:) ! areal age content
	610	u_ice_b(:,:) = u_ice(:,:)
	611	v_ice_b(:,:) = v_ice(:,:)
	612	!
	613	at_i_b (:,:) = SUM( a_i_b(:,:,:), dim=3 )
[7646]	614
[5167]	615	END SUBROUTINE sbc_lim_bef
[921]	616
[6140]	617
[5123]	618	SUBROUTINE sbc_lim_diag0
	619	!!----------------------------------------------------------------------
	620	!! * ROUTINE sbc_lim_diag0 *
[888]	621	!!
[5123]	622	!! ** purpose : set ice-ocean and ice-atm. fluxes to zeros at the beggining
	623	!! of the time step
	624	!!----------------------------------------------------------------------
[7753]	625	sfx (:,:) = 0._wp ;
	626	sfx_bri(:,:) = 0._wp ; sfx_lam(:,:) = 0._wp
	627	sfx_sni(:,:) = 0._wp ; sfx_opw(:,:) = 0._wp
	628	sfx_bog(:,:) = 0._wp ; sfx_dyn(:,:) = 0._wp
	629	sfx_bom(:,:) = 0._wp ; sfx_sum(:,:) = 0._wp
	630	sfx_res(:,:) = 0._wp ; sfx_sub(:,:) = 0._wp
	631	!
	632	wfx_snw(:,:) = 0._wp ; wfx_ice(:,:) = 0._wp
	633	wfx_sni(:,:) = 0._wp ; wfx_opw(:,:) = 0._wp
	634	wfx_bog(:,:) = 0._wp ; wfx_dyn(:,:) = 0._wp
	635	wfx_bom(:,:) = 0._wp ; wfx_sum(:,:) = 0._wp
	636	wfx_res(:,:) = 0._wp ; wfx_sub(:,:) = 0._wp
	637	wfx_spr(:,:) = 0._wp ; wfx_lam(:,:) = 0._wp
[8239]	638	wfx_snw_dyn(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp
	639	wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp
[8233]	640
	641	! MV MP 2016
	642	wfx_pnd(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp
	643	! END MV MP 2016
[7753]	644
	645	hfx_thd(:,:) = 0._wp ;
	646	hfx_snw(:,:) = 0._wp ; hfx_opw(:,:) = 0._wp
	647	hfx_bog(:,:) = 0._wp ; hfx_dyn(:,:) = 0._wp
	648	hfx_bom(:,:) = 0._wp ; hfx_sum(:,:) = 0._wp
	649	hfx_res(:,:) = 0._wp ; hfx_sub(:,:) = 0._wp
	650	hfx_spr(:,:) = 0._wp ; hfx_dif(:,:) = 0._wp
	651	hfx_err(:,:) = 0._wp ; hfx_err_rem(:,:) = 0._wp
	652	hfx_err_dif(:,:) = 0._wp
	653	wfx_err_sub(:,:) = 0._wp
	654	!
	655	afx_tot(:,:) = 0._wp ;
	656	afx_dyn(:,:) = 0._wp ; afx_thd(:,:) = 0._wp
	657	!
	658	diag_heat(:,:) = 0._wp ; diag_smvi(:,:) = 0._wp
	659	diag_vice(:,:) = 0._wp ; diag_vsnw(:,:) = 0._wp
[7646]	660
[8239]	661	! SIMIP diagnostics
	662	diag_dms_dyn(:,:) = 0._wp ; diag_dmi_dyn(:,:) = 0._wp
	663	diag_fc_bo(:,:) = 0._wp ; diag_fc_su(:,:) = 0._wp
	664
[7753]	665	tau_icebfr(:,:) = 0._wp; ! landfast ice param only (clem: important to keep the init here)
[7646]	666
[5123]	667	END SUBROUTINE sbc_lim_diag0
[5407]	668
[7646]	669
[4990]	670	FUNCTION fice_cell_ave ( ptab )
	671	!!--------------------------------------------------------------------------
	672	!! * Compute average over categories, for grid cell (ice covered and free ocean)
	673	!!--------------------------------------------------------------------------
	674	REAL (wp), DIMENSION (jpi,jpj) :: fice_cell_ave
	675	REAL (wp), DIMENSION (jpi,jpj,jpl), INTENT (in) :: ptab
	676	INTEGER :: jl ! Dummy loop index
[7646]	677
	678	fice_cell_ave (:,:) = 0._wp
[4990]	679	DO jl = 1, jpl
[5123]	680	fice_cell_ave (:,:) = fice_cell_ave (:,:) + a_i (:,:,jl) * ptab (:,:,jl)
[4990]	681	END DO
[7646]	682
[4990]	683	END FUNCTION fice_cell_ave
[7646]	684
	685
[4990]	686	FUNCTION fice_ice_ave ( ptab )
	687	!!--------------------------------------------------------------------------
	688	!! * Compute average over categories, for ice covered part of grid cell
	689	!!--------------------------------------------------------------------------
	690	REAL (kind=wp), DIMENSION (jpi,jpj) :: fice_ice_ave
	691	REAL (kind=wp), DIMENSION (jpi,jpj,jpl), INTENT(in) :: ptab
[888]	692
[4990]	693	fice_ice_ave (:,:) = 0.0_wp
[5167]	694	WHERE ( at_i (:,:) > 0.0_wp ) fice_ice_ave (:,:) = fice_cell_ave ( ptab (:,:,:)) / at_i (:,:)
[4990]	695
	696	END FUNCTION fice_ice_ave
	697
[888]	698	#else
	699	!!----------------------------------------------------------------------
	700	!! Default option Dummy module NO LIM 3.0 sea-ice model
	701	!!----------------------------------------------------------------------
	702	CONTAINS
[7777]	703	!
[7646]	704	SUBROUTINE sbc_ice_lim ( kt, ksbc ) ! Dummy routine
[7777]	705	INTEGER, INTENT(in) :: kt, ksbc
[7646]	706	WRITE(,) 'sbc_ice_lim: You should not have seen this print! error?', kt
[888]	707	END SUBROUTINE sbc_ice_lim
[7777]	708	!
[5123]	709	SUBROUTINE sbc_lim_init ! Dummy routine
	710	END SUBROUTINE sbc_lim_init
[888]	711	#endif
	712
	713	!!======================================================================
	714	END MODULE sbcice_lim

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

Download in other formats: