Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

wet_dry.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago
All wrk_alloc removed
File size: 17.5 KB

Rev	Line
[6152]	1	MODULE wet_dry
	2	!!==============================================================================
	3	!! * MODULE wet_dry *
	4	!! Wetting and drying includes initialisation routine and routines to
	5	!! compute and apply flux limiters and preserve water depth positivity
	6	!! only effects if wetting/drying is on (ln_wd == .true.)
	7	!!==============================================================================
[7646]	8	!! History : 3.6 ! 2014-09 ((H.Liu) Original code
[6152]	9	!! ! will add the runoff and periodic BC case later
	10	!!----------------------------------------------------------------------
	11
	12	!!----------------------------------------------------------------------
	13	!! wad_lmt : Compute the horizontal flux limiter and the limited velocity
	14	!! when wetting and drying happens
	15	!!----------------------------------------------------------------------
	16	USE oce ! ocean dynamics and tracers
	17	USE dom_oce ! ocean space and time domain
	18	USE sbc_oce, ONLY : ln_rnf ! surface boundary condition: ocean
	19	USE sbcrnf ! river runoff
	20	USE in_out_manager ! I/O manager
	21	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	22	USE lib_mpp ! MPP library
	23	USE timing ! Timing
	24
	25	IMPLICIT NONE
	26	PRIVATE
	27
	28	!!----------------------------------------------------------------------
	29	!! critical depths,filters, limiters,and masks for Wetting and Drying
	30	!! ---------------------------------------------------------------------
	31
	32	REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wdmask !: u- and v- limiter
[7646]	33	REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: ht_wd !: wetting and drying t-pt depths
	34	! (can include negative depths)
[6152]	35
	36	LOGICAL, PUBLIC :: ln_wd !: Wetting/drying activation switch (T:on,F:off)
	37	REAL(wp), PUBLIC :: rn_wdmin1 !: minimum water depth on dried cells
	38	REAL(wp), PUBLIC :: rn_wdmin2 !: tolerrance of minimum water depth on dried cells
	39	REAL(wp), PUBLIC :: rn_wdld !: land elevation below which wetting/drying
	40	!: will be considered
	41	INTEGER , PUBLIC :: nn_wdit !: maximum number of iteration for W/D limiter
	42
	43	PUBLIC wad_init ! initialisation routine called by step.F90
	44	PUBLIC wad_lmt ! routine called by sshwzv.F90
	45	PUBLIC wad_lmt_bt ! routine called by dynspg_ts.F90
	46
	47	!! * Substitutions
	48	# include "vectopt_loop_substitute.h90"
	49	CONTAINS
	50
	51	SUBROUTINE wad_init
	52	!!----------------------------------------------------------------------
	53	!! * ROUTINE wad_init *
	54	!!
	55	!! ** Purpose : read wetting and drying namelist and print the variables.
	56	!!
	57	!! ** input : - namwad namelist
	58	!!----------------------------------------------------------------------
	59	NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit
	60	INTEGER :: ios ! Local integer output status for namelist read
	61	INTEGER :: ierr ! Local integer status array allocation
	62	!!----------------------------------------------------------------------
	63
	64	REWIND( numnam_ref ) ! Namelist namwad in reference namelist
	65	! : Parameters for Wetting/Drying
	66	READ ( numnam_ref, namwad, IOSTAT = ios, ERR = 905)
	67	905 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in reference namelist', .TRUE.)
	68	REWIND( numnam_cfg ) ! Namelist namwad in configuration namelist
	69	! : Parameters for Wetting/Drying
	70	READ ( numnam_cfg, namwad, IOSTAT = ios, ERR = 906)
	71	906 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in configuration namelist', .TRUE. )
	72	IF(lwm) WRITE ( numond, namwad )
	73
	74	IF(lwp) THEN ! control print
	75	WRITE(numout,*)
	76	WRITE(numout,*) 'wad_init : Wetting and drying initialization through namelist read'
[7646]	77	WRITE(numout,*) '~~~~~~~~'
[6152]	78	WRITE(numout,*) ' Namelist namwad'
	79	WRITE(numout,*) ' Logical activation ln_wd = ', ln_wd
	80	WRITE(numout,*) ' Minimum wet depth on dried cells rn_wdmin1 = ', rn_wdmin1
	81	WRITE(numout,*) ' Tolerance of min wet depth rn_wdmin2 = ', rn_wdmin2
	82	WRITE(numout,*) ' land elevation threshold rn_wdld = ', rn_wdld
	83	WRITE(numout,*) ' Max iteration for W/D limiter nn_wdit = ', nn_wdit
[7646]	84	ENDIF
	85	!
[6152]	86	IF(ln_wd) THEN
[7646]	87	ALLOCATE( wdmask(jpi,jpj), ht_wd(jpi,jpj), STAT=ierr )
[6152]	88	IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error')
	89	ENDIF
[7646]	90	!
[6152]	91	END SUBROUTINE wad_init
	92
[7646]	93
[6152]	94	SUBROUTINE wad_lmt( sshb1, sshemp, z2dt )
	95	!!----------------------------------------------------------------------
	96	!! * ROUTINE wad_lmt *
	97	!!
	98	!! ** Purpose : generate flux limiters for wetting/drying
	99	!!
	100	!! ** Method : - Prevent negative depth occurring (Not ready for Agrif)
	101	!!
	102	!! ** Action : - calculate flux limiter and W/D flag
	103	!!----------------------------------------------------------------------
	104	REAL(wp), DIMENSION(:,:), INTENT(inout) :: sshb1
	105	REAL(wp), DIMENSION(:,:), INTENT(in) :: sshemp
	106	REAL(wp), INTENT(in) :: z2dt
	107	!
	108	INTEGER :: ji, jj, jk, jk1 ! dummy loop indices
[7646]	109	INTEGER :: jflag ! local scalar
[6152]	110	REAL(wp) :: zcoef, zdep1, zdep2 ! local scalars
	111	REAL(wp) :: zzflxp, zzflxn ! local scalars
	112	REAL(wp) :: zdepwd ! local scalar, always wet cell depth
	113	REAL(wp) :: ztmp ! local scalars
[7910]	114	REAL(wp), DIMENSION(jpi,jpj) :: zwdlmtu, zwdlmtv !: W/D flux limiters
	115	REAL(wp), DIMENSION(jpi,jpj) :: zflxp, zflxn ! local 2D workspace
	116	REAL(wp), DIMENSION(jpi,jpj) :: zflxu, zflxv ! local 2D workspace
	117	REAL(wp), DIMENSION(jpi,jpj) :: zflxu1, zflxv1 ! local 2D workspace
[6152]	118	!!----------------------------------------------------------------------
	119	!
	120
	121	IF( nn_timing == 1 ) CALL timing_start('wad_lmt')
	122
	123	IF(ln_wd) THEN
	124
	125	!
	126
	127	!IF(lwp) WRITE(numout,*)
	128	!IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting'
	129
[7646]	130	jflag = 0
[6152]	131	zdepwd = 50._wp !maximum depth on which that W/D could possibly happen
	132
	133
	134	zflxp(:,:) = 0._wp
	135	zflxn(:,:) = 0._wp
	136	zflxu(:,:) = 0._wp
	137	zflxv(:,:) = 0._wp
	138
	139	zwdlmtu(:,:) = 1._wp
	140	zwdlmtv(:,:) = 1._wp
	141
	142	! Horizontal Flux in u and v direction
	143	DO jk = 1, jpkm1
[7646]	144	DO jj = 1, jpj
	145	DO ji = 1, jpi
[6152]	146	zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
	147	zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
	148	END DO
	149	END DO
	150	END DO
	151
	152	zflxu(:,:) = zflxu(:,:) * e2u(:,:)
	153	zflxv(:,:) = zflxv(:,:) * e1v(:,:)
	154
[7646]	155	wdmask(:,:) = 1
	156	DO jj = 2, jpj
	157	DO ji = 2, jpi
[6152]	158
[7646]	159	IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE ! we don't care about land cells
	160	IF( ht_wd(ji,jj) > zdepwd ) CYCLE ! and cells which are unlikely to dry
[6152]	161
	162	zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj), 0._wp) + &
	163	& max(zflxv(ji,jj), 0._wp) - min(zflxv(ji, jj-1), 0._wp)
	164	zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj), 0._wp) + &
	165	& min(zflxv(ji,jj), 0._wp) - max(zflxv(ji, jj-1), 0._wp)
	166
[7646]	167	zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1
	168	IF(zdep2 .le. 0._wp) THEN !add more safty, but not necessary
	169	sshb1(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)
	170	IF(zflxu(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = 0._wp
	171	IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp
	172	IF(zflxv(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = 0._wp
	173	IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp
	174	wdmask(ji,jj) = 0._wp
[6152]	175	END IF
	176	ENDDO
	177	END DO
	178
	179
	180	!! start limiter iterations
	181	DO jk1 = 1, nn_wdit + 1
	182
	183
	184	zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
	185	zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
[7646]	186	jflag = 0 ! flag indicating if any further iterations are needed
[6152]	187
[7646]	188	DO jj = 2, jpj
	189	DO ji = 2, jpi
[6152]	190
[7646]	191	IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE
	192	IF( ht_wd(ji,jj) > zdepwd ) CYCLE
[6152]	193
	194	ztmp = e1e2t(ji,jj)
	195
	196	zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj), 0._wp) + &
	197	& max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji, jj-1), 0._wp)
	198	zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj), 0._wp) + &
	199	& min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji, jj-1), 0._wp)
	200
	201	zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
[7646]	202	zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj)
[6152]	203
[7646]	204	IF( zdep1 > zdep2 ) THEN
	205	wdmask(ji, jj) = 0
[6152]	206	zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
[7646]	207	!zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt )
	208	! flag if the limiter has been used but stop flagging if the only
	209	! changes have zeroed the coefficient since further iterations will
	210	! not change anything
	211	IF( zcoef > 0._wp ) THEN
	212	jflag = 1
	213	ELSE
	214	zcoef = 0._wp
	215	ENDIF
[6152]	216	IF(jk1 > nn_wdit) zcoef = 0._wp
	217	IF(zflxu1(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = zcoef
	218	IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
	219	IF(zflxv1(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = zcoef
[7646]	220	IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef
[6152]	221	END IF
	222	END DO ! ji loop
	223	END DO ! jj loop
	224
	225	CALL lbc_lnk( zwdlmtu, 'U', 1. )
	226	CALL lbc_lnk( zwdlmtv, 'V', 1. )
	227
[7646]	228	IF(lk_mpp) CALL mpp_max(jflag) !max over the global domain
[6152]	229
[7646]	230	IF(jflag == 0) EXIT
[6152]	231
	232	END DO ! jk1 loop
	233
	234	DO jk = 1, jpkm1
	235	un(:,:,jk) = un(:,:,jk) * zwdlmtu(:, :)
	236	vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:, :)
	237	END DO
	238
	239	CALL lbc_lnk( un, 'U', -1. )
	240	CALL lbc_lnk( vn, 'V', -1. )
[7646]	241	!
	242	un_b(:,:) = un_b(:,:) * zwdlmtu(:, :)
	243	vn_b(:,:) = vn_b(:,:) * zwdlmtv(:, :)
	244	CALL lbc_lnk( un_b, 'U', -1. )
	245	CALL lbc_lnk( vn_b, 'V', -1. )
[6152]	246
[7646]	247	IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!'
[6152]	248
	249	!IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field)
	250	!IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (update hdivn field)
	251	!
	252	!
[7646]	253	!
	254	ENDIF
[6152]	255	!
	256	IF( nn_timing == 1 ) CALL timing_stop('wad_lmt')
[7646]	257	!
[6152]	258	END SUBROUTINE wad_lmt
	259
[7646]	260
[6152]	261	SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rdtbt )
	262	!!----------------------------------------------------------------------
	263	!! * ROUTINE wad_lmt *
	264	!!
	265	!! ** Purpose : limiting flux in the barotropic stepping (dynspg_ts)
	266	!!
	267	!! ** Method : - Prevent negative depth occurring (Not ready for Agrif)
	268	!!
	269	!! ** Action : - calculate flux limiter and W/D flag
	270	!!----------------------------------------------------------------------
	271	REAL(wp), INTENT(in) :: rdtbt ! ocean time-step index
	272	REAL(wp), DIMENSION(:,:), INTENT(inout) :: zflxu, zflxv, sshn_e, zssh_frc
	273	!
	274	INTEGER :: ji, jj, jk, jk1 ! dummy loop indices
[7646]	275	INTEGER :: jflag ! local scalar
[6152]	276	REAL(wp) :: z2dt
	277	REAL(wp) :: zcoef, zdep1, zdep2 ! local scalars
	278	REAL(wp) :: zzflxp, zzflxn ! local scalars
	279	REAL(wp) :: zdepwd ! local scalar, always wet cell depth
	280	REAL(wp) :: ztmp ! local scalars
[7910]	281	REAL(wp), DIMENSION(jpi,jpj) :: zwdlmtu, zwdlmtv !: W/D flux limiters
	282	REAL(wp), DIMENSION(jpi,jpj) :: zflxp, zflxn ! local 2D workspace
	283	REAL(wp), DIMENSION(jpi,jpj) :: zflxu1, zflxv1 ! local 2D workspace
[6152]	284	!!----------------------------------------------------------------------
	285	!
	286	IF( nn_timing == 1 ) CALL timing_start('wad_lmt_bt')
	287
	288	IF(ln_wd) THEN
	289
	290	!
	291
	292	!IF(lwp) WRITE(numout,*)
	293	!IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting'
	294
[7646]	295	jflag = 0
[6152]	296	zdepwd = 50._wp !maximum depth that ocean cells can have W/D processes
	297
	298	z2dt = rdtbt
	299
	300	zflxp(:,:) = 0._wp
	301	zflxn(:,:) = 0._wp
	302
	303	zwdlmtu(:,:) = 1._wp
	304	zwdlmtv(:,:) = 1._wp
	305
	306	! Horizontal Flux in u and v direction
	307
[7646]	308	DO jj = 2, jpj
	309	DO ji = 2, jpi
[6152]	310
[7646]	311	IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE ! we don't care about land cells
	312	IF( ht_wd(ji,jj) > zdepwd ) CYCLE ! and cells which are unlikely to dry
[6152]	313
	314	zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj), 0._wp) + &
	315	& max(zflxv(ji,jj), 0._wp) - min(zflxv(ji, jj-1), 0._wp)
	316	zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj), 0._wp) + &
	317	& min(zflxv(ji,jj), 0._wp) - max(zflxv(ji, jj-1), 0._wp)
	318
[7646]	319	zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1
	320	IF(zdep2 .le. 0._wp) THEN !add more safety, but not necessary
	321	sshn_e(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)
	322	IF(zflxu(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = 0._wp
	323	IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp
	324	IF(zflxv(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = 0._wp
	325	IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp
[6152]	326	END IF
	327	ENDDO
	328	END DO
	329
	330
	331	!! start limiter iterations
	332	DO jk1 = 1, nn_wdit + 1
	333
	334
	335	zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
	336	zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
[7646]	337	jflag = 0 ! flag indicating if any further iterations are needed
[6152]	338
[7646]	339	DO jj = 2, jpj
	340	DO ji = 2, jpi
[6152]	341
[7646]	342	IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE
	343	IF( ht_wd(ji,jj) > zdepwd ) CYCLE
[6152]	344
	345	ztmp = e1e2t(ji,jj)
	346
	347	zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj), 0._wp) + &
	348	& max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji, jj-1), 0._wp)
	349	zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj), 0._wp) + &
	350	& min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji, jj-1), 0._wp)
	351
	352	zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
[7646]	353	zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj)
[6152]	354
	355	IF(zdep1 > zdep2) THEN
	356	zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
[7646]	357	!zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt )
	358	! flag if the limiter has been used but stop flagging if the only
	359	! changes have zeroed the coefficient since further iterations will
	360	! not change anything
	361	IF( zcoef > 0._wp ) THEN
	362	jflag = 1
	363	ELSE
	364	zcoef = 0._wp
	365	ENDIF
[6152]	366	IF(jk1 > nn_wdit) zcoef = 0._wp
	367	IF(zflxu1(ji, jj) > 0._wp) zwdlmtu(ji ,jj) = zcoef
	368	IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
	369	IF(zflxv1(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = zcoef
[7646]	370	IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef
[6152]	371	END IF
	372	END DO ! ji loop
	373	END DO ! jj loop
	374
	375	CALL lbc_lnk( zwdlmtu, 'U', 1. )
	376	CALL lbc_lnk( zwdlmtv, 'V', 1. )
	377
[7646]	378	IF(lk_mpp) CALL mpp_max(jflag) !max over the global domain
[6152]	379
[7646]	380	IF(jflag == 0) EXIT
[6152]	381
	382	END DO ! jk1 loop
	383
	384	zflxu(:,:) = zflxu(:,:) * zwdlmtu(:, :)
	385	zflxv(:,:) = zflxv(:,:) * zwdlmtv(:, :)
	386
	387	CALL lbc_lnk( zflxu, 'U', -1. )
	388	CALL lbc_lnk( zflxv, 'V', -1. )
	389
[7646]	390	IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!'
[6152]	391
	392	!IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field)
	393	!IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (update hdivn field)
	394	!
	395	!
[7646]	396	!
[6152]	397	END IF
	398
	399	IF( nn_timing == 1 ) CALL timing_stop('wad_lmt')
	400	END SUBROUTINE wad_lmt_bt
[7646]	401
	402	!!==============================================================================
[6152]	403	END MODULE wet_dry

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90 @ 7910

Download in other formats: