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 @ 5147

Last change on this file since 5147 was 5138, checked in by vancop, 9 years ago
Â«Useless ctl_stop in limtrp.F90 (lim3)
Property svn:keywords set to `Id`
File size: 25.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 dom_ice ! ice domain
	20	USE ice ! ice variables
	21	USE limadv ! ice advection
	22	USE limhdf ! ice horizontal diffusion
[5123]	23	USE limvar !
	24	!
[3625]	25	USE in_out_manager ! I/O manager
	26	USE lbclnk ! lateral boundary conditions -- MPP exchanges
	27	USE lib_mpp ! MPP library
	28	USE wrk_nemo ! work arrays
	29	USE prtctl ! Print control
	30	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[5123]	31	USE timing ! Timing
[4688]	32	USE limcons ! conservation tests
[5123]	33	USE limctl ! control prints
[825]	34
	35	IMPLICIT NONE
	36	PRIVATE
	37
[5123]	38	PUBLIC lim_trp ! called by sbcice_lim
[825]	39
[5123]	40	INTEGER :: ncfl ! number of ice time step with CFL>1/2
	41
[825]	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	!!---------------------------------------------------------------------
[5123]	63	INTEGER, INTENT(in) :: kt ! number of iteration
[2715]	64	!
[5123]	65	INTEGER :: ji, jj, jk, jl, jt ! dummy loop indices
[2715]	66	INTEGER :: initad ! number of sub-timestep for the advection
[4990]	67	REAL(wp) :: zcfl , zusnit ! - -
[5123]	68	CHARACTER(len=80) :: cltmp
[2715]	69	!
[5134]	70	REAL(wp), POINTER, DIMENSION(:,:) :: zsm
	71	REAL(wp), POINTER, DIMENSION(:,:,:) :: z0ice, z0snw, z0ai, z0es , z0smi , z0oi
	72	REAL(wp), POINTER, DIMENSION(:,:,:) :: z0opw
	73	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: z0ei
[5123]	74	REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold, zsmvold ! old ice volume...
	75	REAL(wp), POINTER, DIMENSION(:,:,:) :: zhimax ! old ice thickness
	76	REAL(wp), POINTER, DIMENSION(:,:) :: zatold, zeiold, zesold ! old concentration, enthalpies
	77	REAL(wp) :: zdv, zvi, zvs, zsmv, zes, zei
	78	REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
[2715]	79	!!---------------------------------------------------------------------
[4161]	80	IF( nn_timing == 1 ) CALL timing_start('limtrp')
[825]	81
[5134]	82	CALL wrk_alloc( jpi,jpj, zsm, zatold, zeiold, zesold )
	83	CALL wrk_alloc( jpi,jpj,jpl, z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
	84	CALL wrk_alloc( jpi,jpj,1, z0opw )
	85	CALL wrk_alloc( jpi,jpj,nlay_i+1,jpl, z0ei )
[5123]	86	CALL wrk_alloc( jpi,jpj,jpl, zhimax, zviold, zvsold, zsmvold )
[825]	87
[2715]	88	IF( numit == nstart .AND. lwp ) THEN
	89	WRITE(numout,*)
	90	IF( ln_limdyn ) THEN ; WRITE(numout,*) 'lim_trp : Ice transport '
	91	ELSE ; WRITE(numout,*) 'lim_trp : No ice advection as ln_limdyn = ', ln_limdyn
	92	ENDIF
	93	WRITE(numout,*) '~~~~~~~~~~~~'
[5123]	94	ncfl = 0 ! nb of time step with CFL > 1/2
[2715]	95	ENDIF
[5123]	96
	97	zsm(:,:) = e12t(:,:)
[2715]	98
	99	! !-------------------------------------!
	100	IF( ln_limdyn ) THEN ! Advection of sea ice properties !
	101	! !-------------------------------------!
[4688]	102
	103	! conservation test
[5123]	104	IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[4688]	105
[5123]	106	! mass and salt flux init
[4161]	107	zviold(:,:,:) = v_i(:,:,:)
[5123]	108	zvsold(:,:,:) = v_s(:,:,:)
	109	zsmvold(:,:,:) = smv_i(:,:,:)
	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
[5123]	113	!--- Thickness correction init. -------------------------------
[4161]	114	CALL lim_var_glo2eqv
[5123]	115	zatold(:,:) = SUM( a_i(:,:,:), dim=3 )
[4161]	116	!---------------------------------------------------------------------
	117	! Record max of the surrounding ice thicknesses for correction in limupdate
	118	! in case advection creates ice too thick.
	119	!---------------------------------------------------------------------
[5134]	120	zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:)
[4161]	121	DO jl = 1, jpl
	122	DO jj = 2, jpjm1
	123	DO ji = 2, jpim1
[5134]	124	zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) )
[4161]	125	END DO
	126	END DO
	127	CALL lbc_lnk(zhimax(:,:,jl),'T',1.)
	128	END DO
	129
[5123]	130	!=============================!
	131	!== Prather scheme ==!
	132	!=============================!
	133
	134	! If ice drift field is too fast, use an appropriate time step for advection.
	135	zcfl = MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) ! CFL test for stability
	136	zcfl = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
	137	IF(lk_mpp ) CALL mpp_max( zcfl )
	138
	139	IF( zcfl > 0.5 ) THEN ; initad = 2 ; zusnit = 0.5_wp
	140	ELSE ; initad = 1 ; zusnit = 1.0_wp
	141	ENDIF
	142
	143	IF( zcfl > 0.5_wp .AND. lwp ) ncfl = ncfl + 1
[5138]	144	IF( lwp ) THEN
[5123]	145	IF( ncfl > 0 ) THEN
	146	WRITE(cltmp,'(i6.1)') ncfl
[5138]	147	CALL ctl_warn( 'lim_trp: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ')
[5123]	148	ELSE
[5138]	149	WRITE(numout,*) 'lim_trp : CFL criterion for ice advection is always smaller than 1/2 '
[5123]	150	ENDIF
	151	ENDIF
	152
[825]	153	!-------------------------
[2715]	154	! transported fields
[825]	155	!-------------------------
[5134]	156	z0opw(:,:,1) = ato_i(:,:) * e12t(:,:) ! Open water area
[2715]	157	DO jl = 1, jpl
[5134]	158	z0snw (:,:,jl) = v_s (:,:,jl) * e12t(:,:) ! Snow volume
	159	z0ice(:,:,jl) = v_i (:,:,jl) * e12t(:,:) ! Ice volume
	160	z0ai (:,:,jl) = a_i (:,:,jl) * e12t(:,:) ! Ice area
	161	z0smi (:,:,jl) = smv_i(:,:,jl) * e12t(:,:) ! Salt content
	162	z0oi (:,:,jl) = oa_i (:,:,jl) * e12t(:,:) ! Age content
	163	z0es (:,:,jl) = e_s (:,:,1,jl) * e12t(:,:) ! Snow heat content
[5123]	164	DO jk = 1, nlay_i
[5134]	165	z0ei (:,:,jk,jl) = e_i (:,:,jk,jl) * e12t(:,:) ! Ice heat content
[5123]	166	END DO
[825]	167	END DO
	168
[921]	169
[2715]	170	IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==!
[5123]	171	DO jt = 1, initad
[5134]	172	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area
[5123]	173	& sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )
[5134]	174	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:), &
[5123]	175	& sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )
[825]	176	DO jl = 1, jpl
[5134]	177	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume ---
[2715]	178	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
[5134]	179	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl), &
[2715]	180	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
[5134]	181	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---
[2715]	182	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
[5134]	183	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl), &
[2715]	184	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
[5134]	185	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---
[2715]	186	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
[5134]	187	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl), &
[2715]	188	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
[5134]	189	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---
[2715]	190	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
[5134]	191	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), &
[2715]	192	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
[5134]	193	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---
[2715]	194	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
[5134]	195	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ai (:,:,jl), sxa (:,:,jl), &
[2715]	196	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
[5134]	197	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents ---
[2715]	198	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
[5134]	199	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl), &
[2715]	200	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
[5123]	201	DO jk = 1, nlay_i !--- ice heat contents ---
[5134]	202	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), &
[4870]	203	& sxxe(:,:,jk,jl), sye (:,:,jk,jl), &
	204	& syye(:,:,jk,jl), sxye(:,:,jk,jl) )
[5134]	205	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), &
[4870]	206	& sxxe(:,:,jk,jl), sye (:,:,jk,jl), &
	207	& syye(:,:,jk,jl), sxye(:,:,jk,jl) )
[825]	208	END DO
	209	END DO
	210	END DO
	211	ELSE
[5123]	212	DO jt = 1, initad
[5134]	213	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area
[5123]	214	& sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )
[5134]	215	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:), &
[5123]	216	& sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )
[825]	217	DO jl = 1, jpl
[5134]	218	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume ---
[2715]	219	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
[5134]	220	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl), &
[2715]	221	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
[5134]	222	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---
[2715]	223	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
[5134]	224	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl), &
[2715]	225	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
[5134]	226	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---
[2715]	227	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
[5134]	228	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl), &
[2715]	229	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
	230
[5134]	231	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---
[2715]	232	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
[5134]	233	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), &
[2715]	234	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
[5134]	235	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---
[2715]	236	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
[5134]	237	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ai (:,:,jl), sxa (:,:,jl), &
[2715]	238	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
[5134]	239	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents ---
[2715]	240	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
[5134]	241	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl), &
[2715]	242	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
[4870]	243	DO jk = 1, nlay_i !--- ice heat contents ---
[5134]	244	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), &
[4870]	245	& sxxe(:,:,jk,jl), sye (:,:,jk,jl), &
	246	& syye(:,:,jk,jl), sxye(:,:,jk,jl) )
[5134]	247	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), &
[4870]	248	& sxxe(:,:,jk,jl), sye (:,:,jk,jl), &
	249	& syye(:,:,jk,jl), sxye(:,:,jk,jl) )
[825]	250	END DO
	251	END DO
	252	END DO
	253	ENDIF
	254
	255	!-------------------------------------------
	256	! Recover the properties from their contents
	257	!-------------------------------------------
[5134]	258	ato_i(:,:) = z0opw(:,:,1) * r1_e12t(:,:)
[825]	259	DO jl = 1, jpl
[5134]	260	v_i (:,:,jl) = z0ice(:,:,jl) * r1_e12t(:,:)
	261	v_s (:,:,jl) = z0snw(:,:,jl) * r1_e12t(:,:)
	262	smv_i(:,:,jl) = z0smi(:,:,jl) * r1_e12t(:,:)
	263	oa_i (:,:,jl) = z0oi (:,:,jl) * r1_e12t(:,:)
	264	a_i (:,:,jl) = z0ai (:,:,jl) * r1_e12t(:,:)
	265	e_s (:,:,1,jl) = z0es (:,:,jl) * r1_e12t(:,:)
[5123]	266	DO jk = 1, nlay_i
[5134]	267	e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e12t(:,:)
[5123]	268	END DO
[825]	269	END DO
	270
[5123]	271	at_i(:,:) = a_i(:,:,1) ! total ice fraction
	272	DO jl = 2, jpl
	273	at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
	274	END DO
	275
[921]	276	!------------------------------------------------------------------------------!
[5123]	277	! Diffusion of Ice fields
[921]	278	!------------------------------------------------------------------------------!
[825]	279
[5123]	280	!
[2715]	281	!--------------------------------
	282	! diffusion of open water area
	283	!--------------------------------
	284	! ! Masked eddy diffusivity coefficient at ocean U- and V-points
	285	DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row
	286	DO ji = 1 , fs_jpim1 ! vector opt.
[5123]	287	pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji ,jj) ) ) ) &
	288	& * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj)
	289	pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj ) ) ) ) &
	290	& * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj)
[2715]	291	END DO
	292	END DO
	293	!
[5134]	294	CALL lim_hdf( ato_i (:,:) )
[2715]	295
[921]	296	!------------------------------------
[2715]	297	! Diffusion of other ice variables
[921]	298	!------------------------------------
[825]	299	DO jl = 1, jpl
[2715]	300	! ! Masked eddy diffusivity coefficient at ocean U- and V-points
	301	DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row
	302	DO ji = 1 , fs_jpim1 ! vector opt.
[5123]	303	pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji ,jj,jl) ) ) ) &
	304	& * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)
	305	pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,jj ,jl) ) ) ) &
	306	& * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)
[825]	307	END DO
[921]	308	END DO
[825]	309
[5123]	310	CALL lim_hdf( v_i (:,:, jl) )
	311	CALL lim_hdf( v_s (:,:, jl) )
	312	CALL lim_hdf( smv_i(:,:, jl) )
	313	CALL lim_hdf( oa_i (:,:, jl) )
	314	CALL lim_hdf( a_i (:,:, jl) )
	315	CALL lim_hdf( e_s (:,:,1,jl) )
[825]	316	DO jk = 1, nlay_i
[5123]	317	CALL lim_hdf( e_i(:,:,jk,jl) )
[2715]	318	END DO
	319	END DO
[825]	320
[921]	321	!------------------------------------------------------------------------------!
[5123]	322	! limit ice properties after transport
[921]	323	!------------------------------------------------------------------------------!
[5123]	324	!!gm & cr : MAX should not be active if adv scheme is positive !
[825]	325	DO jl = 1, jpl
	326	DO jj = 1, jpj
	327	DO ji = 1, jpi
[5123]	328	v_s (ji,jj,jl) = MAX( 0._wp, v_s (ji,jj,jl) )
	329	v_i (ji,jj,jl) = MAX( 0._wp, v_i (ji,jj,jl) )
	330	smv_i(ji,jj,jl) = MAX( 0._wp, smv_i(ji,jj,jl) )
	331	oa_i (ji,jj,jl) = MAX( 0._wp, oa_i (ji,jj,jl) )
	332	a_i (ji,jj,jl) = MAX( 0._wp, a_i (ji,jj,jl) )
	333	e_s (ji,jj,1,jl) = MAX( 0._wp, e_s (ji,jj,1,jl) )
[825]	334	END DO
	335	END DO
	336
[4688]	337	DO jk = 1, nlay_i
	338	DO jj = 1, jpj
	339	DO ji = 1, jpi
[5123]	340	e_i(ji,jj,jk,jl) = MAX( 0._wp, e_i(ji,jj,jk,jl) )
	341	END DO
	342	END DO
	343	END DO
	344	END DO
	345	!!gm & cr
[4688]	346
[5123]	347	! zap small areas
	348	CALL lim_var_zapsmall
	349
	350	!--- Thickness correction in case too high --------------------------------------------------------
[4161]	351	CALL lim_var_glo2eqv
	352	DO jl = 1, jpl
	353	DO jj = 1, jpj
	354	DO ji = 1, jpi
	355
	356	IF ( v_i(ji,jj,jl) > 0._wp ) THEN
[5134]	357
[4688]	358	zvi = v_i (ji,jj,jl)
	359	zvs = v_s (ji,jj,jl)
	360	zsmv = smv_i(ji,jj,jl)
	361	zes = e_s (ji,jj,1,jl)
[4990]	362	zei = SUM( e_i(ji,jj,1:nlay_i,jl) )
[5134]	363
	364	zdv = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl)
	365
	366	IF ( ( zdv > 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. &
	367	& ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN
	368
	369	rswitch = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) )
	370	a_i(ji,jj,jl) = rswitch * ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / MAX( zhimax(ji,jj,jl), epsi20 )
	371
	372	! small correction due to *rswitch for a_i
	373	v_i (ji,jj,jl) = rswitch * v_i (ji,jj,jl)
	374	v_s (ji,jj,jl) = rswitch * v_s (ji,jj,jl)
	375	smv_i(ji,jj,jl) = rswitch * smv_i(ji,jj,jl)
	376	e_s(ji,jj,1,jl) = rswitch * e_s(ji,jj,1,jl)
	377	e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl)
	378
	379	! Update mass fluxes
	380	wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice
	381	wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
	382	sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice
	383	hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
	384	hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * r1_rdtice ! W.m-2 <0
	385
[4161]	386	ENDIF
	387
	388	ENDIF
[5123]	389
[825]	390	END DO
	391	END DO
	392	END DO
[4688]	393	! -------------------------------------------------
[5123]	394
	395	!--------------------------------------
	396	! Impose a_i < amax in mono-category
	397	!--------------------------------------
	398	!
	399	IF ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) THEN ! simple conservative piling, comparable with LIM2
	400	DO jj = 1, jpj
	401	DO ji = 1, jpi
	402	a_i(ji,jj,1) = MIN( a_i(ji,jj,1), rn_amax )
	403	END DO
	404	END DO
	405	ENDIF
[825]	406
[4688]	407	! --- diags ---
[825]	408	DO jj = 1, jpj
	409	DO ji = 1, jpi
[5123]	410	diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) * r1_rdtice
	411	diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) * r1_rdtice
[4990]	412
[5123]	413	diag_trp_vi (ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice
	414	diag_trp_vs (ji,jj) = SUM( v_s(ji,jj,:) - zvsold(ji,jj,:) ) * r1_rdtice
	415	diag_trp_smv(ji,jj) = SUM( smv_i(ji,jj,:) - zsmvold(ji,jj,:) ) * r1_rdtice
[825]	416	END DO
	417	END DO
	418
[4990]	419	! --- agglomerate variables -----------------
[4688]	420	vt_i (:,:) = 0._wp
	421	vt_s (:,:) = 0._wp
	422	at_i (:,:) = 0._wp
[825]	423	DO jl = 1, jpl
	424	DO jj = 1, jpj
	425	DO ji = 1, jpi
[5134]	426	vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl)
	427	vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl)
	428	at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl)
[4688]	429	END DO
	430	END DO
	431	END DO
[825]	432
[5134]	433	! --- open water = 1 if at_i=0 --------------------------------
[4161]	434	DO jj = 1, jpj
	435	DO ji = 1, jpi
[4990]	436	rswitch = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
[5123]	437	ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
[4161]	438	END DO
[4688]	439	END DO
[4161]	440
[4688]	441	! conservation test
	442	IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[4161]	443
[825]	444	ENDIF
	445
[5123]	446	CALL lim_var_glo2eqv ! equivalent variables, requested for rafting
[4161]	447
[5123]	448	! -------------------------------------------------
	449	! control prints
	450	! -------------------------------------------------
[5128]	451	IF( ln_icectl ) CALL lim_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' )
[2715]	452	!
[5134]	453	CALL wrk_dealloc( jpi,jpj, zsm, zatold, zeiold, zesold )
	454	CALL wrk_dealloc( jpi,jpj,jpl, z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
	455	CALL wrk_dealloc( jpi,jpj,1, z0opw )
	456	CALL wrk_dealloc( jpi,jpj,nlay_i+1,jpl, z0ei )
[5123]	457	CALL wrk_dealloc( jpi,jpj,jpl, zviold, zvsold, zhimax, zsmvold )
	458	!
[4161]	459	IF( nn_timing == 1 ) CALL timing_stop('limtrp')
[5123]	460
[825]	461	END SUBROUTINE lim_trp
	462
	463	#else
	464	!!----------------------------------------------------------------------
	465	!! Default option Empty Module No sea-ice model
	466	!!----------------------------------------------------------------------
	467	CONTAINS
	468	SUBROUTINE lim_trp ! Empty routine
	469	END SUBROUTINE lim_trp
	470	#endif
	471	!!======================================================================
	472	END MODULE limtrp

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

Download in other formats: