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

Context Navigation

source: trunk/NEMO/LIM_SRC_3/limtrp.F90 @ 973

Last change on this file since 973 was 921, checked in by rblod, 16 years ago
Correct indentation and print for debug in LIM3, see ticket #134, step I
File size: 28.1 KB

Rev	Line
[825]	1	MODULE limtrp
	2	!!======================================================================
	3	!! * MODULE limtrp *
	4	!! LIM transport ice model : sea-ice advection/diffusion
	5	!!======================================================================
	6	#if defined key_lim3
	7	!!----------------------------------------------------------------------
[834]	8	!! 'key_lim3' LIM3 sea-ice model
[825]	9	!!----------------------------------------------------------------------
	10	!! lim_trp : advection/diffusion process of sea ice
	11	!! lim_trp_init : initialization and namelist read
	12	!!----------------------------------------------------------------------
	13	!! * Modules used
	14	USE phycst
	15	USE dom_oce
	16	USE daymod
	17	USE in_out_manager ! I/O manager
	18	USE ice_oce ! ice variables
[888]	19	USE sbc_oce ! Surface boundary condition: ocean fields
[825]	20	USE dom_ice
	21	USE ice
	22	USE iceini
	23	USE limistate
	24	USE limadv
	25	USE limhdf
	26	USE lbclnk
	27	USE lib_mpp
	28	USE par_ice
[863]	29	USE prtctl ! Print control
[825]	30
	31	IMPLICIT NONE
	32	PRIVATE
	33
	34	!! * Routine accessibility
	35	PUBLIC lim_trp ! called by ice_step
	36
	37	!! * Shared module variables
	38	REAL(wp), PUBLIC :: & !:
	39	bound = 0.e0 !: boundary condit. (0.0 no-slip, 1.0 free-slip)
	40
	41	!! * Module variables
	42	REAL(wp) :: & ! constant values
	43	epsi06 = 1.e-06 , &
	44	epsi03 = 1.e-03 , &
	45	epsi16 = 1.e-16 , &
	46	rzero = 0.e0 , &
	47	rone = 1.e0 , &
	48	zeps10 = 1.e-10
	49
	50	!! * Substitution
	51	# include "vectopt_loop_substitute.h90"
	52	!!----------------------------------------------------------------------
[834]	53	!! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (2008)
[888]	54	!! $ Id: $
	55	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[825]	56	!!----------------------------------------------------------------------
	57	CONTAINS
	58
[921]	59	SUBROUTINE lim_trp( kt )
[825]	60	!!-------------------------------------------------------------------
	61	!! * ROUTINE lim_trp *
	62	!!
	63	!! ** purpose : advection/diffusion process of sea ice
	64	!!
	65	!! ** method : variables included in the process are scalar,
	66	!! other values are considered as second order.
	67	!! For advection, a second order Prather scheme is used.
	68	!!
	69	!! ** action :
	70	!!
	71	!! History :
	72	!! 1.0 ! 00-01 (M.A. Morales Maqueda, H. Goosse, and T. Fichefet) Original code
	73	!! ! 01-05 (G. Madec, R. Hordoir) opa norm
	74	!! 2.0 ! 04-01 (G. Madec, C. Ethe) F90, mpp
	75	!! 3.0 ! 05-11 (M. Vancoppenolle) Multi-layer sea ice, salinity variations
	76	!!---------------------------------------------------------------------
[921]	77	INTEGER, INTENT(in) :: kt ! number of iteration
[825]	78	!! * Local Variables
[834]	79	INTEGER :: ji, jj, jk, jl, layer, & ! dummy loop indices
[921]	80	initad ! number of sub-timestep for the advection
[825]	81	INTEGER :: ji_maxu, ji_maxv, jj_maxu, jj_maxv
	82
	83	REAL(wp) :: &
	84	zindb , &
[834]	85	zindsn , &
	86	zindic , &
	87	zusvosn, &
	88	zusvoic, &
	89	zvbord , &
	90	zcfl , &
	91	zusnit , &
	92	zrtt, zsal, zage, &
	93	zbigval, ze, &
[825]	94	zmaxu, zmaxv
	95
	96	REAL(wp), DIMENSION(jpi,jpj) :: & ! temporary workspace
	97	zui_u , zvi_v , zsm , &
	98	zs0at, zs0ow
	99
	100	REAL(wp), DIMENSION(jpi,jpj,jpl):: & ! temporary workspace
	101	zs0ice, zs0sn, zs0a , &
[834]	102	zs0c0 , &
[825]	103	zs0sm , zs0oi
	104
[921]	105	! MHE Multilayer heat content
[825]	106	REAL(wp), DIMENSION(jpi,jpj,jkmax,jpl) :: & ! temporary workspace
	107	zs0e
	108
	109	!---------------------------------------------------------------------
	110
	111	IF( numit == nstart ) CALL lim_trp_init ! Initialization (first time-step only)
	112
	113	zsm(:,:) = area(:,:)
	114
[921]	115	IF( ln_limdyn .AND. lwp ) THEN
	116	IF( kt == nit000 ) THEN
	117	WRITE(numout,*) ' lim_trp : Ice Advection'
	118	WRITE(numout,*) ' ~~~~~~~'
	119	ENDIF
[825]	120
[921]	121	!-----------------------------------------------------------------------------!
	122	! 1) CFL Test
	123	!-----------------------------------------------------------------------------!
	124
[825]	125	!------------------------------------------
	126	! ice velocities at ocean U- and V-points
	127	!------------------------------------------
[921]	128
[825]	129	! zvbord factor between 1 and 2 to take into account slip or no-slip boundary conditions.
	130	zvbord = 1.0 + ( 1.0 - bound )
	131	DO jj = 1, jpjm1
[868]	132	DO ji = 1, fs_jpim1
[825]	133	zui_u(ji,jj) = u_ice(ji,jj)
	134	zvi_v(ji,jj) = v_ice(ji,jj)
	135	END DO
	136	END DO
	137
	138	! Lateral boundary conditions
	139	CALL lbc_lnk( zui_u, 'U', -1. )
	140	CALL lbc_lnk( zvi_v, 'V', -1. )
	141
	142	!-------------------------
	143	! CFL test for stability
	144	!-------------------------
	145
	146	zcfl = 0.e0
	147	zcfl = MAX( zcfl, MAXVAL( ABS( zui_u(1:jpim1, : ) ) * rdt_ice / e1u(1:jpim1, : ) ) )
	148	zcfl = MAX( zcfl, MAXVAL( ABS( zvi_v( : ,1:jpjm1) ) * rdt_ice / e2v( : ,1:jpjm1) ) )
	149
	150	zmaxu = 0.0
	151	zmaxv = 0.0
	152	DO ji = 1, jpim1
	153	DO jj = 1, jpjm1
	154	IF ( (ABS(zui_u(ji,jj)) .GT. zmaxu) ) THEN
	155	zmaxu = MAX(zui_u(ji,jj), zmaxu )
	156	ji_maxu = ji
	157	jj_maxu = jj
	158	ENDIF
	159	IF ( (ABS(zvi_v(ji,jj)) .GT. zmaxv) ) THEN
	160	zmaxv = MAX(zvi_v(ji,jj), zmaxv )
	161	ji_maxv = ji
	162	jj_maxv = jj
	163	ENDIF
	164	END DO
	165	END DO
	166
	167	IF (lk_mpp ) CALL mpp_max(zcfl)
	168
	169	IF ( zcfl > 0.5 .AND. lwp ) &
[921]	170	WRITE(numout,*) 'lim_trp : violation of cfl criterion the ',nday,'th day, cfl = ',zcfl
[825]	171
[921]	172	!-----------------------------------------------------------------------------!
	173	! 2) Computation of transported fields
	174	!-----------------------------------------------------------------------------!
[825]	175
	176	!------------------------------------------------------
	177	! 1.1) Snow vol, ice vol, salt and age contents, area
	178	!------------------------------------------------------
	179
	180	zs0ow (:,:) = ato_i(:,:) * area(:,:) ! Open water area
	181	DO jl = 1, jpl !sum over thickness categories
	182	! area -> is the unmasked and masked area of T-grid cell
	183	zs0sn (:,:,jl) = v_s(:,:,jl) * area(:,:) ! Snow volume.
	184	zs0ice(:,:,jl) = v_i(:,:,jl) * area(:,:) ! Ice volume.
	185	zs0a (:,:,jl) = a_i(:,:,jl) * area(:,:) ! Ice area
	186	zs0sm (:,:,jl) = smv_i(:,:,jl) * area(:,:) ! Salt content
	187	zs0oi (:,:,jl) = oa_i (:,:,jl) * area(:,:) ! Age content
	188
[921]	189	!----------------------------------
	190	! 1.2) Ice and snow heat contents
	191	!----------------------------------
[825]	192
	193	zs0c0 (:,:,jl) = e_s(:,:,1,jl) ! Snow heat cont.
	194	DO jk = 1, nlay_i
	195	zs0e(:,:,jk,jl) = e_i(:,:,jk,jl) ! Ice heat content
	196	END DO
	197	END DO
	198
[921]	199	!-----------------------------------------------------------------------------!
	200	! 3) Advection of Ice fields
	201	!-----------------------------------------------------------------------------!
[825]	202
	203	! If ice drift field is too fast, use an appropriate time step for advection.
	204	initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )
	205	zusnit = 1.0 / REAL( initad )
[921]	206
[825]	207	IF ( MOD( nday , 2 ) == 0) THEN
	208	DO jk = 1,initad
	209	!--- ice open water area
	210	CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0ow(:,:) , sxopw(:,:) , &
[921]	211	sxxopw(:,:), syopw(:,:) , &
	212	syyopw(:,:), sxopw(:,:) )
[825]	213	CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0ow(:,:) , sxopw (:,:) , &
[921]	214	sxxopw(:,:), syopw (:,:) , &
	215	syyopw(:,:), sxyopw(:,:) )
[825]	216	DO jl = 1, jpl
	217	!--- ice volume ---
	218	CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0ice(:,:,jl) , sxice (:,:,jl) , &
[921]	219	sxxice(:,:,jl) , syice (:,:,jl) , &
	220	syyice(:,:,jl) , sxyice(:,:,jl) )
[825]	221	CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0ice(:,:,jl) , sxice (:,:,jl) , &
[921]	222	sxxice(:,:,jl) , syice (:,:,jl) , &
	223	syyice(:,:,jl) , sxyice(:,:,jl) )
[825]	224	!--- snow volume ---
	225	CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0sn (:,:,jl) , sxsn (:,:,jl) , &
[921]	226	sxxsn (:,:,jl) , sysn (:,:,jl) , &
	227	syysn (:,:,jl) , sxysn (:,:,jl) )
[825]	228	CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0sn (:,:,jl) , sxsn (:,:,jl) , &
[921]	229	sxxsn (:,:,jl) , sysn (:,:,jl) , &
	230	syysn (:,:,jl) , sxysn (:,:,jl) )
[825]	231	!--- ice salinity ---
	232	CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0sm (:,:,jl) , sxsal (:,:,jl) , &
[921]	233	sxxsal(:,:,jl) , sysal (:,:,jl) , &
	234	syysal(:,:,jl) , sxysal(:,:,jl) )
[825]	235	CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0sm (:,:,jl) , sxsal (:,:,jl) , &
[921]	236	sxxsal(:,:,jl) , sysal (:,:,jl) , &
	237	syysal(:,:,jl) , sxysal(:,:,jl) )
[825]	238	!--- ice age ---
	239	CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0oi (:,:,jl) , sxage (:,:,jl) , &
[921]	240	sxxage(:,:,jl) , syage (:,:,jl) , &
	241	syyage(:,:,jl) , sxyage(:,:,jl) )
[825]	242	CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0oi (:,:,jl) , sxage (:,:,jl) , &
[921]	243	sxxage(:,:,jl) , syage (:,:,jl) , &
	244	syyage(:,:,jl) , sxyage(:,:,jl) )
[825]	245	!--- ice concentrations ---
	246	CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0a (:,:,jl) , sxa (:,:,jl) , &
[921]	247	sxxa (:,:,jl) , sya (:,:,jl) , &
	248	syya (:,:,jl) , sxya (:,:,jl) )
[825]	249	CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0a (:,:,jl) , sxa (:,:,jl) , &
[921]	250	sxxa (:,:,jl) , sya (:,:,jl) , &
	251	syya (:,:,jl) , sxya (:,:,jl) )
[825]	252	!--- ice / snow thermal energetic contents ---
	253	CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0c0 (:,:,jl) , sxc0 (:,:,jl) , &
[921]	254	sxxc0 (:,:,jl) , syc0 (:,:,jl) , &
	255	syyc0 (:,:,jl) , sxyc0 (:,:,jl) )
[825]	256	CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0c0 (:,:,jl) , sxc0 (:,:,jl) , &
[921]	257	sxxc0 (:,:,jl) , syc0 (:,:,jl) , &
	258	syyc0 (:,:,jl) , sxyc0 (:,:,jl) )
[825]	259	DO layer = 1, nlay_i
	260	CALL lim_adv_x( zusnit, zui_u, rone , zsm, &
[921]	261	zs0e(:,:,layer,jl) , sxe (:,:,layer,jl) , &
	262	sxxe(:,:,layer,jl) , sye (:,:,layer,jl) , &
	263	syye(:,:,layer,jl) , sxye(:,:,layer,jl) )
[825]	264	CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, &
[921]	265	zs0e(:,:,layer,jl) , sxe (:,:,layer,jl) , &
	266	sxxe(:,:,layer,jl) , sye (:,:,layer,jl) , &
	267	syye(:,:,layer,jl) , sxye(:,:,layer,jl) )
[825]	268	END DO
	269	END DO
	270	END DO
	271	ELSE
	272	DO jk = 1, initad
	273	!--- ice volume ---
	274	CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0ow (:,:) , sxopw (:,:) , &
[921]	275	sxxopw(:,:) , syopw (:,:) , &
	276	syyopw(:,:) , sxyopw(:,:) )
[825]	277	CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0ow (:,:) , sxopw (:,:) , &
[921]	278	sxxopw(:,:) , syopw (:,:) , &
	279	syyopw(:,:) , sxyopw(:,:) )
[825]	280	DO jl = 1, jpl
	281	!--- ice volume ---
	282	CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0ice(:,:,jl) , sxice (:,:,jl) , &
[921]	283	sxxice(:,:,jl) , syice (:,:,jl) , &
	284	syyice(:,:,jl) , sxyice(:,:,jl) )
[825]	285	CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0ice(:,:,jl) , sxice (:,:,jl) , &
[921]	286	sxxice(:,:,jl) , syice (:,:,jl) , &
	287	syyice(:,:,jl) , sxyice(:,:,jl) )
[825]	288	!--- snow volume ---
	289	CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0sn (:,:,jl) , sxsn (:,:,jl) , &
[921]	290	sxxsn (:,:,jl) , sysn (:,:,jl) , &
	291	syysn (:,:,jl) , sxysn (:,:,jl) )
[825]	292	CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0sn (:,:,jl) , sxsn (:,:,jl) , &
[921]	293	sxxsn (:,:,jl) , sysn (:,:,jl) , &
	294	syysn (:,:,jl) , sxysn (:,:,jl) )
[825]	295	!--- ice salinity ---
	296	CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0sm (:,:,jl) , sxsal (:,:,jl) , &
[921]	297	sxxsal(:,:,jl) , sysal (:,:,jl) , &
	298	syysal(:,:,jl) , sxysal(:,:,jl) )
[825]	299	CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0sm (:,:,jl) , sxsal (:,:,jl) , &
[921]	300	sxxsal(:,:,jl) , sysal (:,:,jl) , &
	301	syysal(:,:,jl) , sxysal(:,:,jl) )
[825]	302	!--- ice age ---
	303	CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0oi (:,:,jl) , sxage (:,:,jl) , &
[921]	304	sxxage(:,:,jl) , syage (:,:,jl) , &
	305	syyage(:,:,jl) , sxyage(:,:,jl) )
[825]	306	CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0oi (:,:,jl) , sxage (:,:,jl) , &
[921]	307	sxxage(:,:,jl) , syage (:,:,jl) , &
	308	syyage(:,:,jl) , sxyage(:,:,jl) )
[825]	309	!--- ice concentration ---
	310	CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0a (:,:,jl) , sxa (:,:,jl) , &
[921]	311	sxxa (:,:,jl) , sya (:,:,jl) , &
	312	syya (:,:,jl) , sxya (:,:,jl) )
[825]	313	CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0a (:,:,jl) , sxa (:,:,jl) , &
[921]	314	sxxa (:,:,jl) , sya (:,:,jl) , &
	315	syya (:,:,jl) , sxya (:,:,jl) )
[825]	316	!--- ice / snow thermal energetic contents ---
	317	CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0c0 (:,:,jl) , sxc0 (:,:,jl) , &
[921]	318	sxxc0 (:,:,jl) , syc0 (:,:,jl) , &
	319	syyc0 (:,:,jl) , sxyc0 (:,:,jl) )
[825]	320	CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0c0 (:,:,jl) , sxc0 (:,:,jl) , &
[921]	321	sxxc0 (:,:,jl) , syc0 (:,:,jl) , &
	322	syyc0 (:,:,jl) , sxyc0 (:,:,jl) )
[825]	323	DO layer = 1, nlay_i
	324	CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0e(:,:,layer,jl) , &
[921]	325	sxe (:,:,layer,jl) , sxxe (:,:,layer,jl) , sye (:,:,layer,jl) , &
	326	syye (:,:,layer,jl), sxye (:,:,layer,jl) )
[825]	327	CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0e(:,:,layer,jl) , &
[921]	328	sxe (:,:,layer,jl) , sxxe (:,:,layer,jl) , sye (:,:,layer,jl) , &
	329	syye (:,:,layer,jl), sxye (:,:,layer,jl) )
[825]	330	END DO
	331
	332	END DO
	333	END DO
	334	ENDIF
	335
	336	!-------------------------------------------
	337	! Recover the properties from their contents
	338	!-------------------------------------------
	339
	340	zs0ow (:,:) = zs0ow(:,:) / area(:,:)
	341	DO jl = 1, jpl
	342	zs0ice(:,:,jl) = zs0ice(:,:,jl) / area(:,:)
	343	zs0sn (:,:,jl) = zs0sn (:,:,jl) / area(:,:)
	344	zs0sm (:,:,jl) = zs0sm (:,:,jl) / area(:,:)
	345	zs0oi (:,:,jl) = zs0oi (:,:,jl) / area(:,:)
	346	zs0a (:,:,jl) = zs0a (:,:,jl) / area(:,:)
	347	zs0c0 (:,:,jl) = zs0c0 (:,:,jl) / area(:,:)
	348	DO jk = 1, nlay_i
	349	zs0e(:,:,jk,jl) = zs0e(:,:,jk,jl) / area(:,:)
	350	END DO
	351	END DO
	352
[921]	353	!------------------------------------------------------------------------------!
	354	! 4) Diffusion of Ice fields
	355	!------------------------------------------------------------------------------!
[825]	356
[921]	357	!------------------------------------
	358	! 4.1) diffusion of open water area
	359	!------------------------------------
[825]	360
	361	! Compute total ice fraction
	362	zs0at(:,:) = 0.0
	363	DO jl = 1, jpl
	364	DO jj = 1, jpj
	365	DO ji = 1, jpi
	366	zs0at (ji,jj) = zs0at(ji,jj) + zs0a(ji,jj,jl) !
	367	END DO
[921]	368	END DO
[825]	369	END DO
	370
	371	! Masked eddy diffusivity coefficient at ocean U- and V-points
	372	DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row
	373	DO ji = 1 , fs_jpim1 ! vector opt.
	374	pahu(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0at(ji ,jj) ) ) ) &
	375	& * ( 1.0 - MAX( rzero, SIGN( rone, -zs0at(ji+1,jj) ) ) ) * ahiu(ji,jj)
	376	pahv(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0at(ji,jj ) ) ) ) &
	377	& * ( 1.0 - MAX( rzero, SIGN( rone,- zs0at(ji,jj+1) ) ) ) * ahiv(ji,jj)
	378	END DO !jj
	379	END DO ! ji
	380
	381	! Diffusion
	382	CALL lim_hdf( zs0ow (:,:) )
	383
[921]	384	!----------------------------------------
	385	! 4.2) Diffusion of other ice variables
	386	!----------------------------------------
[825]	387	DO jl = 1, jpl
	388
[921]	389	! Masked eddy diffusivity coefficient at ocean U- and V-points
[825]	390	DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row
	391	DO ji = 1 , fs_jpim1 ! vector opt.
	392	pahu(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji ,jj,jl) ) ) ) &
	393	& * ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)
	394	pahv(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji,jj,jl ) ) ) ) &
	395	& * ( 1.0 - MAX( rzero, SIGN( rone,- zs0a(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)
	396	END DO !jj
	397	END DO ! ji
	398
	399	CALL lim_hdf( zs0ice (:,:,jl) )
	400	CALL lim_hdf( zs0sn (:,:,jl) )
	401	CALL lim_hdf( zs0sm (:,:,jl) )
	402	CALL lim_hdf( zs0oi (:,:,jl) )
	403	CALL lim_hdf( zs0a (:,:,jl) )
	404	CALL lim_hdf( zs0c0 (:,:,jl) )
	405	DO jk = 1, nlay_i
	406	CALL lim_hdf( zs0e (:,:,jk,jl) )
	407	END DO ! jk
	408	END DO !jl
	409
[921]	410	!-----------------------------------------
	411	! 4.3) Remultiply everything by ice area
	412	!-----------------------------------------
[825]	413	zs0ow(:,:) = MAX(rzero, zs0ow(:,:) * area(:,:) )
	414	DO jl = 1, jpl
	415	zs0ice(:,:,jl) = MAX( rzero, zs0ice(:,:,jl) * area(:,:) ) !!bug: est-ce utile
	416	zs0sn (:,:,jl) = MAX( rzero, zs0sn (:,:,jl) * area(:,:) ) !!bug: cf /area juste apres
	417	zs0sm (:,:,jl) = MAX( rzero, zs0sm (:,:,jl) * area(:,:) ) !!bug: cf /area juste apres
	418	zs0oi (:,:,jl) = MAX( rzero, zs0oi (:,:,jl) * area(:,:) )
	419	zs0a (:,:,jl) = MAX( rzero, zs0a (:,:,jl) * area(:,:) ) !! suppress both change le resultat
	420	zs0c0 (:,:,jl) = MAX( rzero, zs0c0 (:,:,jl) * area(:,:) )
	421	DO jk = 1, nlay_i
	422	zs0e(:,:,jk,jl) = MAX( rzero, zs0e (:,:,jk,jl) * area(:,:) )
	423	END DO ! jk
	424	END DO ! jl
	425
[921]	426	!------------------------------------------------------------------------------!
	427	! 5) Update and limit ice properties after transport
	428	!------------------------------------------------------------------------------!
[825]	429
[921]	430	!--------------------------------------------------
	431	! 5.1) Recover mean values over the grid squares.
	432	!--------------------------------------------------
[825]	433
	434	DO jl = 1, jpl
	435	DO jk = 1, nlay_i
	436	DO jj = 1, jpj
	437	DO ji = 1, jpi
	438	zs0e (ji,jj,jk,jl) = &
[921]	439	MAX( rzero, zs0e (ji,jj,jk,jl) / area(ji,jj) )
[825]	440	END DO
	441	END DO
	442	END DO
	443	END DO
	444
	445	DO jj = 1, jpj
	446	DO ji = 1, jpi
	447	zs0ow (ji,jj) = MAX( rzero, zs0ow (ji,jj) / area(ji,jj) )
	448	END DO
	449	END DO
[921]	450
[825]	451	zs0at(:,:) = 0.0
	452	DO jl = 1, jpl
	453	DO jj = 1, jpj
	454	DO ji = 1, jpi
	455	zs0sn (ji,jj,jl) = MAX( rzero, zs0sn (ji,jj,jl)/area(ji,jj) )
	456	zs0ice(ji,jj,jl) = MAX( rzero, zs0ice(ji,jj,jl)/area(ji,jj) )
	457	zs0sm (ji,jj,jl) = MAX( rzero, zs0sm (ji,jj,jl)/area(ji,jj) )
	458	zs0oi (ji,jj,jl) = MAX( rzero, zs0oi (ji,jj,jl)/area(ji,jj) )
	459	zs0a (ji,jj,jl) = MAX( rzero, zs0a (ji,jj,jl)/area(ji,jj) )
	460	zs0c0 (ji,jj,jl) = MAX( rzero, zs0c0 (ji,jj,jl)/area(ji,jj) )
	461	zs0at (ji,jj) = zs0at(ji,jj) + zs0a(ji,jj,jl)
	462	END DO
	463	END DO
	464	END DO
	465
[921]	466	!---------------------------------------------------------
	467	! 5.2) Snow thickness, Ice thickness, Ice concentrations
	468	!---------------------------------------------------------
[825]	469
	470	DO jj = 1, jpj
	471	DO ji = 1, jpi
	472	zindb = MAX( 0.0 , SIGN( 1.0, zs0at(ji,jj) - zeps10) )
	473	zs0ow(ji,jj) = (1.0 - zindb) + zindb*MAX( zs0ow(ji,jj), 0.00 )
	474	ato_i(ji,jj) = zs0ow(ji,jj)
	475	END DO
	476	END DO
	477
	478	! Remove very small areas
	479	DO jl = 1, jpl
	480	DO jj = 1, jpj
	481	DO ji = 1, jpi
	482	zindb = MAX( 0.0 , SIGN( 1.0, zs0a(ji,jj,jl) - zeps10) )
	483
	484	zs0a(ji,jj,jl) = zindb * MIN( zs0a(ji,jj,jl), 0.99 )
	485	v_s(ji,jj,jl) = zindb * zs0sn (ji,jj,jl)
	486	v_i(ji,jj,jl) = zindb * zs0ice(ji,jj,jl)
	487
	488	zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - zeps10 ) )
	489	zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) )
	490	zindb = MAX( zindsn, zindic )
	491	zs0a (ji,jj,jl) = zindb * zs0a(ji,jj,jl) !ice concentration
	492	a_i (ji,jj,jl) = zs0a(ji,jj,jl)
	493	v_s (ji,jj,jl) = zindsn * v_s(ji,jj,jl)
	494	v_i (ji,jj,jl) = zindic * v_i(ji,jj,jl)
	495	END DO
	496	END DO
	497	END DO
	498
	499	DO jj = 1, jpj
	500	DO ji = 1, jpi
	501	zs0at(ji,jj) = SUM( zs0a(ji,jj,1:jpl) )
	502	END DO
	503	END DO
	504
[921]	505	!----------------------
	506	! 5.3) Ice properties
	507	!----------------------
[825]	508
	509	zbigval = 1.0d+13
	510
	511	DO jl = 1, jpl
	512	DO jj = 1, jpj
	513	DO ji = 1, jpi
	514
	515	! Switches and dummy variables
	516	zusvosn = 1.0/MAX( v_s(ji,jj,jl) , epsi16 )
	517	zusvoic = 1.0/MAX( v_i(ji,jj,jl) , epsi16 )
	518	zrtt = 173.15 * rone
	519	zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - zeps10 ) )
	520	zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) )
	521	zindb = MAX( zindsn, zindic )
	522
	523	! Ice salinity and age
[888]	524	zsal = MAX( MIN( (rhoic-rhosn)/rhoic*sss_m(ji,jj) , &
[921]	525	zusvoic * zs0sm(ji,jj,jl) ), s_i_min ) * &
	526	v_i(ji,jj,jl)
[825]	527	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &
	528	smv_i(ji,jj,jl) = zindiczsal + (1.0-zindic)0.0
	529
	530	zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / &
[921]	531	MAX( a_i(ji,jj,jl), epsi16 ) ), 0.0 ) * &
	532	a_i(ji,jj,jl)
[825]	533	oa_i (ji,jj,jl) = zindic*zage
	534
	535	! Snow heat content
	536	ze = MIN( MAX( 0.0, zs0c0(ji,jj,jl)*area(ji,jj) ), zbigval )
	537	e_s(ji,jj,1,jl) = zindsn * ze + (1.0 - zindsn) * 0.0
	538
	539	END DO !ji
	540	END DO !jj
	541	END DO ! jl
	542
	543	DO jl = 1, jpl
	544	DO jk = 1, nlay_i
	545	DO jj = 1, jpj
	546	DO ji = 1, jpi
	547	! Ice heat content
	548	zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) )
	549	ze = MIN( MAX( 0.0, zs0e(ji,jj,jk,jl)*area(ji,jj) ), zbigval )
	550	e_i(ji,jj,jk,jl) = zindic * ze + ( 1.0 - zindic ) * 0.0
	551	END DO !ji
	552	END DO ! jj
	553	END DO ! jk
	554	END DO ! jl
	555
	556	ENDIF
	557
[863]	558	IF(ln_ctl) THEN ! Control print
[867]	559	CALL prt_ctl_info(' ')
	560	CALL prt_ctl_info(' - Cell values : ')
	561	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
[863]	562	CALL prt_ctl(tab2d_1=area , clinfo1=' lim_trp : cell area :')
	563	CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_trp : at_i :')
	564	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_trp : vt_i :')
	565	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_trp : vt_s :')
	566	DO jl = 1, jpl
[867]	567	CALL prt_ctl_info(' ')
[863]	568	CALL prt_ctl_info(' - Category : ', ivar1=jl)
	569	CALL prt_ctl_info(' ~~~~~~~~~~')
	570	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_trp : a_i : ')
	571	CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_trp : ht_i : ')
	572	CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_trp : ht_s : ')
	573	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_trp : v_i : ')
	574	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_trp : v_s : ')
	575	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_trp : e_s : ')
	576	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_trp : t_su : ')
	577	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_trp : t_snow : ')
	578	CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_trp : sm_i : ')
	579	CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_trp : smv_i : ')
	580	DO jk = 1, nlay_i
[867]	581	CALL prt_ctl_info(' ')
[863]	582	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
	583	CALL prt_ctl_info(' ~~~~~~~')
	584	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_trp : t_i : ')
	585	CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_trp : e_i : ')
	586	END DO
	587	END DO
	588	ENDIF
	589
[825]	590	END SUBROUTINE lim_trp
	591
	592
	593	SUBROUTINE lim_trp_init
	594	!!-------------------------------------------------------------------
	595	!! * ROUTINE lim_trp_init *
	596	!!
	597	!! ** Purpose : initialization of ice advection parameters
	598	!!
	599	!! ** Method : Read the namicetrp namelist and check the parameter
	600	!! values called at the first timestep (nit000)
	601	!!
	602	!! ** input : Namelist namicetrp
	603	!!
	604	!! history :
	605	!! 2.0 ! 03-08 (C. Ethe) Original code
	606	!!-------------------------------------------------------------------
	607	NAMELIST/namicetrp/ bound
	608	!!-------------------------------------------------------------------
	609
	610	! Read Namelist namicetrp
	611	REWIND ( numnam_ice )
	612	READ ( numnam_ice , namicetrp )
	613	IF(lwp) THEN
	614	WRITE(numout,*)
	615	WRITE(numout,*) 'lim_trp_init : Ice parameters for advection '
	616	WRITE(numout,*) '~~~~~~~~~~~~'
	617	WRITE(numout,*) ' boundary conditions (0.0 no-slip, 1.0 free-slip) bound = ', bound
	618	WRITE(numout,*)
	619	ENDIF
[921]	620
[825]	621	END SUBROUTINE lim_trp_init
	622
	623	#else
	624	!!----------------------------------------------------------------------
	625	!! Default option Empty Module No sea-ice model
	626	!!----------------------------------------------------------------------
	627	CONTAINS
	628	SUBROUTINE lim_trp ! Empty routine
	629	END SUBROUTINE lim_trp
	630	#endif
	631
	632	!!======================================================================
	633	END MODULE limtrp

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

Download in other formats: