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

Last change on this file since 3594 was 3558, checked in by rblod, 12 years ago
Fix issues when using key_nosignedzeo, see ticket #996
Property svn:keywords set to `Id`
File size: 26.2 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	!!----------------------------------------------------------------------
[2715]	16	USE phycst ! physical constant
	17	USE dom_oce ! ocean domain
	18	USE sbc_oce ! ocean surface boundary condition
	19	USE par_ice ! LIM-3 parameter
	20	USE dom_ice ! LIM-3 domain
	21	USE ice ! LIM-3 variables
	22	USE limadv ! LIM-3 advection
	23	USE limhdf ! LIM-3 horizontal diffusion
[825]	24	USE in_out_manager ! I/O manager
[2715]	25	USE lbclnk ! lateral boundary conditions -- MPP exchanges
	26	USE lib_mpp ! MPP library
[3294]	27	USE wrk_nemo ! work arrays
[863]	28	USE prtctl ! Print control
[3558]	29	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[825]	30
	31	IMPLICIT NONE
	32	PRIVATE
	33
[2715]	34	PUBLIC lim_trp ! called by ice_step
[825]	35
[2715]	36	REAL(wp) :: epsi06 = 1.e-06_wp ! constant values
	37	REAL(wp) :: epsi03 = 1.e-03_wp
	38	REAL(wp) :: zeps10 = 1.e-10_wp
	39	REAL(wp) :: epsi16 = 1.e-16_wp
	40	REAL(wp) :: rzero = 0._wp
	41	REAL(wp) :: rone = 1._wp
[825]	42
[2777]	43	REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:) :: zs0e
	44
[825]	45	!! * Substitution
	46	# include "vectopt_loop_substitute.h90"
	47	!!----------------------------------------------------------------------
[2715]	48	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
[1156]	49	!! $Id$
[2715]	50	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	51	!!----------------------------------------------------------------------
	52	CONTAINS
	53
[921]	54	SUBROUTINE lim_trp( kt )
[825]	55	!!-------------------------------------------------------------------
	56	!! * ROUTINE lim_trp *
	57	!!
	58	!! ** purpose : advection/diffusion process of sea ice
	59	!!
	60	!! ** method : variables included in the process are scalar,
	61	!! other values are considered as second order.
	62	!! For advection, a second order Prather scheme is used.
	63	!!
	64	!! ** action :
	65	!!---------------------------------------------------------------------
[2715]	66	INTEGER, INTENT(in) :: kt ! number of iteration
	67	!
	68	INTEGER :: ji, jj, jk, jl, layer ! dummy loop indices
	69	INTEGER :: initad ! number of sub-timestep for the advection
[2777]	70	INTEGER :: ierr ! error status
[2715]	71	REAL(wp) :: zindb , zindsn , zindic ! local scalar
	72	REAL(wp) :: zusvosn, zusvoic, zbigval ! - -
	73	REAL(wp) :: zcfl , zusnit , zrtt ! - -
	74	REAL(wp) :: ze , zsal , zage ! - -
	75	!
[3294]	76	REAL(wp), POINTER, DIMENSION(:,:) :: zui_u, zvi_v, zsm, zs0at, zs0ow
	77	REAL(wp), POINTER, DIMENSION(:,:,:) :: zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi
	78	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zs0e
[2715]	79	!!---------------------------------------------------------------------
[825]	80
[3294]	81	CALL wrk_alloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow )
	82	CALL wrk_alloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi )
	83	CALL wrk_alloc( jpi, jpj, jkmax, jpl, zs0e )
[825]	84
[2715]	85	IF( numit == nstart .AND. lwp ) THEN
	86	WRITE(numout,*)
	87	IF( ln_limdyn ) THEN ; WRITE(numout,*) 'lim_trp : Ice transport '
	88	ELSE ; WRITE(numout,*) 'lim_trp : No ice advection as ln_limdyn = ', ln_limdyn
	89	ENDIF
	90	WRITE(numout,*) '~~~~~~~~~~~~'
	91	ENDIF
	92
[825]	93	zsm(:,:) = area(:,:)
	94
[2715]	95	! !-------------------------------------!
	96	IF( ln_limdyn ) THEN ! Advection of sea ice properties !
	97	! !-------------------------------------!
	98	!
[825]	99
	100	!-------------------------
[2715]	101	! transported fields
[825]	102	!-------------------------
[2715]	103	! Snow vol, ice vol, salt and age contents, area
	104	zs0ow(:,:) = ato_i(:,:) * area(:,:) ! Open water area
	105	DO jl = 1, jpl
	106	zs0sn (:,:,jl) = v_s (:,:,jl) * area(:,:) ! Snow volume
	107	zs0ice(:,:,jl) = v_i (:,:,jl) * area(:,:) ! Ice volume
	108	zs0a (:,:,jl) = a_i (:,:,jl) * area(:,:) ! Ice area
	109	zs0sm (:,:,jl) = smv_i(:,:,jl) * area(:,:) ! Salt content
	110	zs0oi (:,:,jl) = oa_i (:,:,jl) * area(:,:) ! Age content
	111	zs0c0 (:,:,jl) = e_s (:,:,1,jl) ! Snow heat content
	112	zs0e (:,:,:,jl) = e_i (:,:,:,jl) ! Ice heat content
[825]	113	END DO
	114
[2715]	115	!--------------------------
	116	! Advection of Ice fields (Prather scheme)
	117	!--------------------------
[825]	118	! If ice drift field is too fast, use an appropriate time step for advection.
[2715]	119	! CFL test for stability
	120	zcfl = MAXVAL( ABS( u_ice(:,:) ) * rdt_ice / e1u(:,:) )
	121	zcfl = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice / e2v(:,:) ) )
	122	IF(lk_mpp ) CALL mpp_max( zcfl )
	123	!!gm more readability:
	124	! IF( zcfl > 0.5 ) THEN ; initad = 2 ; zusnit = 0.5_wp
	125	! ELSE ; initad = 1 ; zusnit = 1.0_wp
	126	! ENDIF
	127	!!gm end
[825]	128	initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )
	129	zusnit = 1.0 / REAL( initad )
[2715]	130	IF( zcfl > 0.5 .AND. lwp ) &
	131	WRITE(numout,*) 'lim_trp_2 : CFL violation at day ', nday, ', cfl = ', zcfl, &
	132	& ': the ice time stepping is split in two'
[921]	133
[2715]	134	IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==!
[825]	135	DO jk = 1,initad
[2715]	136	CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area
	137	& sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) )
	138	CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:), &
	139	& sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) )
[825]	140	DO jl = 1, jpl
[2715]	141	CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume ---
	142	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
	143	CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl), &
	144	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
	145	CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---
	146	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
	147	CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), &
	148	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
	149	CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---
	150	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
	151	CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), &
	152	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
	153	CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---
	154	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
	155	CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl), &
	156	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
	157	CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---
	158	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
	159	CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0a (:,:,jl), sxa (:,:,jl), &
	160	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
	161	CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents ---
	162	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
	163	CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), &
	164	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
	165	DO layer = 1, nlay_i !--- ice heat contents ---
	166	CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &
	167	& sxxe(:,:,layer,jl), sye (:,:,layer,jl), &
	168	& syye(:,:,layer,jl), sxye(:,:,layer,jl) )
	169	CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &
	170	& sxxe(:,:,layer,jl), sye (:,:,layer,jl), &
	171	& syye(:,:,layer,jl), sxye(:,:,layer,jl) )
[825]	172	END DO
	173	END DO
	174	END DO
	175	ELSE
	176	DO jk = 1, initad
[3554]	177	CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0ow (:,:), sxopw(:,:), & !--- ice open water area
[2715]	178	& sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) )
[3554]	179	CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:), &
[2715]	180	& sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:) )
[825]	181	DO jl = 1, jpl
[3554]	182	CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume ---
[2715]	183	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
[3554]	184	CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl), &
[2715]	185	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
[3554]	186	CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---
[2715]	187	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
[3554]	188	CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), &
[2715]	189	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
[3554]	190	CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---
[2715]	191	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
[3554]	192	CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), &
[2715]	193	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
	194
[3554]	195	CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---
[2715]	196	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
[3554]	197	CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl), &
[2715]	198	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
[3554]	199	CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---
[2715]	200	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
[3554]	201	CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0a (:,:,jl), sxa (:,:,jl), &
[2715]	202	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
[3554]	203	CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents ---
[2715]	204	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
[3554]	205	CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), &
[2715]	206	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
	207	DO layer = 1, nlay_i !--- ice heat contents ---
[3554]	208	CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &
[2715]	209	& sxxe(:,:,layer,jl), sye (:,:,layer,jl), &
	210	& syye(:,:,layer,jl), sxye(:,:,layer,jl) )
[3554]	211	CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl), &
[2715]	212	& sxxe(:,:,layer,jl), sye (:,:,layer,jl), &
	213	& syye(:,:,layer,jl), sxye(:,:,layer,jl) )
[825]	214	END DO
	215	END DO
	216	END DO
	217	ENDIF
	218
	219	!-------------------------------------------
	220	! Recover the properties from their contents
	221	!-------------------------------------------
[2715]	222	zs0ow(:,:) = zs0ow(:,:) / area(:,:)
[825]	223	DO jl = 1, jpl
	224	zs0ice(:,:,jl) = zs0ice(:,:,jl) / area(:,:)
	225	zs0sn (:,:,jl) = zs0sn (:,:,jl) / area(:,:)
	226	zs0sm (:,:,jl) = zs0sm (:,:,jl) / area(:,:)
	227	zs0oi (:,:,jl) = zs0oi (:,:,jl) / area(:,:)
	228	zs0a (:,:,jl) = zs0a (:,:,jl) / area(:,:)
	229	zs0c0 (:,:,jl) = zs0c0 (:,:,jl) / area(:,:)
	230	DO jk = 1, nlay_i
	231	zs0e(:,:,jk,jl) = zs0e(:,:,jk,jl) / area(:,:)
	232	END DO
	233	END DO
	234
[921]	235	!------------------------------------------------------------------------------!
	236	! 4) Diffusion of Ice fields
	237	!------------------------------------------------------------------------------!
[825]	238
[2715]	239	!--------------------------------
	240	! diffusion of open water area
	241	!--------------------------------
	242	zs0at(:,:) = zs0a(:,:,1) ! total ice fraction
	243	DO jl = 2, jpl
	244	zs0at(:,:) = zs0at(:,:) + zs0a(:,:,jl)
	245	END DO
	246	!
	247	! ! Masked eddy diffusivity coefficient at ocean U- and V-points
	248	DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row
	249	DO ji = 1 , fs_jpim1 ! vector opt.
	250	pahu(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji ,jj) ) ) ) &
	251	& * ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji+1,jj) ) ) ) * ahiu(ji,jj)
	252	pahv(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji,jj ) ) ) ) &
	253	& * ( 1._wp - MAX( rzero, SIGN( rone,- zs0at(ji,jj+1) ) ) ) * ahiv(ji,jj)
	254	END DO
	255	END DO
	256	!
	257	CALL lim_hdf( zs0ow (:,:) ) ! Diffusion
	258
[921]	259	!------------------------------------
[2715]	260	! Diffusion of other ice variables
[921]	261	!------------------------------------
[825]	262	DO jl = 1, jpl
[2715]	263	! ! Masked eddy diffusivity coefficient at ocean U- and V-points
	264	DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row
	265	DO ji = 1 , fs_jpim1 ! vector opt.
	266	pahu(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji ,jj,jl) ) ) ) &
	267	& * ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)
	268	pahv(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji,jj ,jl) ) ) ) &
	269	& * ( 1._wp - MAX( rzero, SIGN( rone,- zs0a(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)
[825]	270	END DO
[921]	271	END DO
[825]	272
	273	CALL lim_hdf( zs0ice (:,:,jl) )
	274	CALL lim_hdf( zs0sn (:,:,jl) )
	275	CALL lim_hdf( zs0sm (:,:,jl) )
	276	CALL lim_hdf( zs0oi (:,:,jl) )
	277	CALL lim_hdf( zs0a (:,:,jl) )
	278	CALL lim_hdf( zs0c0 (:,:,jl) )
	279	DO jk = 1, nlay_i
	280	CALL lim_hdf( zs0e (:,:,jk,jl) )
[2715]	281	END DO
	282	END DO
[825]	283
[921]	284	!-----------------------------------------
[2715]	285	! Remultiply everything by ice area
[921]	286	!-----------------------------------------
[2715]	287	zs0ow(:,:) = MAX( rzero, zs0ow(:,:) * area(:,:) )
[825]	288	DO jl = 1, jpl
	289	zs0ice(:,:,jl) = MAX( rzero, zs0ice(:,:,jl) * area(:,:) ) !!bug: est-ce utile
	290	zs0sn (:,:,jl) = MAX( rzero, zs0sn (:,:,jl) * area(:,:) ) !!bug: cf /area juste apres
	291	zs0sm (:,:,jl) = MAX( rzero, zs0sm (:,:,jl) * area(:,:) ) !!bug: cf /area juste apres
	292	zs0oi (:,:,jl) = MAX( rzero, zs0oi (:,:,jl) * area(:,:) )
	293	zs0a (:,:,jl) = MAX( rzero, zs0a (:,:,jl) * area(:,:) ) !! suppress both change le resultat
	294	zs0c0 (:,:,jl) = MAX( rzero, zs0c0 (:,:,jl) * area(:,:) )
	295	DO jk = 1, nlay_i
	296	zs0e(:,:,jk,jl) = MAX( rzero, zs0e (:,:,jk,jl) * area(:,:) )
	297	END DO ! jk
	298	END DO ! jl
	299
[921]	300	!------------------------------------------------------------------------------!
	301	! 5) Update and limit ice properties after transport
	302	!------------------------------------------------------------------------------!
[825]	303
[921]	304	!--------------------------------------------------
	305	! 5.1) Recover mean values over the grid squares.
	306	!--------------------------------------------------
[825]	307
	308	DO jl = 1, jpl
	309	DO jk = 1, nlay_i
	310	DO jj = 1, jpj
	311	DO ji = 1, jpi
[2715]	312	zs0e(ji,jj,jk,jl) = MAX( rzero, zs0e(ji,jj,jk,jl) / area(ji,jj) )
[825]	313	END DO
	314	END DO
	315	END DO
	316	END DO
	317
	318	DO jj = 1, jpj
	319	DO ji = 1, jpi
[2715]	320	zs0ow(ji,jj) = MAX( rzero, zs0ow (ji,jj) / area(ji,jj) )
[825]	321	END DO
	322	END DO
[921]	323
[2715]	324	zs0at(:,:) = 0._wp
[825]	325	DO jl = 1, jpl
	326	DO jj = 1, jpj
	327	DO ji = 1, jpi
	328	zs0sn (ji,jj,jl) = MAX( rzero, zs0sn (ji,jj,jl)/area(ji,jj) )
	329	zs0ice(ji,jj,jl) = MAX( rzero, zs0ice(ji,jj,jl)/area(ji,jj) )
	330	zs0sm (ji,jj,jl) = MAX( rzero, zs0sm (ji,jj,jl)/area(ji,jj) )
	331	zs0oi (ji,jj,jl) = MAX( rzero, zs0oi (ji,jj,jl)/area(ji,jj) )
	332	zs0a (ji,jj,jl) = MAX( rzero, zs0a (ji,jj,jl)/area(ji,jj) )
	333	zs0c0 (ji,jj,jl) = MAX( rzero, zs0c0 (ji,jj,jl)/area(ji,jj) )
	334	zs0at (ji,jj) = zs0at(ji,jj) + zs0a(ji,jj,jl)
	335	END DO
	336	END DO
	337	END DO
	338
[921]	339	!---------------------------------------------------------
	340	! 5.2) Snow thickness, Ice thickness, Ice concentrations
	341	!---------------------------------------------------------
[825]	342	DO jj = 1, jpj
	343	DO ji = 1, jpi
[2715]	344	zindb = MAX( 0._wp , SIGN( 1.0, zs0at(ji,jj) - zeps10) )
	345	zs0ow(ji,jj) = ( 1._wp - zindb ) + zindb * MAX( zs0ow(ji,jj), 0._wp )
	346	ato_i(ji,jj) = zs0ow(ji,jj)
[825]	347	END DO
	348	END DO
	349
[2715]	350	DO jl = 1, jpl ! Remove very small areas
[825]	351	DO jj = 1, jpj
	352	DO ji = 1, jpi
	353	zindb = MAX( 0.0 , SIGN( 1.0, zs0a(ji,jj,jl) - zeps10) )
[2715]	354	!
	355	zs0a(ji,jj,jl) = zindb * MIN( zs0a(ji,jj,jl), 0.99 )
	356	v_s(ji,jj,jl) = zindb * zs0sn (ji,jj,jl)
	357	v_i(ji,jj,jl) = zindb * zs0ice(ji,jj,jl)
	358	!
	359	zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - zeps10 ) )
	360	zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) )
	361	zindb = MAX( zindsn, zindic )
	362	zs0a(ji,jj,jl) = zindb * zs0a(ji,jj,jl) !ice concentration
	363	a_i (ji,jj,jl) = zs0a(ji,jj,jl)
	364	v_s (ji,jj,jl) = zindsn * v_s(ji,jj,jl)
	365	v_i (ji,jj,jl) = zindic * v_i(ji,jj,jl)
[825]	366	END DO
	367	END DO
	368	END DO
	369
	370	DO jj = 1, jpj
	371	DO ji = 1, jpi
[2715]	372	zs0at(ji,jj) = SUM( zs0a(ji,jj,1:jpl) )
[825]	373	END DO
	374	END DO
	375
[921]	376	!----------------------
	377	! 5.3) Ice properties
	378	!----------------------
[825]	379
[2715]	380	zbigval = 1.d+13
[825]	381
	382	DO jl = 1, jpl
	383	DO jj = 1, jpj
	384	DO ji = 1, jpi
	385
	386	! Switches and dummy variables
	387	zusvosn = 1.0/MAX( v_s(ji,jj,jl) , epsi16 )
	388	zusvoic = 1.0/MAX( v_i(ji,jj,jl) , epsi16 )
	389	zrtt = 173.15 * rone
	390	zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - zeps10 ) )
	391	zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) )
	392	zindb = MAX( zindsn, zindic )
	393
	394	! Ice salinity and age
[2715]	395	zsal = MAX( MIN( (rhoic-rhosn)/rhoic*sss_m(ji,jj) , &
	396	zusvoic * zs0sm(ji,jj,jl) ), s_i_min ) * v_i(ji,jj,jl)
[825]	397	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &
	398	smv_i(ji,jj,jl) = zindiczsal + (1.0-zindic)0.0
	399
[2715]	400	zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / &
	401	MAX( a_i(ji,jj,jl), epsi16 ) ), 0.0 ) * a_i(ji,jj,jl)
[825]	402	oa_i (ji,jj,jl) = zindic*zage
	403
	404	! Snow heat content
	405	ze = MIN( MAX( 0.0, zs0c0(ji,jj,jl)*area(ji,jj) ), zbigval )
	406	e_s(ji,jj,1,jl) = zindsn * ze + (1.0 - zindsn) * 0.0
	407
	408	END DO !ji
	409	END DO !jj
	410	END DO ! jl
	411
	412	DO jl = 1, jpl
	413	DO jk = 1, nlay_i
	414	DO jj = 1, jpj
	415	DO ji = 1, jpi
	416	! Ice heat content
	417	zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) )
	418	ze = MIN( MAX( 0.0, zs0e(ji,jj,jk,jl)*area(ji,jj) ), zbigval )
	419	e_i(ji,jj,jk,jl) = zindic * ze + ( 1.0 - zindic ) * 0.0
	420	END DO !ji
	421	END DO ! jj
	422	END DO ! jk
	423	END DO ! jl
	424
	425	ENDIF
	426
[863]	427	IF(ln_ctl) THEN ! Control print
[867]	428	CALL prt_ctl_info(' ')
	429	CALL prt_ctl_info(' - Cell values : ')
	430	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
[863]	431	CALL prt_ctl(tab2d_1=area , clinfo1=' lim_trp : cell area :')
	432	CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_trp : at_i :')
	433	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_trp : vt_i :')
	434	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_trp : vt_s :')
	435	DO jl = 1, jpl
[867]	436	CALL prt_ctl_info(' ')
[863]	437	CALL prt_ctl_info(' - Category : ', ivar1=jl)
	438	CALL prt_ctl_info(' ~~~~~~~~~~')
	439	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_trp : a_i : ')
	440	CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_trp : ht_i : ')
	441	CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_trp : ht_s : ')
	442	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_trp : v_i : ')
	443	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_trp : v_s : ')
	444	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_trp : e_s : ')
	445	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_trp : t_su : ')
	446	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_trp : t_snow : ')
	447	CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_trp : sm_i : ')
	448	CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_trp : smv_i : ')
	449	DO jk = 1, nlay_i
[867]	450	CALL prt_ctl_info(' ')
[863]	451	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
	452	CALL prt_ctl_info(' ~~~~~~~')
	453	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_trp : t_i : ')
	454	CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_trp : e_i : ')
	455	END DO
	456	END DO
	457	ENDIF
[2715]	458	!
[3294]	459	CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow )
	460	CALL wrk_dealloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi )
	461	CALL wrk_dealloc( jpi, jpj, jkmax, jpl, zs0e )
[2715]	462	!
[825]	463	END SUBROUTINE lim_trp
	464
	465	#else
	466	!!----------------------------------------------------------------------
	467	!! Default option Empty Module No sea-ice model
	468	!!----------------------------------------------------------------------
	469	CONTAINS
	470	SUBROUTINE lim_trp ! Empty routine
	471	END SUBROUTINE lim_trp
	472	#endif
	473
	474	!!======================================================================
	475	END MODULE limtrp

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

Download in other formats: