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 trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 4873

Last change on this file since 4873 was 4873, checked in by clem, 9 years ago
LIM3: removing useless jkmax (=nlay_i+1) parameter
Property svn:keywords set to `Id`
File size: 29.6 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
[4688]	32	USE limcons ! conservation tests
[825]	33
	34	IMPLICIT NONE
	35	PRIVATE
	36
[2715]	37	PUBLIC lim_trp ! called by ice_step
[825]	38
[4161]	39	REAL(wp) :: epsi10 = 1.e-10_wp
[4688]	40	REAL(wp) :: epsi20 = 1.e-20_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	!
[4870]	65	INTEGER :: ji, jj, jk, jl ! dummy loop indices
[2715]	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
[4205]	69	REAL(wp) :: zcfl , zusnit ! - -
[4688]	70	REAL(wp) :: zsal , zage ! - -
[2715]	71	!
[3294]	72	REAL(wp), POINTER, DIMENSION(:,:) :: zui_u, zvi_v, zsm, zs0at, zs0ow
	73	REAL(wp), POINTER, DIMENSION(:,:,:) :: zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi
	74	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zs0e
[4161]	75	! mass and salt flux (clem)
[4688]	76	REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold ! old ice volume...
[4161]	77	REAL(wp), POINTER, DIMENSION(:,:,:) :: zaiold, zhimax ! old ice concentration and thickness
[4688]	78	REAL(wp), POINTER, DIMENSION(:,:) :: zeiold, zesold ! old enthalpies
	79	REAL(wp) :: zdv, zda, zvi, zvs, zsmv, zes, zei
	80	!
	81	REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
[2715]	82	!!---------------------------------------------------------------------
[4161]	83	IF( nn_timing == 1 ) CALL timing_start('limtrp')
[825]	84
[4688]	85	CALL wrk_alloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold )
[3294]	86	CALL wrk_alloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi )
[4873]	87	CALL wrk_alloc( jpi, jpj, nlay_i+1, jpl, zs0e )
[825]	88
[4688]	89	CALL wrk_alloc( jpi, jpj, jpl, zaiold, zhimax, zviold, zvsold ) ! clem
[4161]	90
[2715]	91	IF( numit == nstart .AND. lwp ) THEN
	92	WRITE(numout,*)
	93	IF( ln_limdyn ) THEN ; WRITE(numout,*) 'lim_trp : Ice transport '
	94	ELSE ; WRITE(numout,*) 'lim_trp : No ice advection as ln_limdyn = ', ln_limdyn
	95	ENDIF
	96	WRITE(numout,*) '~~~~~~~~~~~~'
	97	ENDIF
	98
[825]	99	zsm(:,:) = area(:,:)
	100
[2715]	101	! !-------------------------------------!
	102	IF( ln_limdyn ) THEN ! Advection of sea ice properties !
	103	! !-------------------------------------!
[4688]	104
	105	! conservation test
	106	IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
	107
[4161]	108	! mass and salt flux init (clem)
	109	zviold(:,:,:) = v_i(:,:,:)
[4688]	110	zeiold(:,:) = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )
	111	zesold(:,:) = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )
[825]	112
[4161]	113	!--- Thickness correction init. (clem) -------------------------------
	114	CALL lim_var_glo2eqv
	115	zaiold(:,:,:) = a_i(:,:,:)
	116	!---------------------------------------------------------------------
	117	! Record max of the surrounding ice thicknesses for correction in limupdate
	118	! in case advection creates ice too thick.
	119	!---------------------------------------------------------------------
	120	zhimax(:,:,:) = ht_i(:,:,:)
	121	DO jl = 1, jpl
	122	DO jj = 2, jpjm1
	123	DO ji = 2, jpim1
	124	zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) )
	125	!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) &
	126	! & + ht_i(ji-1,jj ,jl) * tmask(ji-1,jj ,1) + ht_i(ji ,jj-1,jl) * tmask(ji ,jj-1,1) &
	127	! & + ht_i(ji+1,jj ,jl) * tmask(ji+1,jj ,1) + ht_i(ji ,jj+1,jl) * tmask(ji ,jj+1,1) &
	128	! & + 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) )
	129	END DO
	130	END DO
	131	CALL lbc_lnk(zhimax(:,:,jl),'T',1.)
	132	END DO
	133
[825]	134	!-------------------------
[2715]	135	! transported fields
[825]	136	!-------------------------
[2715]	137	! Snow vol, ice vol, salt and age contents, area
	138	zs0ow(:,:) = ato_i(:,:) * area(:,:) ! Open water area
	139	DO jl = 1, jpl
	140	zs0sn (:,:,jl) = v_s (:,:,jl) * area(:,:) ! Snow volume
	141	zs0ice(:,:,jl) = v_i (:,:,jl) * area(:,:) ! Ice volume
	142	zs0a (:,:,jl) = a_i (:,:,jl) * area(:,:) ! Ice area
	143	zs0sm (:,:,jl) = smv_i(:,:,jl) * area(:,:) ! Salt content
	144	zs0oi (:,:,jl) = oa_i (:,:,jl) * area(:,:) ! Age content
	145	zs0c0 (:,:,jl) = e_s (:,:,1,jl) ! Snow heat content
	146	zs0e (:,:,:,jl) = e_i (:,:,:,jl) ! Ice heat content
[825]	147	END DO
	148
[2715]	149	!--------------------------
	150	! Advection of Ice fields (Prather scheme)
	151	!--------------------------
[825]	152	! If ice drift field is too fast, use an appropriate time step for advection.
[2715]	153	! CFL test for stability
	154	zcfl = MAXVAL( ABS( u_ice(:,:) ) * rdt_ice / e1u(:,:) )
	155	zcfl = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice / e2v(:,:) ) )
	156	IF(lk_mpp ) CALL mpp_max( zcfl )
	157	!!gm more readability:
	158	! IF( zcfl > 0.5 ) THEN ; initad = 2 ; zusnit = 0.5_wp
	159	! ELSE ; initad = 1 ; zusnit = 1.0_wp
	160	! ENDIF
	161	!!gm end
[4688]	162	initad = 1 + NINT( MAX( 0._wp, SIGN( 1._wp, zcfl-0.5 ) ) )
[825]	163	zusnit = 1.0 / REAL( initad )
[2715]	164	IF( zcfl > 0.5 .AND. lwp ) &
[3625]	165	WRITE(numout,*) 'lim_trp : CFL violation at day ', nday, ', cfl = ', zcfl, &
[2715]	166	& ': the ice time stepping is split in two'
[921]	167
[2715]	168	IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==!
[825]	169	DO jk = 1,initad
[4688]	170	CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area
[2715]	171	& sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) )
[4688]	172	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ow (:,:), sxopw(:,:), &
[2715]	173	& sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) )
[825]	174	DO jl = 1, jpl
[4688]	175	CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume ---
[2715]	176	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
[4688]	177	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl), &
[2715]	178	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
[4688]	179	CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---
[2715]	180	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
[4688]	181	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), &
[2715]	182	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
[4688]	183	CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---
[2715]	184	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
[4688]	185	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), &
[2715]	186	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
[4688]	187	CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---
[2715]	188	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
[4688]	189	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl), &
[2715]	190	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
[4688]	191	CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---
[2715]	192	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
[4688]	193	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0a (:,:,jl), sxa (:,:,jl), &
[2715]	194	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
[4688]	195	CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents ---
[2715]	196	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
[4688]	197	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), &
[2715]	198	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
[4870]	199	DO jk = 1, nlay_i !--- ice heat contents ---
	200	CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), &
	201	& sxxe(:,:,jk,jl), sye (:,:,jk,jl), &
	202	& syye(:,:,jk,jl), sxye(:,:,jk,jl) )
	203	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), &
	204	& sxxe(:,:,jk,jl), sye (:,:,jk,jl), &
	205	& syye(:,:,jk,jl), sxye(:,:,jk,jl) )
[825]	206	END DO
	207	END DO
	208	END DO
	209	ELSE
	210	DO jk = 1, initad
[4688]	211	CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area
[2715]	212	& sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) )
[4688]	213	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ow (:,:), sxopw(:,:), &
[2715]	214	& sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) )
[825]	215	DO jl = 1, jpl
[4688]	216	CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume ---
[2715]	217	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
[4688]	218	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl), &
[2715]	219	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
[4688]	220	CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---
[2715]	221	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
[4688]	222	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), &
[2715]	223	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
[4688]	224	CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---
[2715]	225	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
[4688]	226	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), &
[2715]	227	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
	228
[4688]	229	CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---
[2715]	230	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
[4688]	231	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl), &
[2715]	232	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
[4688]	233	CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---
[2715]	234	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
[4688]	235	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0a (:,:,jl), sxa (:,:,jl), &
[2715]	236	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
[4688]	237	CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents ---
[2715]	238	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
[4688]	239	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), &
[2715]	240	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
[4870]	241	DO jk = 1, nlay_i !--- ice heat contents ---
	242	CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), &
	243	& sxxe(:,:,jk,jl), sye (:,:,jk,jl), &
	244	& syye(:,:,jk,jl), sxye(:,:,jk,jl) )
	245	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), &
	246	& sxxe(:,:,jk,jl), sye (:,:,jk,jl), &
	247	& syye(:,:,jk,jl), sxye(:,:,jk,jl) )
[825]	248	END DO
	249	END DO
	250	END DO
	251	ENDIF
	252
	253	!-------------------------------------------
	254	! Recover the properties from their contents
	255	!-------------------------------------------
[2715]	256	zs0ow(:,:) = zs0ow(:,:) / area(:,:)
[825]	257	DO jl = 1, jpl
	258	zs0ice(:,:,jl) = zs0ice(:,:,jl) / area(:,:)
	259	zs0sn (:,:,jl) = zs0sn (:,:,jl) / area(:,:)
	260	zs0sm (:,:,jl) = zs0sm (:,:,jl) / area(:,:)
	261	zs0oi (:,:,jl) = zs0oi (:,:,jl) / area(:,:)
	262	zs0a (:,:,jl) = zs0a (:,:,jl) / area(:,:)
[4688]	263	!
[825]	264	END DO
	265
[921]	266	!------------------------------------------------------------------------------!
	267	! 4) Diffusion of Ice fields
	268	!------------------------------------------------------------------------------!
[825]	269
[2715]	270	!--------------------------------
	271	! diffusion of open water area
	272	!--------------------------------
	273	zs0at(:,:) = zs0a(:,:,1) ! total ice fraction
	274	DO jl = 2, jpl
	275	zs0at(:,:) = zs0at(:,:) + zs0a(:,:,jl)
	276	END DO
	277	!
	278	! ! Masked eddy diffusivity coefficient at ocean U- and V-points
	279	DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row
	280	DO ji = 1 , fs_jpim1 ! vector opt.
[4688]	281	pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0at(ji ,jj) ) ) ) &
	282	& * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0at(ji+1,jj) ) ) ) * ahiu(ji,jj)
	283	pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0at(ji,jj ) ) ) ) &
	284	& * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- zs0at(ji,jj+1) ) ) ) * ahiv(ji,jj)
[2715]	285	END DO
	286	END DO
	287	!
	288	CALL lim_hdf( zs0ow (:,:) ) ! Diffusion
	289
[921]	290	!------------------------------------
[2715]	291	! Diffusion of other ice variables
[921]	292	!------------------------------------
[825]	293	DO jl = 1, jpl
[2715]	294	! ! Masked eddy diffusivity coefficient at ocean U- and V-points
	295	DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row
	296	DO ji = 1 , fs_jpim1 ! vector opt.
[4688]	297	pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0a(ji ,jj,jl) ) ) ) &
	298	& * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0a(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)
	299	pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0a(ji,jj ,jl) ) ) ) &
	300	& * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- zs0a(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)
[825]	301	END DO
[921]	302	END DO
[825]	303
	304	CALL lim_hdf( zs0ice (:,:,jl) )
	305	CALL lim_hdf( zs0sn (:,:,jl) )
	306	CALL lim_hdf( zs0sm (:,:,jl) )
	307	CALL lim_hdf( zs0oi (:,:,jl) )
	308	CALL lim_hdf( zs0a (:,:,jl) )
	309	CALL lim_hdf( zs0c0 (:,:,jl) )
	310	DO jk = 1, nlay_i
	311	CALL lim_hdf( zs0e (:,:,jk,jl) )
[2715]	312	END DO
	313	END DO
[825]	314
[921]	315	!------------------------------------------------------------------------------!
	316	! 5) Update and limit ice properties after transport
	317	!------------------------------------------------------------------------------!
[825]	318
[921]	319	!--------------------------------------------------
	320	! 5.1) Recover mean values over the grid squares.
	321	!--------------------------------------------------
[2715]	322	zs0at(:,:) = 0._wp
[825]	323	DO jl = 1, jpl
	324	DO jj = 1, jpj
	325	DO ji = 1, jpi
[4688]	326	zs0sn (ji,jj,jl) = MAX( 0._wp, zs0sn (ji,jj,jl) )
	327	zs0ice(ji,jj,jl) = MAX( 0._wp, zs0ice(ji,jj,jl) )
	328	zs0sm (ji,jj,jl) = MAX( 0._wp, zs0sm (ji,jj,jl) )
	329	zs0oi (ji,jj,jl) = MAX( 0._wp, zs0oi (ji,jj,jl) )
	330	zs0a (ji,jj,jl) = MAX( 0._wp, zs0a (ji,jj,jl) )
	331	zs0c0 (ji,jj,jl) = MAX( 0._wp, zs0c0 (ji,jj,jl) )
[825]	332	zs0at (ji,jj) = zs0at(ji,jj) + zs0a(ji,jj,jl)
	333	END DO
	334	END DO
	335	END DO
	336
[921]	337	!---------------------------------------------------------
[4688]	338	! 5.2) Update and mask variables
[921]	339	!---------------------------------------------------------
[4688]	340	DO jl = 1, jpl
[825]	341	DO jj = 1, jpj
	342	DO ji = 1, jpi
[4688]	343	zindb= MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) )
	344
	345	zvi = zs0ice(ji,jj,jl)
	346	zvs = zs0sn (ji,jj,jl)
	347	zes = zs0c0 (ji,jj,jl)
	348	zsmv = zs0sm (ji,jj,jl)
[2715]	349	!
[4688]	350	! Remove very small areas
	351	v_s(ji,jj,jl) = zindb * zs0sn (ji,jj,jl)
	352	v_i(ji,jj,jl) = zindb * zs0ice(ji,jj,jl)
	353	a_i(ji,jj,jl) = zindb * zs0a (ji,jj,jl)
	354	e_s(ji,jj,1,jl) = zindb * zs0c0 (ji,jj,jl)
	355	! Ice salinity and age
	356	IF( num_sal == 2 ) THEN
	357	smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), zsmv ), s_i_min * v_i(ji,jj,jl) )
	358	ENDIF
	359	oa_i(ji,jj,jl) = MAX( zindb * zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ), 0._wp ) * a_i(ji,jj,jl)
	360
	361	! Update fluxes
	362	wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice
	363	wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
	364	sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice
	365	hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
[4161]	366	END DO
	367	END DO
	368	END DO
	369
[4688]	370	DO jl = 1, jpl
	371	DO jk = 1, nlay_i
	372	DO jj = 1, jpj
	373	DO ji = 1, jpi
	374	zindb = MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) )
	375	zei = zs0e(ji,jj,jk,jl)
	376	e_i(ji,jj,jk,jl) = zindb * MAX( 0._wp, zs0e(ji,jj,jk,jl) )
	377	! Update fluxes
	378	hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
	379	END DO !ji
	380	END DO ! jj
	381	END DO ! jk
	382	END DO ! jl
	383
[4161]	384	!--- Thickness correction in case too high (clem) --------------------------------------------------------
	385	CALL lim_var_glo2eqv
	386	DO jl = 1, jpl
	387	DO jj = 1, jpj
	388	DO ji = 1, jpi
	389
	390	IF ( v_i(ji,jj,jl) > 0._wp ) THEN
[4688]	391	zvi = v_i (ji,jj,jl)
	392	zvs = v_s (ji,jj,jl)
	393	zsmv = smv_i(ji,jj,jl)
	394	zes = e_s (ji,jj,1,jl)
	395	zei = SUM( e_i(ji,jj,:,jl) )
	396	zdv = v_i(ji,jj,jl) - zviold(ji,jj,jl)
[4161]	397	!zda = a_i(ji,jj,jl) - zaiold(ji,jj,jl)
	398
	399	zindh = 1._wp
	400	IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. SUM( zaiold(ji,jj,1:jpl) ) < 0.80 ) .OR. &
	401	& ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN
	402	ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) )
[4688]	403	zindh = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) )
	404	a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 )
[4161]	405	ELSE
	406	ht_i(ji,jj,jl) = MAX( MIN( ht_i(ji,jj,jl), hi_max(jl) ), hi_max(jl-1) )
[4688]	407	zindh = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) )
	408	a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 )
[4161]	409	ENDIF
	410
	411	! small correction due to *zindh for a_i
[4688]	412	v_i (ji,jj,jl) = zindh * v_i (ji,jj,jl)
	413	v_s (ji,jj,jl) = zindh * v_s (ji,jj,jl)
	414	smv_i(ji,jj,jl) = zindh * smv_i(ji,jj,jl)
	415	e_s(ji,jj,1,jl) = zindh * e_s(ji,jj,1,jl)
	416	e_i(ji,jj,:,jl) = zindh * e_i(ji,jj,:,jl)
[4161]	417
	418	! Update mass fluxes
[4688]	419	wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice
	420	wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
	421	sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice
	422	hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
	423	hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,:,jl) ) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
[4161]	424
	425	ENDIF
	426
	427	diag_trp_vi(ji,jj) = diag_trp_vi(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice
[4688]	428	diag_trp_vs(ji,jj) = diag_trp_vs(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * r1_rdtice
[4161]	429
[825]	430	END DO
	431	END DO
	432	END DO
[4688]	433	! -------------------------------------------------
[825]	434
[4688]	435	! --- diags ---
[825]	436	DO jj = 1, jpj
	437	DO ji = 1, jpi
[4688]	438	diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice
	439	diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice
[825]	440	END DO
	441	END DO
	442
[4688]	443	! --- agglomerate variables (clem) -----------------
	444	vt_i (:,:) = 0._wp
	445	vt_s (:,:) = 0._wp
	446	at_i (:,:) = 0._wp
	447	!
[825]	448	DO jl = 1, jpl
	449	DO jj = 1, jpj
	450	DO ji = 1, jpi
[4688]	451	!
	452	vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume
	453	vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume
	454	at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration
	455	END DO
	456	END DO
	457	END DO
	458	! -------------------------------------------------
[825]	459
[4688]	460	! open water
[4161]	461	DO jj = 1, jpj
	462	DO ji = 1, jpi
[4688]	463	! open water = 1 if at_i=0
	464	zindb = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
	465	ato_i(ji,jj) = zindb + (1._wp - zindb ) * zs0ow(ji,jj)
[4161]	466	END DO
[4688]	467	END DO
[4161]	468
[4688]	469	! conservation test
	470	IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[4161]	471
[825]	472	ENDIF
	473
[863]	474	IF(ln_ctl) THEN ! Control print
[867]	475	CALL prt_ctl_info(' ')
	476	CALL prt_ctl_info(' - Cell values : ')
	477	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
[863]	478	CALL prt_ctl(tab2d_1=area , clinfo1=' lim_trp : cell area :')
	479	CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_trp : at_i :')
	480	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_trp : vt_i :')
	481	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_trp : vt_s :')
	482	DO jl = 1, jpl
[867]	483	CALL prt_ctl_info(' ')
[863]	484	CALL prt_ctl_info(' - Category : ', ivar1=jl)
	485	CALL prt_ctl_info(' ~~~~~~~~~~')
	486	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_trp : a_i : ')
	487	CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_trp : ht_i : ')
	488	CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_trp : ht_s : ')
	489	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_trp : v_i : ')
	490	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_trp : v_s : ')
	491	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_trp : e_s : ')
	492	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_trp : t_su : ')
	493	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_trp : t_snow : ')
	494	CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_trp : sm_i : ')
	495	CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_trp : smv_i : ')
	496	DO jk = 1, nlay_i
[867]	497	CALL prt_ctl_info(' ')
[863]	498	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
	499	CALL prt_ctl_info(' ~~~~~~~')
	500	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_trp : t_i : ')
	501	CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_trp : e_i : ')
	502	END DO
	503	END DO
	504	ENDIF
[2715]	505	!
[4688]	506	CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold )
[3294]	507	CALL wrk_dealloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi )
[4873]	508	CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, zs0e )
[4161]	509
[4688]	510	CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zaiold, zhimax ) ! clem
[2715]	511	!
[4161]	512	IF( nn_timing == 1 ) CALL timing_stop('limtrp')
[825]	513	END SUBROUTINE lim_trp
	514
	515	#else
	516	!!----------------------------------------------------------------------
	517	!! Default option Empty Module No sea-ice model
	518	!!----------------------------------------------------------------------
	519	CONTAINS
	520	SUBROUTINE lim_trp ! Empty routine
	521	END SUBROUTINE lim_trp
	522	#endif
	523
	524	!!======================================================================
	525	END MODULE limtrp

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

Download in other formats: