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.

limtrp.F90 in branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 4359

Last change on this file since 4359 was 4333, checked in by clem, 11 years ago
remove remaining bugs in LIM3, so that it can run in both regional and global config
Property svn:keywords set to `Id`
File size: 31.9 KB

Rev	Line
[825]	1	MODULE limtrp
	2	!!======================================================================
	3	!! * MODULE limtrp *
	4	!! LIM transport ice model : sea-ice advection/diffusion
	5	!!======================================================================
[2715]	6	!! History : LIM-2 ! 2000-01 (M.A. Morales Maqueda, H. Goosse, and T. Fichefet) Original code
	7	!! 3.0 ! 2005-11 (M. Vancoppenolle) Multi-layer sea ice, salinity variations
	8	!! 4.0 ! 2011-02 (G. Madec) dynamical allocation
	9	!!----------------------------------------------------------------------
[825]	10	#if defined key_lim3
	11	!!----------------------------------------------------------------------
[834]	12	!! 'key_lim3' LIM3 sea-ice model
[825]	13	!!----------------------------------------------------------------------
	14	!! lim_trp : advection/diffusion process of sea ice
	15	!!----------------------------------------------------------------------
[3625]	16	USE phycst ! physical constant
	17	USE dom_oce ! ocean domain
	18	USE sbc_oce ! ocean surface boundary condition
	19	USE par_ice ! ice parameter
	20	USE dom_ice ! ice domain
	21	USE ice ! ice variables
	22	USE limadv ! ice advection
	23	USE limhdf ! ice horizontal diffusion
	24	USE in_out_manager ! I/O manager
	25	USE lbclnk ! lateral boundary conditions -- MPP exchanges
	26	USE lib_mpp ! MPP library
	27	USE wrk_nemo ! work arrays
	28	USE prtctl ! Print control
	29	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[4161]	30	USE limvar ! clem for ice thickness correction
	31	USE timing ! Timing
[825]	32
	33	IMPLICIT NONE
	34	PRIVATE
	35
[2715]	36	PUBLIC lim_trp ! called by ice_step
[825]	37
[4161]	38	REAL(wp) :: epsi10 = 1.e-10_wp
[2715]	39	REAL(wp) :: rzero = 0._wp
	40	REAL(wp) :: rone = 1._wp
[825]	41
	42	!! * Substitution
	43	# include "vectopt_loop_substitute.h90"
	44	!!----------------------------------------------------------------------
[4161]	45	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
[1156]	46	!! $Id$
[2715]	47	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	48	!!----------------------------------------------------------------------
	49	CONTAINS
	50
[921]	51	SUBROUTINE lim_trp( kt )
[825]	52	!!-------------------------------------------------------------------
	53	!! * ROUTINE lim_trp *
	54	!!
	55	!! ** purpose : advection/diffusion process of sea ice
	56	!!
	57	!! ** method : variables included in the process are scalar,
	58	!! other values are considered as second order.
	59	!! For advection, a second order Prather scheme is used.
	60	!!
	61	!! ** action :
	62	!!---------------------------------------------------------------------
[2715]	63	INTEGER, INTENT(in) :: kt ! number of iteration
	64	!
	65	INTEGER :: ji, jj, jk, jl, layer ! dummy loop indices
	66	INTEGER :: initad ! number of sub-timestep for the advection
[2777]	67	INTEGER :: ierr ! error status
[4161]	68	REAL(wp) :: zindb , zindsn , zindic, zindh, zinda ! local scalar
[2715]	69	REAL(wp) :: zusvosn, zusvoic, zbigval ! - -
[4205]	70	REAL(wp) :: zcfl , zusnit ! - -
[2715]	71	REAL(wp) :: ze , zsal , zage ! - -
	72	!
[3294]	73	REAL(wp), POINTER, DIMENSION(:,:) :: zui_u, zvi_v, zsm, zs0at, zs0ow
	74	REAL(wp), POINTER, DIMENSION(:,:,:) :: zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi
	75	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zs0e
[4161]	76	REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)
	77	REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax, zchk_umax ! Check errors (C Rousset)
	78	! mass and salt flux (clem)
	79	REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold ! old ice volume...
	80	! correct ice thickness (clem)
	81	REAL(wp), POINTER, DIMENSION(:,:,:) :: zaiold, zhimax ! old ice concentration and thickness
	82	REAL(wp) :: zdv, zda, zvi, zvs, zsmv
[2715]	83	!!---------------------------------------------------------------------
[4161]	84	IF( nn_timing == 1 ) CALL timing_start('limtrp')
[825]	85
[3294]	86	CALL wrk_alloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow )
	87	CALL wrk_alloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi )
	88	CALL wrk_alloc( jpi, jpj, jkmax, jpl, zs0e )
[825]	89
[4161]	90	CALL wrk_alloc( jpi,jpj,jpl,zviold ) ! clem
	91	CALL wrk_alloc( jpi,jpj,jpl,zaiold, zhimax ) ! clem
	92
	93	! -------------------------------
	94	!- check conservation (C Rousset)
	95	IF( ln_limdiahsb ) THEN
	96	zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
	97	zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
	98	zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) )
	99	zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) )
	100	ENDIF
	101	!- check conservation (C Rousset)
	102	! -------------------------------
	103
[2715]	104	IF( numit == nstart .AND. lwp ) THEN
	105	WRITE(numout,*)
	106	IF( ln_limdyn ) THEN ; WRITE(numout,*) 'lim_trp : Ice transport '
	107	ELSE ; WRITE(numout,*) 'lim_trp : No ice advection as ln_limdyn = ', ln_limdyn
	108	ENDIF
	109	WRITE(numout,*) '~~~~~~~~~~~~'
	110	ENDIF
	111
[825]	112	zsm(:,:) = area(:,:)
	113
[2715]	114	! !-------------------------------------!
	115	IF( ln_limdyn ) THEN ! Advection of sea ice properties !
	116	! !-------------------------------------!
[4161]	117	! mass and salt flux init (clem)
	118	zviold(:,:,:) = v_i(:,:,:)
[825]	119
[4161]	120	!--- Thickness correction init. (clem) -------------------------------
	121	CALL lim_var_glo2eqv
	122	zaiold(:,:,:) = a_i(:,:,:)
	123	!---------------------------------------------------------------------
	124	! Record max of the surrounding ice thicknesses for correction in limupdate
	125	! in case advection creates ice too thick.
	126	!---------------------------------------------------------------------
	127	zhimax(:,:,:) = ht_i(:,:,:)
	128	DO jl = 1, jpl
	129	DO jj = 2, jpjm1
	130	DO ji = 2, jpim1
	131	zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) )
	132	!zhimax(ji,jj,jl) = ( ht_i(ji ,jj ,jl) * tmask(ji, jj ,1) + ht_i(ji-1,jj-1,jl) * tmask(ji-1,jj-1,1) + ht_i(ji+1,jj+1,jl) * tmask(ji+1,jj+1,1) &
	133	! & + ht_i(ji-1,jj ,jl) * tmask(ji-1,jj ,1) + ht_i(ji ,jj-1,jl) * tmask(ji ,jj-1,1) &
	134	! & + ht_i(ji+1,jj ,jl) * tmask(ji+1,jj ,1) + ht_i(ji ,jj+1,jl) * tmask(ji ,jj+1,1) &
	135	! & + ht_i(ji-1,jj+1,jl) * tmask(ji-1,jj+1,1) + ht_i(ji+1,jj-1,jl) * tmask(ji+1,jj-1,1) )
	136	END DO
	137	END DO
	138	CALL lbc_lnk(zhimax(:,:,jl),'T',1.)
	139	END DO
	140
[825]	141	!-------------------------
[2715]	142	! transported fields
[825]	143	!-------------------------
[2715]	144	! Snow vol, ice vol, salt and age contents, area
	145	zs0ow(:,:) = ato_i(:,:) * area(:,:) ! Open water area
	146	DO jl = 1, jpl
	147	zs0sn (:,:,jl) = v_s (:,:,jl) * area(:,:) ! Snow volume
	148	zs0ice(:,:,jl) = v_i (:,:,jl) * area(:,:) ! Ice volume
	149	zs0a (:,:,jl) = a_i (:,:,jl) * area(:,:) ! Ice area
	150	zs0sm (:,:,jl) = smv_i(:,:,jl) * area(:,:) ! Salt content
	151	zs0oi (:,:,jl) = oa_i (:,:,jl) * area(:,:) ! Age content
	152	zs0c0 (:,:,jl) = e_s (:,:,1,jl) ! Snow heat content
	153	zs0e (:,:,:,jl) = e_i (:,:,:,jl) ! Ice heat content
[825]	154	END DO
	155
[2715]	156	!--------------------------
	157	! Advection of Ice fields (Prather scheme)
	158	!--------------------------
[825]	159	! If ice drift field is too fast, use an appropriate time step for advection.
[2715]	160	! CFL test for stability
	161	zcfl = MAXVAL( ABS( u_ice(:,:) ) * rdt_ice / e1u(:,:) )
	162	zcfl = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice / e2v(:,:) ) )
	163	IF(lk_mpp ) CALL mpp_max( zcfl )
	164	!!gm more readability:
	165	! IF( zcfl > 0.5 ) THEN ; initad = 2 ; zusnit = 0.5_wp
	166	! ELSE ; initad = 1 ; zusnit = 1.0_wp
	167	! ENDIF
	168	!!gm end
[4161]	169	initad = 1 + NINT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )
[825]	170	zusnit = 1.0 / REAL( initad )
[2715]	171	IF( zcfl > 0.5 .AND. lwp ) &
[3625]	172	WRITE(numout,*) 'lim_trp : CFL violation at day ', nday, ', cfl = ', zcfl, &
[2715]	173	& ': the ice time stepping is split in two'
[921]	174
[2715]	175	IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==!
[825]	176	DO jk = 1,initad
[2715]	177	CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area
	178	& sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) )
	179	CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:), &
	180	& sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) )
[825]	181	DO jl = 1, jpl
[2715]	182	CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume ---
	183	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
	184	CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl), &
	185	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
	186	CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---
	187	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
	188	CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), &
	189	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
	190	CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---
	191	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
	192	CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), &
	193	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
	194	CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---
	195	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
	196	CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl), &
	197	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
	198	CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---
	199	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
	200	CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0a (:,:,jl), sxa (:,:,jl), &
	201	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
	202	CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents ---
	203	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
	204	CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), &
	205	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
	206	DO layer = 1, nlay_i !--- ice heat contents ---
	207	CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &
	208	& sxxe(:,:,layer,jl), sye (:,:,layer,jl), &
	209	& syye(:,:,layer,jl), sxye(:,:,layer,jl) )
	210	CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &
	211	& sxxe(:,:,layer,jl), sye (:,:,layer,jl), &
	212	& syye(:,:,layer,jl), sxye(:,:,layer,jl) )
[825]	213	END DO
	214	END DO
	215	END DO
	216	ELSE
	217	DO jk = 1, initad
[3625]	218	CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area
[2715]	219	& sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) )
[3625]	220	CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:), &
[2715]	221	& sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) )
[825]	222	DO jl = 1, jpl
[3625]	223	CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume ---
[2715]	224	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
[3625]	225	CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl), &
[2715]	226	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
[3625]	227	CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---
[2715]	228	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
[3625]	229	CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), &
[2715]	230	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
[3625]	231	CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---
[2715]	232	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
[3625]	233	CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), &
[2715]	234	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
	235
[3625]	236	CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---
[2715]	237	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
[3625]	238	CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl), &
[2715]	239	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
[3625]	240	CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---
[2715]	241	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
[3625]	242	CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0a (:,:,jl), sxa (:,:,jl), &
[2715]	243	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
[3625]	244	CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents ---
[2715]	245	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
[3625]	246	CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), &
[2715]	247	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
	248	DO layer = 1, nlay_i !--- ice heat contents ---
[3625]	249	CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &
[2715]	250	& sxxe(:,:,layer,jl), sye (:,:,layer,jl), &
	251	& syye(:,:,layer,jl), sxye(:,:,layer,jl) )
[3625]	252	CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &
[2715]	253	& sxxe(:,:,layer,jl), sye (:,:,layer,jl), &
	254	& syye(:,:,layer,jl), sxye(:,:,layer,jl) )
[825]	255	END DO
	256	END DO
	257	END DO
	258	ENDIF
	259
	260	!-------------------------------------------
	261	! Recover the properties from their contents
	262	!-------------------------------------------
[2715]	263	zs0ow(:,:) = zs0ow(:,:) / area(:,:)
[825]	264	DO jl = 1, jpl
	265	zs0ice(:,:,jl) = zs0ice(:,:,jl) / area(:,:)
	266	zs0sn (:,:,jl) = zs0sn (:,:,jl) / area(:,:)
	267	zs0sm (:,:,jl) = zs0sm (:,:,jl) / area(:,:)
	268	zs0oi (:,:,jl) = zs0oi (:,:,jl) / area(:,:)
	269	zs0a (:,:,jl) = zs0a (:,:,jl) / area(:,:)
	270	zs0c0 (:,:,jl) = zs0c0 (:,:,jl) / area(:,:)
	271	DO jk = 1, nlay_i
	272	zs0e(:,:,jk,jl) = zs0e(:,:,jk,jl) / area(:,:)
	273	END DO
	274	END DO
	275
[921]	276	!------------------------------------------------------------------------------!
	277	! 4) Diffusion of Ice fields
	278	!------------------------------------------------------------------------------!
[825]	279
[2715]	280	!--------------------------------
	281	! diffusion of open water area
	282	!--------------------------------
	283	zs0at(:,:) = zs0a(:,:,1) ! total ice fraction
	284	DO jl = 2, jpl
	285	zs0at(:,:) = zs0at(:,:) + zs0a(:,:,jl)
	286	END DO
	287	!
	288	! ! Masked eddy diffusivity coefficient at ocean U- and V-points
	289	DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row
	290	DO ji = 1 , fs_jpim1 ! vector opt.
	291	pahu(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji ,jj) ) ) ) &
	292	& * ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji+1,jj) ) ) ) * ahiu(ji,jj)
	293	pahv(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji,jj ) ) ) ) &
	294	& * ( 1._wp - MAX( rzero, SIGN( rone,- zs0at(ji,jj+1) ) ) ) * ahiv(ji,jj)
	295	END DO
	296	END DO
	297	!
	298	CALL lim_hdf( zs0ow (:,:) ) ! Diffusion
	299
[921]	300	!------------------------------------
[2715]	301	! Diffusion of other ice variables
[921]	302	!------------------------------------
[825]	303	DO jl = 1, jpl
[2715]	304	! ! Masked eddy diffusivity coefficient at ocean U- and V-points
	305	DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row
	306	DO ji = 1 , fs_jpim1 ! vector opt.
	307	pahu(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji ,jj,jl) ) ) ) &
	308	& * ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)
	309	pahv(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji,jj ,jl) ) ) ) &
	310	& * ( 1._wp - MAX( rzero, SIGN( rone,- zs0a(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)
[825]	311	END DO
[921]	312	END DO
[825]	313
	314	CALL lim_hdf( zs0ice (:,:,jl) )
	315	CALL lim_hdf( zs0sn (:,:,jl) )
	316	CALL lim_hdf( zs0sm (:,:,jl) )
	317	CALL lim_hdf( zs0oi (:,:,jl) )
	318	CALL lim_hdf( zs0a (:,:,jl) )
	319	CALL lim_hdf( zs0c0 (:,:,jl) )
	320	DO jk = 1, nlay_i
	321	CALL lim_hdf( zs0e (:,:,jk,jl) )
[2715]	322	END DO
	323	END DO
[825]	324
[921]	325	!------------------------------------------------------------------------------!
	326	! 5) Update and limit ice properties after transport
	327	!------------------------------------------------------------------------------!
[825]	328
[921]	329	!--------------------------------------------------
	330	! 5.1) Recover mean values over the grid squares.
	331	!--------------------------------------------------
[2715]	332	zs0at(:,:) = 0._wp
[825]	333	DO jl = 1, jpl
	334	DO jj = 1, jpj
	335	DO ji = 1, jpi
[4161]	336	zs0sn (ji,jj,jl) = MAX( rzero, zs0sn (ji,jj,jl) )
	337	zs0ice(ji,jj,jl) = MAX( rzero, zs0ice(ji,jj,jl) )
	338	zs0sm (ji,jj,jl) = MAX( rzero, zs0sm (ji,jj,jl) )
	339	zs0oi (ji,jj,jl) = MAX( rzero, zs0oi (ji,jj,jl) )
	340	zs0a (ji,jj,jl) = MAX( rzero, zs0a (ji,jj,jl) )
	341	zs0c0 (ji,jj,jl) = MAX( rzero, zs0c0 (ji,jj,jl) )
[825]	342	zs0at (ji,jj) = zs0at(ji,jj) + zs0a(ji,jj,jl)
	343	END DO
	344	END DO
	345	END DO
	346
[921]	347	!---------------------------------------------------------
	348	! 5.2) Snow thickness, Ice thickness, Ice concentrations
	349	!---------------------------------------------------------
[825]	350	DO jj = 1, jpj
	351	DO ji = 1, jpi
[4161]	352	zindb = MAX( 0._wp , SIGN( 1.0, zs0at(ji,jj) - epsi10) )
[2715]	353	zs0ow(ji,jj) = ( 1._wp - zindb ) + zindb * MAX( zs0ow(ji,jj), 0._wp )
	354	ato_i(ji,jj) = zs0ow(ji,jj)
[825]	355	END DO
	356	END DO
	357
[2715]	358	DO jl = 1, jpl ! Remove very small areas
[825]	359	DO jj = 1, jpj
	360	DO ji = 1, jpi
[4161]	361	zvi = zs0ice(ji,jj,jl)
	362	zvs = zs0sn(ji,jj,jl)
[2715]	363	!
[4161]	364	zindb = MAX( 0.0 , SIGN( 1.0, zs0a(ji,jj,jl) - epsi10) )
	365	!
[2715]	366	v_s(ji,jj,jl) = zindb * zs0sn (ji,jj,jl)
	367	v_i(ji,jj,jl) = zindb * zs0ice(ji,jj,jl)
	368	!
[4161]	369	zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) )
	370	zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )
[2715]	371	zindb = MAX( zindsn, zindic )
[4161]	372	!
[2715]	373	zs0a(ji,jj,jl) = zindb * zs0a(ji,jj,jl) !ice concentration
	374	a_i (ji,jj,jl) = zs0a(ji,jj,jl)
	375	v_s (ji,jj,jl) = zindsn * v_s(ji,jj,jl)
	376	v_i (ji,jj,jl) = zindic * v_i(ji,jj,jl)
[4161]	377	!
	378	! Update mass fluxes (clem)
	379	rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic
	380	rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn
	381	END DO
	382	END DO
	383	END DO
	384
	385	!--- Thickness correction in case too high (clem) --------------------------------------------------------
	386	CALL lim_var_glo2eqv
	387	DO jl = 1, jpl
	388	DO jj = 1, jpj
	389	DO ji = 1, jpi
	390
	391	IF ( v_i(ji,jj,jl) > 0._wp ) THEN
	392	zvi = v_i(ji,jj,jl)
	393	zvs = v_s(ji,jj,jl)
	394	zdv = v_i(ji,jj,jl) - zviold(ji,jj,jl)
	395	!zda = a_i(ji,jj,jl) - zaiold(ji,jj,jl)
	396
	397	zindh = 1._wp
	398	IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. SUM( zaiold(ji,jj,1:jpl) ) < 0.80 ) .OR. &
	399	& ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN
	400	ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) )
	401	zindh = MAX( rzero, SIGN( rone, ht_i(ji,jj,jl) - epsi10 ) )
	402	a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi10 )
	403	ELSE
	404	ht_i(ji,jj,jl) = MAX( MIN( ht_i(ji,jj,jl), hi_max(jl) ), hi_max(jl-1) )
	405	zindh = MAX( rzero, SIGN( rone, ht_i(ji,jj,jl) - epsi10 ) )
	406	a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi10 )
	407	ENDIF
	408
	409	! small correction due to *zindh for a_i
	410	v_i(ji,jj,jl) = zindh * v_i(ji,jj,jl)
	411	v_s(ji,jj,jl) = zindh * v_s(ji,jj,jl)
	412
	413	! Update mass fluxes
	414	rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic
	415	rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn
	416
	417	ENDIF
	418
	419	diag_trp_vi(ji,jj) = diag_trp_vi(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice
	420
[825]	421	END DO
	422	END DO
	423	END DO
	424
[4161]	425	! ---
[825]	426	DO jj = 1, jpj
	427	DO ji = 1, jpi
[4161]	428	zs0at(ji,jj) = SUM( zs0a(ji,jj,1:jpl) ) ! clem@useless??
[825]	429	END DO
	430	END DO
	431
[921]	432	!----------------------
	433	! 5.3) Ice properties
	434	!----------------------
[825]	435
[4161]	436	zbigval = 1.e+13
[825]	437
	438	DO jl = 1, jpl
	439	DO jj = 1, jpj
	440	DO ji = 1, jpi
[4161]	441	zsmv = zs0sm(ji,jj,jl)
[825]	442
	443	! Switches and dummy variables
[4333]	444	zusvosn = 1.0/MAX( v_s(ji,jj,jl) , epsi10 )
	445	zusvoic = 1.0/MAX( v_i(ji,jj,jl) , epsi10 )
[4161]	446	zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) )
	447	zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )
[825]	448	zindb = MAX( zindsn, zindic )
	449
	450	! Ice salinity and age
[4161]	451	!clem zsal = MAX( MIN( (rhoic-rhosn)/rhoicsss_m(ji,jj), zusvoic zs0sm(ji,jj,jl) ), s_i_min ) * v_i(ji,jj,jl)
	452	IF( num_sal == 2 ) THEN
	453	smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), zsmv ), s_i_min * v_i(ji,jj,jl) )
	454	ENDIF
[825]	455
[4333]	456	zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) ), 0._wp ) * a_i(ji,jj,jl)
[3625]	457	oa_i (ji,jj,jl) = zindic * zage
[825]	458
	459	! Snow heat content
	460	ze = MIN( MAX( 0.0, zs0c0(ji,jj,jl)*area(ji,jj) ), zbigval )
[4161]	461	e_s(ji,jj,1,jl) = zindsn * ze
[825]	462
[4161]	463	! Update salt fluxes (clem)
	464	sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice
[825]	465	END DO !ji
	466	END DO !jj
	467	END DO ! jl
	468
	469	DO jl = 1, jpl
	470	DO jk = 1, nlay_i
	471	DO jj = 1, jpj
	472	DO ji = 1, jpi
	473	! Ice heat content
[4161]	474	zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )
[825]	475	ze = MIN( MAX( 0.0, zs0e(ji,jj,jk,jl)*area(ji,jj) ), zbigval )
[4161]	476	e_i(ji,jj,jk,jl) = zindic * ze
[825]	477	END DO !ji
	478	END DO ! jj
	479	END DO ! jk
	480	END DO ! jl
	481
[4161]	482
	483	! --- agglomerate variables (clem) -----------------
	484	vt_i (:,:) = 0._wp
	485	vt_s (:,:) = 0._wp
	486	at_i (:,:) = 0._wp
	487	!
	488	DO jl = 1, jpl
	489	DO jj = 1, jpj
	490	DO ji = 1, jpi
	491	!
	492	vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume
	493	vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume
	494	at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration
	495	!
[4333]	496	zinda = MAX( rzero , SIGN( rone , at_i(ji,jj) - epsi10 ) )
	497	icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * zinda ! ice thickness
[4161]	498	END DO
	499	END DO
	500	END DO
	501	! -------------------------------------------------
	502
	503
	504
[825]	505	ENDIF
	506
[863]	507	IF(ln_ctl) THEN ! Control print
[867]	508	CALL prt_ctl_info(' ')
	509	CALL prt_ctl_info(' - Cell values : ')
	510	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
[863]	511	CALL prt_ctl(tab2d_1=area , clinfo1=' lim_trp : cell area :')
	512	CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_trp : at_i :')
	513	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_trp : vt_i :')
	514	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_trp : vt_s :')
	515	DO jl = 1, jpl
[867]	516	CALL prt_ctl_info(' ')
[863]	517	CALL prt_ctl_info(' - Category : ', ivar1=jl)
	518	CALL prt_ctl_info(' ~~~~~~~~~~')
	519	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_trp : a_i : ')
	520	CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_trp : ht_i : ')
	521	CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_trp : ht_s : ')
	522	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_trp : v_i : ')
	523	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_trp : v_s : ')
	524	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_trp : e_s : ')
	525	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_trp : t_su : ')
	526	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_trp : t_snow : ')
	527	CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_trp : sm_i : ')
	528	CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_trp : smv_i : ')
	529	DO jk = 1, nlay_i
[867]	530	CALL prt_ctl_info(' ')
[863]	531	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
	532	CALL prt_ctl_info(' ~~~~~~~')
	533	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_trp : t_i : ')
	534	CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_trp : e_i : ')
	535	END DO
	536	END DO
	537	ENDIF
[4161]	538	! -------------------------------
	539	!- check conservation (C Rousset)
	540	IF( ln_limdiahsb ) THEN
	541	zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
	542	zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b
	543
	544	zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice
	545	zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic )
	546
	547	zchk_vmin = glob_min(v_i)
	548	zchk_amax = glob_max(SUM(a_i,dim=3))
	549	zchk_amin = glob_min(a_i)
	550	zchk_umax = glob_max(SQRT(u_ice2 + v_ice2))
	551
	552	IF(lwp) THEN
	553	IF ( ABS( zchk_v_i ) > 1.e-5 ) THEN
	554	WRITE(numout,) 'violation volume [m3/day] (limtrp) = ',(zchk_v_i rday)
	555	WRITE(numout,*) 'u_ice max [m/s] (limtrp) = ',zchk_umax
	556	WRITE(numout,*) 'number of time steps (limtrp) =',kt
	557	ENDIF
	558	IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,) 'violation saline [psum3/day] (limtrp) = ',(zchk_smv * rday)
	559	IF ( zchk_vmin < 0. ) WRITE(numout,) 'violation v_i<0 [mm] (limtrp) = ',(zchk_vmin 1.e-3)
	560	IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limtrp) = ',zchk_amin
	561	ENDIF
	562	ENDIF
	563	!- check conservation (C Rousset)
	564	! -------------------------------
[2715]	565	!
[3294]	566	CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow )
	567	CALL wrk_dealloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi )
	568	CALL wrk_dealloc( jpi, jpj, jkmax, jpl, zs0e )
[4161]	569
	570	CALL wrk_dealloc( jpi,jpj,jpl,zaiold, zhimax ) ! clem
[2715]	571	!
[4161]	572	IF( nn_timing == 1 ) CALL timing_stop('limtrp')
[825]	573	END SUBROUTINE lim_trp
	574
	575	#else
	576	!!----------------------------------------------------------------------
	577	!! Default option Empty Module No sea-ice model
	578	!!----------------------------------------------------------------------
	579	CONTAINS
	580	SUBROUTINE lim_trp ! Empty routine
	581	END SUBROUTINE lim_trp
	582	#endif
	583
	584	!!======================================================================
	585	END MODULE limtrp

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

Download in other formats: