Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_initialize_real.F @ 3547

Last change on this file since 3547 was 11, checked in by aslmd, 14 years ago
spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 53.3 KB

Line
1	!REAL:MODEL_LAYER:INITIALIZATION
2
3	! This MODULE holds the routines which are used to perform various initializations
4	! for the individual domains, specifically for the Eulerian, mass-based coordinate.
5
6	!-----------------------------------------------------------------------
7
8	MODULE module_initialize
9
10	USE module_bc
11	USE module_configure
12	USE module_domain
13	USE module_io_domain
14	USE module_model_constants
15	! USE module_si_io_nmm
16	USE module_state_description
17	USE module_timing
18	USE module_soil_pre
19	#ifdef DM_PARALLEL
20	USE module_dm
21	#endif
22
23
24	CONTAINS
25
26	!-------------------------------------------------------------------
27
28	SUBROUTINE init_domain ( grid )
29
30	IMPLICIT NONE
31
32	! Input space and data. No gridded meteorological data has been stored, though.
33
34	! TYPE (domain), POINTER :: grid
35	TYPE (domain) :: grid
36
37	! Local data.
38
39	INTEGER :: dyn_opt
40	INTEGER :: idum1, idum2
41
42	#ifdef DEREF_KLUDGE
43	! see http://www.mmm.ucar.edu/wrf/WG2/topics/deref_kludge.htm
44	INTEGER :: sm31 , em31 , sm32 , em32 , sm33 , em33
45	INTEGER :: sm31x, em31x, sm32x, em32x, sm33x, em33x
46	INTEGER :: sm31y, em31y, sm32y, em32y, sm33y, em33y
47	#endif
48
49	#include "deref_kludge.h"
50
51	CALL nl_get_dyn_opt ( head_grid%id, dyn_opt )
52
53	CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )
54
55	IF ( dyn_opt .eq. 1 &
56	.or. dyn_opt .eq. 2 &
57	.or. dyn_opt .eq. 3 &
58	) THEN
59	CALL wrf_error_fatal ( "no RK version within dyn_nmm, dyn_opt wrong in namelist, wrf_error_fataling" )
60
61	ELSEIF ( dyn_opt .eq. 4 ) THEN
62
63	CALL init_domain_nmm (grid &
64	!
65	#include <nmm_actual_args.inc>
66	!
67	)
68
69	ELSE
70	WRITE(0,*)' init_domain: unknown or unimplemented dyn_opt = ',dyn_opt
71	CALL wrf_error_fatal ( "ERROR-dyn_opt-wrong-in-namelist" )
72	ENDIF
73
74	END SUBROUTINE init_domain
75
76	!-------------------------------------------------------------------
77	!---------------------------------------------------------------------
78	SUBROUTINE init_domain_nmm ( grid &
79	!
80	# include <nmm_dummy_args.inc>
81	!
82	)
83
84	USE module_optional_si_input
85	IMPLICIT NONE
86
87	! Input space and data. No gridded meteorological data has been stored, though.
88
89	! TYPE (domain), POINTER :: grid
90	TYPE (domain) :: grid
91
92	# include <nmm_dummy_decl.inc>
93
94	TYPE (grid_config_rec_type) :: config_flags
95
96	! Local domain indices and counters.
97
98	INTEGER :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
99
100	INTEGER :: &
101	ids, ide, jds, jde, kds, kde, &
102	ims, ime, jms, jme, kms, kme, &
103	its, ite, jts, jte, kts, kte, &
104	i, j, k, NNXP, NNYP, ICOUNT, &
105	FLAG_RUC_SNOW, FLAG_NAM_SNOW
106
107	! Local data
108
109	CHARACTER(LEN=19):: start_date
110
111	#ifdef DM_PARALLEL
112
113	LOGICAL,EXTERNAL :: WRF_DM_ON_MONITOR
114
115	! INTEGER :: DOMDESC
116	REAL,ALLOCATABLE :: SICE_G(:,:), SM_G(:,:)
117	INTEGER, ALLOCATABLE:: IHE_G(:),IHW_G(:)
118	#endif
119
120
121	CHARACTER (LEN=132) :: message
122
123	INTEGER :: error
124	REAL :: p_surf, p_level
125	REAL :: cof1, cof2
126	REAL :: qvf , qvf1 , qvf2 , pd_surf
127	REAL :: p00 , t00 , a
128	REAL :: hold_znw, rmin,rmax
129
130	LOGICAL :: stretch_grid, dry_sounding, debug, log_flag_sst
131
132	REAL, ALLOCATABLE,DIMENSION(:,:):: ADUM2D,SNOWC,HT,TG_ALT
133
134	INTEGER, ALLOCATABLE, DIMENSION(:):: KHL2,KVL2,KHH2,KVH2, &
135	KHLA,KHHA,KVLA,KVHA
136
137	! INTEGER, ALLOCATABLE, DIMENSION(:,:):: LU_INDEX
138
139	REAL, ALLOCATABLE, DIMENSION(:):: DXJ,WPDARJ,CPGFUJ,CURVJ, &
140	FCPJ,FDIVJ,EMJ,EMTJ,FADJ, &
141	HDACJ,DDMPUJ,DDMPVJ
142
143	!-- Carsel and Parrish [1988]
144	REAL , DIMENSION(100) :: lqmi
145	integer iicount
146
147	REAL:: TPH0D,TLM0D
148	REAL:: TPH0,WB,SB,DLM,DPH,TDLM,TDPH
149	REAL:: WBI,SBI,EBI,ANBI,STPH0,CTPH0
150	REAL:: TSPH,DTAD,DTCF
151	REAL:: ACDT,CDDAMP,TPH,DXP,TLM,FP
152	REAL:: CTPH,STPH
153	REAL:: WBD,SBD
154	REAL:: RSNOW,SNOFAC
155	REAL, PARAMETER:: SALP=2.60
156	REAL, PARAMETER:: SNUP=0.040
157	REAL:: SMCSUM,STCSUM,SEAICESUM,FISX
158	REAL:: cur_smc, aposs_smc
159
160	REAL:: TERM1,APH
161
162	INTEGER:: KHH,KVH,JAM,JA, IHL, IHH, L
163	INTEGER:: II,JJ,ISRCH,ISUM,ITER
164
165	REAL, PARAMETER:: DTR=0.01745329
166	REAL, PARAMETER:: W_NMM=0.08
167	REAL, PARAMETER:: COAC=1.6
168	REAL, PARAMETER:: CODAMP=6.4
169	REAL, PARAMETER:: TWOM=.00014584
170	REAL, PARAMETER:: CP=1004.6
171	REAL, PARAMETER:: DFC=1.0
172	REAL, PARAMETER:: DDFC=8.0
173	REAL, PARAMETER:: ROI=916.6
174	REAL, PARAMETER:: R=287.04
175	REAL, PARAMETER:: CI=2060.0
176	REAL, PARAMETER:: ROS=1500.
177	REAL, PARAMETER:: CS=1339.2
178	REAL, PARAMETER:: DS=0.050
179	REAL, PARAMETER:: AKS=.0000005
180	REAL, PARAMETER:: DZG=2.85
181	REAL, PARAMETER:: DI=.1000
182	REAL, PARAMETER:: AKI=0.000001075
183	REAL, PARAMETER:: DZI=2.0
184	REAL, PARAMETER:: THL=210.
185	REAL, PARAMETER:: PLQ=70000.
186	REAL, PARAMETER:: ERAD=6371200.
187	REAL, PARAMETER:: TG0=258.16
188	REAL, PARAMETER:: TGA=30.0
189
190
191	#ifdef DEREF_KLUDGE
192	! see http://www.mmm.ucar.edu/wrf/WG2/topics/deref_kludge.htm
193	INTEGER :: sm31 , em31 , sm32 , em32 , sm33 , em33
194	INTEGER :: sm31x, em31x, sm32x, em32x, sm33x, em33x
195	INTEGER :: sm31y, em31y, sm32y, em32y, sm33y, em33y
196	#endif
197
198	#include "deref_kludge.h"
199
200	if (ALLOCATED(ADUM2D)) DEALLOCATE(ADUM2D)
201	if (ALLOCATED(TG_ALT)) DEALLOCATE(TG_ALT)
202
203	#define COPY_IN
204	#include <nmm_scalar_derefs.inc>
205	#ifdef DM_PARALLEL
206	# include <nmm_data_calls.inc>
207	#endif
208
209	SELECT CASE ( model_data_order )
210	CASE ( DATA_ORDER_ZXY )
211	kds = grid%sd31 ; kde = grid%ed31 ;
212	ids = grid%sd32 ; ide = grid%ed32 ;
213	jds = grid%sd33 ; jde = grid%ed33 ;
214
215	kms = grid%sm31 ; kme = grid%em31 ;
216	ims = grid%sm32 ; ime = grid%em32 ;
217	jms = grid%sm33 ; jme = grid%em33 ;
218
219	kts = grid%sp31 ; kte = grid%ep31 ; ! tile is entire patch
220	its = grid%sp32 ; ite = grid%ep32 ; ! tile is entire patch
221	jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch
222
223	CASE ( DATA_ORDER_XYZ )
224	ids = grid%sd31 ; ide = grid%ed31 ;
225	jds = grid%sd32 ; jde = grid%ed32 ;
226	kds = grid%sd33 ; kde = grid%ed33 ;
227
228	ims = grid%sm31 ; ime = grid%em31 ;
229	jms = grid%sm32 ; jme = grid%em32 ;
230	kms = grid%sm33 ; kme = grid%em33 ;
231
232	its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch
233	jts = grid%sp32 ; jte = grid%ep32 ; ! tile is entire patch
234	kts = grid%sp33 ; kte = grid%ep33 ; ! tile is entire patch
235
236	CASE ( DATA_ORDER_XZY )
237	ids = grid%sd31 ; ide = grid%ed31 ;
238	kds = grid%sd32 ; kde = grid%ed32 ;
239	jds = grid%sd33 ; jde = grid%ed33 ;
240
241	ims = grid%sm31 ; ime = grid%em31 ;
242	kms = grid%sm32 ; kme = grid%em32 ;
243	jms = grid%sm33 ; jme = grid%em33 ;
244
245	its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch
246	kts = grid%sp32 ; kte = grid%ep32 ; ! tile is entire patch
247	jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch
248
249	END SELECT
250
251
252	grid%DT=float(grid%TIME_STEP)
253
254	NNXP=min(ITE,IDE-1)
255	NNYP=min(JTE,JDE-1)
256
257	write(0,*) 'nnxp,nnyp: ', nnxp,nnyp
258	write(0,*) 'IDE, JDE: ', IDE,JDE
259
260	JAM=6+2*(JDE-JDS-10)
261
262	ALLOCATE(ADUM2D(grid%sm31:grid%em31,jms:jme))
263	ALLOCATE(KHL2(JTS:NNYP),KVL2(JTS:NNYP),KHH2(JTS:NNYP),KVH2(JTS:NNYP))
264	ALLOCATE(DXJ(JTS:NNYP),WPDARJ(JTS:NNYP),CPGFUJ(JTS:NNYP),CURVJ(JTS:NNYP))
265	ALLOCATE(FCPJ(JTS:NNYP),FDIVJ(JTS:NNYP),&
266	FADJ(JTS:NNYP))
267	ALLOCATE(HDACJ(JTS:NNYP),DDMPUJ(JTS:NNYP),DDMPVJ(JTS:NNYP))
268	ALLOCATE(KHLA(JAM),KHHA(JAM))
269	ALLOCATE(KVLA(JAM),KVHA(JAM))
270
271
272	CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
273
274	write(0,*) 'cen_lat: ', config_flags%cen_lat
275	write(0,*) 'cen_lon: ', config_flags%cen_lon
276	write(0,*) 'dx: ', config_flags%dx
277	write(0,*) 'dy: ', config_flags%dy
278	write(0,*) 'config_flags%start_year: ', config_flags%start_year
279	write(0,*) 'config_flags%start_month: ', config_flags%start_month
280	write(0,*) 'config_flags%start_day: ', config_flags%start_day
281	write(0,*) 'config_flags%start_hour: ', config_flags%start_hour
282
283	write(0,*) 'writing to start_date'
284
285	write(start_date,435) config_flags%start_year, config_flags%start_month, &
286	config_flags%start_day, config_flags%start_hour
287	435 format(I4,'-',I2.2,'-',I2.2,'_',I2.2,':00:00')
288
289	dlmd=config_flags%dx
290	dphd=config_flags%dy
291	tph0d=config_flags%cen_lat
292	tlm0d=config_flags%cen_lon
293
294	write(0,*) 'dlmd, dphd, tph0d, tlm0d: ', dlmd, dphd, tph0d, tlm0d
295
296	!==========================================================================
297
298	!!
299
300	! Check to see if the boundary conditions are set
301	! properly in the namelist file.
302	! This checks for sufficiency and redundancy.
303
304	CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )
305
306	! Some sort of "this is the first time" initialization. Who knows.
307
308	grid%itimestep=0
309
310	! Pull in the info in the namelist to compare it to the input data.
311
312	grid%real_data_init_type = model_config_rec%real_data_init_type
313
314	!!! WEASD has "snow water equivalent" in mm
315
316	FLAG_RUC_SNOW=0
317	FLAG_NAM_SNOW=0
318	IF(maxval(snow).gt.0.) FLAG_RUC_SNOW=1
319	IF(maxval(weasd).gt.0.) FLAG_NAM_SNOW=1
320
321	DO j = jts, MIN(jte,jde-1)
322	DO i = its, MIN(ite,ide-1)
323
324	IF(SM(I,J).GT.0.9) THEN
325
326	IF (XICE(I,J) .gt. 0) then
327	SI(I,J)=1.0
328	ENDIF
329
330	! SEA
331	EPSR(I,J)=.97
332	GFFC(I,J)=0.
333	ALBEDO(I,J)=.06
334	ALBASE(I,J)=.06
335	IF(SI (I,J).GT.0. ) THEN
336	! SEA-ICE
337	SM(I,J)=0.
338	SI(I,J)=0.
339	SICE(I,J)=1.
340	GFFC(I,J)=0. ! just leave zero as irrelevant
341	ALBEDO(I,J)=.60
342	ALBASE(I,J)=.60
343	ENDIF
344	ELSE
345
346	! LAND
347	IF(FLAG_RUC_SNOW .eq. 1) THEN
348	! RUC variable for snow is SNOW
349	SI(I,J)=SNOW(I,J)/RHOSN(I,J)
350	EPSR(I,J)=1.0
351	GFFC(I,J)=0.0 ! just leave zero as irrelevant
352	SICE(I,J)=0.
353	SNO(I,J)=SNOW(I,J)/1000.
354
355	ELSEIF (FLAG_NAM_SNOW .eq. 1) THEN
356	! NAM variable for snow is WEASD
357	SI(I,J)=5.0*WEASD(I,J)/1000.
358	EPSR(I,J)=1.0
359	GFFC(I,J)=0.0 ! just leave zero as irrelevant
360	SICE(I,J)=0.
361	SNO(I,J)=WEASD(I,J)/1000.
362	ELSE
363	SI(I,J)=0.
364	EPSR(I,J)=1.0
365	GFFC(I,J)=0.0 ! just leave zero as irrelevant
366	SICE(I,J)=0.
367	SNO(I,J)=0.
368	ENDIF
369	ENDIF
370	ENDDO
371	ENDDO
372
373	! DETERMINE ALBEDO OVER LAND
374	DO j = jts, MIN(jte,jde-1)
375	DO i = its, MIN(ite,ide-1)
376	IF(SM(I,J).LT.0.9.AND.SICE(I,J).LT.0.9) THEN
377	! SNOWFREE ALBEDO
378	IF ( (SNO(I,J) .EQ. 0.0) .OR. &
379	(ALBASE(I,J) .GE. MXSNAL(I,J) ) ) THEN
380	ALBEDO(I,J) = ALBASE(I,J)
381	ELSE
382	! MODIFY ALBEDO IF SNOWCOVER:
383	! BELOW SNOWDEPTH THRESHOLD...
384	IF (SNO(I,J) .LT. SNUP) THEN
385	RSNOW = SNO(I,J)/SNUP
386	SNOFAC = 1. - ( EXP(-SALPRSNOW) - RSNOWEXP(-SALP))
387	! ABOVE SNOWDEPTH THRESHOLD...
388	ELSE
389	SNOFAC = 1.0
390	ENDIF
391	! CALCULATE ALBEDO ACCOUNTING FOR SNOWDEPTH AND VGFRCK
392	ALBEDO(I,J) = ALBASE(I,J) &
393	+ (1.0-VEGFRA(I,J))SNOFAC(MXSNAL(I,J)-ALBASE(I,J))
394	ENDIF
395	END IF
396	IF(FLAG_RUC_SNOW .eq. 1) then
397	SI(I,J)=SNOW(I,J)/RHOSN(I,J)*1000.
398	SNO(I,J)=SNOW(I,J)
399	! SNO(I,J)=SNOW(I,J)
400	! SNO(I,J)=WEASD(I,J)
401	ELSEIF (FLAG_NAM_SNOW .eq. 1) then
402
403	SI(I,J)=5.0*WEASD(I,J)
404	SNO(I,J)=WEASD(I,J)
405	ELSE
406	SI(I,J)=0.
407	SNO(I,J)=0.
408	ENDIF
409	ENDDO
410	ENDDO
411
412	#ifdef DM_PARALLEL
413
414	ALLOCATE(SM_G(IDS:IDE,JDS:JDE),SICE_G(IDS:IDE,JDS:JDE))
415
416	CALL WRF_PATCH_TO_GLOBAL_REAL( SICE(IMS,JMS) &
417	&, SICE_G,grid%DOMDESC &
418	&, 'z','xy' &
419	&, IDS,IDE-1,JDS,JDE-1,1,1 &
420	&, IMS,IME,JMS,JME,1,1 &
421	&, ITS,ITE,JTS,JTE,1,1 )
422
423	CALL WRF_PATCH_TO_GLOBAL_REAL( SM(IMS,JMS) &
424	&, SM_G,grid%DOMDESC &
425	&, 'z','xy' &
426	&, IDS,IDE-1,JDS,JDE-1,1,1 &
427	&, IMS,IME,JMS,JME,1,1 &
428	&, ITS,ITE,JTS,JTE,1,1 )
429
430
431	IF (WRF_DM_ON_MONITOR()) THEN
432
433	637 format(40(f3.0,1x))
434
435	allocate(IHE_G(JDS:JDE-1),IHW_G(JDS:JDE-1))
436	DO j = JDS, JDE-1
437	IHE_G(J)=MOD(J+1,2)
438	IHW_G(J)=IHE_G(J)-1
439	ENDDO
440
441	DO ITER=1,10
442	DO j = jds+1, (jde-1)-1
443	DO i = ids+1, (ide-1)-1
444
445	! any sea ice around point in question?
446
447	IF (SM_G(I,J) .eq. 1.) THEN
448	SEAICESUM=SICE_G(I+IHE_G(J),J+1)+SICE_G(I+IHW_G(J),J+1)+ &
449	SICE_G(I+IHE_G(J),J-1)+SICE_G(I+IHW_G(J),J-1)
450	IF (SEAICESUM .ge. 1. .and. SEAICESUM .lt. 3.) THEN
451	IF ((SICE_G(I+IHE_G(J),J+1).eq.0 .and. SM_G(I+IHE_G(J),J+1).eq.0) .OR. &
452	(SICE_G(I+IHW_G(J),J+1).eq.0 .and. SM_G(I+IHW_G(J),J+1).eq.0) .OR. &
453	(SICE_G(I+IHE_G(J),J-1).eq.0 .and. SM_G(I+IHE_G(J),J-1).eq.0) .OR. &
454	(SICE_G(I+IHW_G(J),J-1).eq.0 .and. SM_G(I+IHW_G(J),J-1).eq.0)) THEN
455
456	! HAVE SEA ICE AND A SURROUNDING LAND POINT - CONVERT TO SEA ICE
457	! write(0,*) 'making seaice (1): ', I,J
458	SICE_G(I,J)=1.0
459	SM_G(I,J)=0.
460	ENDIF
461	ELSEIF (SEAICESUM .ge. 3) THEN
462	! WATER POINT SURROUNDED BY ICE - CONVERT TO SEA ICE
463	! write(0,*) 'making seaice (2): ', I,J
464	SICE_G(I,J)=1.0
465	SM_G(I,J)=0.
466	ENDIF
467	ENDIF
468
469	ENDDO
470	ENDDO
471	ENDDO
472
473	ENDIF
474
475	CALL WRF_GLOBAL_TO_PATCH_REAL( SICE_G, SICE &
476	&, grid%DOMDESC &
477	&, 'z','xy' &
478	&, IDS,IDE-1,JDS,JDE-1,1,1 &
479	&, IMS,IME,JMS,JME,1,1 &
480	&, ITS,ITE,JTS,JTE,1,1 )
481
482	CALL WRF_GLOBAL_TO_PATCH_REAL( SM_G,SM &
483	&, grid%DOMDESC &
484	&, 'z','xy' &
485	&, IDS,IDE-1,JDS,JDE-1,1,1 &
486	&, IMS,IME,JMS,JME,1,1 &
487	&, ITS,ITE,JTS,JTE,1,1 )
488
489	IF (WRF_DM_ON_MONITOR()) THEN
490
491	DEALLOCATE(SM_G,SICE_G)
492	DEALLOCATE(IHE_G,IHW_G)
493
494	ENDIF
495
496	! write(0,*) 'revised sea ice on patch'
497	! do J=JTE,JTS,-JTE/25
498	! write(0,637) (SICE(I,J),I=ITS,ITE,ITE/20)
499	! enddo
500
501	! write(0,*) 'revised sea mask on patch'
502	! do J=JTE,JTS,-JTE/25
503	! write(0,637) (SM(I,J),I=ITS,ITE,ITE/20)
504	! enddo
505
506
507	#else
508	! serial sea ice reprocessing
509
510	DO j = jts, MIN(jte,jde-1)
511	IHE(J)=MOD(J+1,2)
512	IHW(J)=IHE(J)-1
513	ENDDO
514
515	DO ITER=1,10
516	DO j = jts+1, MIN(jte,jde-1)-1
517	DO i = its+1, MIN(ite,ide-1)-1
518
519	! any sea ice around point in question?
520
521	IF (SM(I,J) .eq. 1) THEN
522	SEAICESUM=SICE(I+IHE(J),J+1)+SICE(I+IHW(J),J+1)+ &
523	SICE(I+IHE(J),J-1)+SICE(I+IHW(J),J-1)
524	IF (SEAICESUM .ge. 1. .and. SEAICESUM .lt. 3.) THEN
525	IF ((SICE(I+IHE(J),J+1).eq.0 .and. SM(I+IHE(J),J+1).eq.0) .OR. &
526	(SICE(I+IHW(J),J+1).eq.0 .and. SM(I+IHW(J),J+1).eq.0) .OR. &
527	(SICE(I+IHE(J),J-1).eq.0 .and. SM(I+IHE(J),J-1).eq.0) .OR. &
528	(SICE(I+IHW(J),J-1).eq.0 .and. SM(I+IHW(J),J-1).eq.0)) THEN
529
530	! HAVE SEA ICE AND A SURROUNDING LAND POINT - CONVERT TO SEA ICE
531	SICE(I,J)=1.0
532	SM(I,J)=0.
533	ENDIF
534	ELSEIF (SEAICESUM .ge. 3) THEN
535	! WATER POINT SURROUNDED BY ICE - CONVERT TO SEA ICE
536	SICE(I,J)=1.0
537	SM(I,J)=0.
538	ENDIF
539	ENDIF
540
541	ENDDO
542	ENDDO
543	ENDDO
544
545	#endif
546
547	! this block meant to guarantee land/sea agreement between SM and landmask
548
549	DO j = jts, MIN(jte,jde-1)
550	DO i = its, MIN(ite,ide-1)
551
552	if (SM(I,J) .gt. 0.5) then
553	landmask(I,J)=0.0
554	elseif (SM(I,J) .eq. 0 .and. SICE(I,J) .eq. 1) then
555	landmask(I,J)=0.0
556	elseif (SM(I,J) .lt. 0.5 .and. SICE(I,J) .eq. 0) then
557	landmask(I,J)=1.0
558	else
559	write(0,*) 'missed point in landmask definition ' , I,J
560	landmask(I,J)=0.0
561	endif
562
563	ENDDO
564	ENDDO
565
566	! For sf_surface_physics = 1, we want to use close to a 10 cm value
567	! for the bottom level of the soil temps.
568
569	IF ( ( model_config_rec%sf_surface_physics(grid%id) .EQ. 1 ) .AND. &
570	( flag_st000010 .EQ. 1 ) ) THEN
571	DO j = jts , MIN(jde-1,jte)
572	DO i = its , MIN(ide-1,ite)
573	soiltb(i,j) = st000010(i,j)
574	END DO
575	END DO
576	END IF
577
578	! Adjust the various soil temperature values depending on the difference in
579	! in elevation between the current model's elevation and the incoming data's
580	! orography.
581
582	write(0,*) 'flag_toposoil= ', flag_toposoil
583
584	IF ( ( flag_toposoil .EQ. 1 ) ) THEN
585
586	ALLOCATE(HT(ims:ime,jms:jme))
587
588	DO J=jms,jme
589	DO I=ims,ime
590	HT(I,J)=FIS(I,J)/9.81
591	END DO
592	END DO
593
594	! if (maxval(toposoil) .gt. 100.) then
595	!
596	! Being avoided. Something to revisit eventually.
597	!
598	!1219 might be simply a matter of including TOPOSOIL
599	!
600	! CODE NOT TESTED AT NCEP USING THIS FUNCTIONALITY,
601	! SO TO BE SAFE WILL AVOID FOR RETRO RUNS.
602	!
603	! write(0,*) 'calling adjust_soil_temp_new'
604	! CALL adjust_soil_temp_new ( soiltb , 2 , &
605	! nmm_tsk , ht , toposoil , landmask, flag_toposoil , &
606	! st000010 , st010040 , st040100 , st100200 , st010200 , &
607	! flag_st000010 , flag_st010040 , flag_st040100 , &
608	! flag_st100200 , flag_st010200 , &
609	! soilt000 , soilt005 , soilt020 , soilt040 , &
610	! soilt160 , soilt300 , &
611	! flag_soilt000 , flag_soilt005 , flag_soilt020 , &
612	! flag_soilt040 , flag_soilt160 , flag_soilt300 , &
613	! ids , ide , jds , jde , kds , kde , &
614	! ims , ime , jms , jme , kms , kme , &
615	! its , ite , jts , jte , kts , kte )
616	! endif
617
618	END IF
619
620	! Process the LSM data.
621
622	IF ( grid%real_data_init_type .EQ. 1 ) THEN
623
624	num_veg_cat = SIZE ( landusef , DIM=2 )
625	num_soil_top_cat = SIZE ( soilctop , DIM=2 )
626	num_soil_bot_cat = SIZE ( soilcbot , DIM=2 )
627
628	! sm (1=water, 0=land)
629	! landmask(0=water, 1=land)
630
631	CALL process_percent_cat_new ( landmask , &
632	landusef , soilctop , soilcbot , &
633	isltyp , ivgtyp , &
634	num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
635	ids , ide , jds , jde , kds , kde , &
636	ims , ime , jms , jme , kms , kme , &
637	its , ite , jts , jte , kts , kte , &
638	model_config_rec%iswater(grid%id) )
639
640	DO j = jts, MIN(jde-1,jte)
641	DO i = its, MIN(ide-1,ite)
642
643	IF (SICE(I,J) .eq. 0) THEN
644
645	if (landmask(I,J) .gt. 0.5 .and. sm(I,J) .eq. 1.0) then
646	write(0,*) 'land mask and SM both > 0.5: ', &
647	I,J,landmask(I,J),sm(I,J)
648
649	SM(I,J)=0.
650
651	elseif (landmask(I,J) .lt. 0.5 .and. sm(I,J) .eq. 0.0) then
652	write(0,*) 'land mask and SM both < 0.5: ', &
653	I,J, landmask(I,J),sm(I,J)
654
655	SM(I,J)=1.
656
657	endif
658
659	ELSE
660
661	if (landmask(I,J) .gt. 0.5 .and. SM(I,J)+SICE(I,J) .eq. 1) then
662	write(0,*) 'landmask says LAND, SM/SICE say SEAICE: ', I,J
663	endif
664
665	ENDIF
666
667	ENDDO
668	ENDDO
669
670	DO j = jts, MIN(jde-1,jte)
671	DO i = its, MIN(ide-1,ite)
672
673	if (SICE(I,J) .eq. 1.0) then
674	!!!! change vegtyp and sltyp to fit seaice (desireable??)
675	ISLTYP(I,J)=16
676	IVGTYP(I,J)=24
677	endif
678
679	ENDDO
680	ENDDO
681
682
683	! MOVE HERE
684
685	write(0,*) 'flag_sst before define is: ', flag_sst
686	FLAG_SST=1
687
688	DO j = jts, MIN(jde-1,jte)
689	DO i = its, MIN(ide-1,ite)
690
691	if (SM(I,J) .lt. 0.5) then
692	SST(I,J)=0.
693	endif
694
695	if (SM(I,J) .gt. 0.5) then
696	if (SST(I,J) .eq. 0) then
697	SST(I,J)=NMM_TSK(I,J)
698	endif
699	NMM_TSK(I,J)=0.
700	endif
701
702
703	if ( (NMM_TSK(I,J)+SST(I,J)) .lt. 200. .or. &
704	(NMM_TSK(I,J)+SST(I,J)) .gt. 350. ) then
705	write(0,*) 'TSK, SST trouble at : ', I,J
706	write(0,*) 'SM= ', SM(I,J)
707	write(0,*) 'NMM_TSK(I,J), SST(I,J): ', NMM_TSK(I,J), SST(I,J)
708	endif
709
710	ENDDO
711	ENDDO
712
713	write(0,*) 'SM'
714	do J=min(jde-1,jte),jts,-(jte-jts)/15
715	write(0,635) (sm(i,J),I=its,ite,(ite-its)/10)
716	enddo
717
718	! write(0,*) 'SST/NMM_TSK'
719	! do J=min(jde-1,jte),jts,-(jte-jts)/15
720	! write(0,635) (SST(I,J)+NMM_TSK(I,J),I=ITS,min(ide-1,ite),(ite-its)/10)
721	! enddo
722
723	635 format(20(f5.1,1x))
724
725	DO j = jts, MIN(jde-1,jte)
726	DO i = its, MIN(ide-1,ite)
727	IF ( ( landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) ) THEN
728	soiltb(i,j) = sst(i,j)
729	!curious ELSE IF ( landmask(i,j) .LT. 0.5 ) THEN
730	ELSE IF ( landmask(i,j) .GT. 0.5 ) THEN
731	soiltb(i,j) = nmm_tsk(i,j)
732	END IF
733	END DO
734	END DO
735
736	! END IF
737
738	! END MOVE HERE
739
740	! Land use categories, dominant soil and vegetation types (if available).
741
742	! allocate(lu_index(ims:ime,jms:jme))
743
744	DO j = jts, MIN(jde-1,jte)
745	DO i = its, MIN(ide-1,ite)
746	lu_index(i,j) = ivgtyp(i,j)
747	END DO
748	END DO
749
750	END IF
751
752	if (flag_sst .eq. 1) log_flag_sst=.true.
753	if (flag_sst .eq. 0) log_flag_sst=.false.
754
755	write(0,*) 'st_input dimensions: ', size(st_input,dim=1),size(st_input,dim=2),size(st_input,dim=3)
756
757	write(0,*) 'maxval st_input(1): ', maxval(st_input(:,1,:))
758	write(0,*) 'maxval st_input(2): ', maxval(st_input(:,2,:))
759	write(0,*) 'maxval st_input(3): ', maxval(st_input(:,3,:))
760	write(0,*) 'maxval st_input(4): ', maxval(st_input(:,4,:))
761
762	!!!!!!!!!!!!!!!!!!!!!!!!!
763
764	ALLOCATE(TG_ALT(grid%sm31:grid%em31,jms:jme))
765
766	TPH0=TPH0D*DTR
767	WBD=-(((ide-1)-1)*DLMD)
768	WB= WBD*DTR
769	SBD=-(((jde-1)/2)*DPHD)
770	SB= SBD*DTR
771	DLM=DLMD*DTR
772	DPH=DPHD*DTR
773	TDLM=DLM+DLM
774	TDPH=DPH+DPH
775	WBI=WB+TDLM
776	SBI=SB+TDPH
777	EBI=WB+(ide-2)*TDLM
778	ANBI=SB+(jde-2)*DPH
779	STPH0=SIN(TPH0)
780	CTPH0=COS(TPH0)
781	TSPH=3600./GRID%DT
782	DO J=JTS,min(JTE,JDE-1)
783	TLM=WB-TDLM+MOD(J,2)*DLM !For velocity points on the E grid
784	TPH=SB+float(J-1)*DPH
785	STPH=SIN(TPH)
786	CTPH=COS(TPH)
787	DO I=ITS,MIN(ITE,IDE-1)
788
789	if (I .eq. ITS) THEN
790	TLM=TLM+TDLM*ITS
791	else
792	TLM=TLM+TDLM
793	endif
794
795	TERM1=(STPH0CTPHCOS(TLM)+CTPH0*STPH)
796	FP=TWOM*(TERM1)
797	F(I,J)=0.5GRID%DTFP
798	ENDDO
799	ENDDO
800	DO J=JTS,min(JTE,JDE-1)
801	TLM=WB-TDLM+MOD(J+1,2)*DLM !For mass points on the E grid
802	TPH=SB+float(J-1)*DPH
803	STPH=SIN(TPH)
804	CTPH=COS(TPH)
805	DO I=ITS,MIN(ITE,IDE-1)
806
807	if (I .eq. ITS) THEN
808	TLM=TLM+TDLM*ITS
809	else
810	TLM=TLM+TDLM
811	endif
812
813	TERM1=(STPH0CTPHCOS(TLM)+CTPH0*STPH)
814	APH=ASIN(TERM1)
815	TG_ALT(I,J)=TG0+TGA*COS(APH)-FIS(I,J)/3333.
816	ENDDO
817	ENDDO
818
819	DO j = jts, MIN(jde-1,jte)
820	DO i = its, MIN(ide-1,ite)
821	! IF ( ( landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) .AND. &
822	! SICE(I,J) .eq. 0. ) THEN
823	! TG(i,j) = sst(i,j)
824	! ELSEIF (SICE(I,J) .eq. 1) THEN
825	! TG(i,j) = 271.16
826	! END IF
827
828	if (TG(I,J) .lt. 200.) then ! only use default TG_ALT definition if
829	! not getting TGROUND from SI
830	TG(I,J)=TG_ALT(I,J)
831	endif
832
833	if (TG(I,J) .lt. 200. .or. TG(I,J) .gt. 320.) then
834	write(message,*) 'problematic TG point at : ', I,J
835	CALL wrf_message( message )
836	endif
837
838	adum2d(i,j)=nmm_tsk(I,J)+sst(I,J)
839
840	END DO
841	END DO
842
843	DEALLOCATE(TG_ALT)
844
845	CALL process_soil_real ( adum2d, TG , &
846	landmask, sst, &
847	st_input, sm_input, sw_input, &
848	st_levels_input , sm_levels_input , &
849	sw_levels_input , &
850	sldpth , dzsoil , stc , smc , sh2o, &
851	flag_sst , flag_soilt000, flag_soilm000, &
852	ids , ide , jds , jde , kds , kde , &
853	ims , ime , jms , jme , kms , kme , &
854	its , ite , jts , jte , kts , kte , &
855	model_config_rec%sf_surface_physics(grid%id) , &
856	model_config_rec%num_soil_layers , &
857	model_config_rec%real_data_init_type , &
858	num_st_levels_input , num_sm_levels_input , &
859	num_sw_levels_input , &
860	num_st_levels_alloc , num_sm_levels_alloc , &
861	num_sw_levels_alloc )
862	write(0,*)' its=',its,' ite=',ite,' jts=',jts,' jte=',jte
863	write(0,*)' ide=',ide,' jde=',jde
864
865	! Minimum soil values, residual, from RUC LSM scheme. For input from Noah and using
866	! RUC LSM scheme, this must be subtracted from the input total soil moisture. For
867	! input RUC data and using the Noah LSM scheme, this value must be added to the soil
868	! moisture input.
869
870	lqmi(1:num_soil_top_cat) = &
871	(/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
872	0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
873	0.004, 0.065 /) !dusan , 0.020, 0.004, 0.008 /)
874
875	! At the initial time we care about values of soil moisture and temperature, other times are
876	! ignored by the model, so we ignore them, too.
877
878	account_for_zero_soil_moisture : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
879
880	CASE ( LSMSCHEME , NMMLSMSCHEME)
881	iicount = 0
882	IF ( FLAG_SM000010 .EQ. 1 ) THEN
883	DO j = jts, MIN(jde-1,jte)
884	DO i = its, MIN(ide-1,ite)
885	IF ( (landmask(i,j).gt.0.5) .and. ( stc(i,1,j) .gt. 200 ) .and. &
886	( stc(i,1,j) .lt. 400 ) .and. ( smc(i,1,j) .lt. 0.005 ) ) then
887	print *,'Noah > Noah: bad soil moisture at i,j = ',i,j,smc(i,:,j)
888	iicount = iicount + 1
889	smc(i,:,j) = 0.005
890	END IF
891	END DO
892	END DO
893	IF ( iicount .GT. 0 ) THEN
894	print *,'Noah -> Noah: total number of small soil moisture locations = ',iicount
895	END IF
896	ELSE IF ( FLAG_SOILM000 .EQ. 1 ) THEN
897	DO j = jts, MIN(jde-1,jte)
898	DO i = its, MIN(ide-1,ite)
899	smc(i,:,j) = smc(i,:,j) + lqmi(isltyp(i,j))
900	END DO
901	END DO
902	DO j = jts, MIN(jde-1,jte)
903	DO i = its, MIN(ide-1,ite)
904	IF ( (landmask(i,j).gt.0.5) .and. ( stc(i,1,j) .gt. 200 ) .and. &
905	( stc(i,1,j) .lt. 400 ) .and. ( smc(i,1,j) .lt. 0.004 ) ) then
906	print *,'RUC -> Noah: bad soil moisture at i,j = ',i,j,smc(i,:,j)
907	iicount = iicount + 1
908	smc(i,:,j) = 0.004
909	END IF
910	END DO
911	END DO
912	IF ( iicount .GT. 0 ) THEN
913	print *,'RUC -> Noah: total number of small soil moisture locations = ',iicount
914	END IF
915	END IF
916	CASE ( RUCLSMSCHEME )
917	iicount = 0
918	IF ( FLAG_SM000010 .EQ. 1 ) THEN
919	DO j = jts, MIN(jde-1,jte)
920	DO i = its, MIN(ide-1,ite)
921	smc(i,:,j) = MAX ( smc(i,:,j) - lqmi(isltyp(i,j)) , 0. )
922	END DO
923	END DO
924	ELSE IF ( FLAG_SOILM000 .EQ. 1 ) THEN
925	! no op
926	END IF
927
928	END SELECT account_for_zero_soil_moisture
929
930	!!! zero out NMM_TSK at water points again
931
932	DO j = jts, MIN(jde-1,jte)
933	DO i = its, MIN(ide-1,ite)
934	if (SM(I,J) .gt. 0.5) then
935	NMM_TSK(I,J)=0.
936	endif
937	END DO
938	END DO
939
940	!! check on STC
941
942	DO j = jts, MIN(jde-1,jte)
943	DO i = its, MIN(ide-1,ite)
944
945	IF (SICE(I,J) .eq. 1.0) then
946	DO L = 1, grid%num_soil_layers
947	STC(I,L,J)=271.16 ! TG value used by Eta/NMM
948	!tgs - initialize sea ice temp profile if needed
949	! mid_point_depth=(3./grid%num_soil_layers)/2. + &
950	! (l-1)*(3./grid%num_soil_layers)
951	! stc(i,l,j) = ( (3.-mid_point_depth)*nmm_tsk(i,j) + &
952	! mid_point_depth*271.16 ) / 3.
953
954	END DO
955	END IF
956
957	END DO
958	END DO
959
960	DO j = jts, MIN(jde-1,jte)
961	DO i = its, MIN(ide-1,ite)
962
963	if (STC(I,1,J) .eq. 0) then
964	write(0,*) 'troublesome STC,SMC value: ', I,J, stc(I,1,J),smc(I,1,J)
965	do JJ=J-1,J+1
966	do L=1, grid%num_soil_layers
967	do II=I-1,I+1
968
969	if (II .ge. its .and. II .le. MIN(ide-1,ite) .and. &
970	JJ .ge. jts .and. JJ .le. MIN(jde-1,jte)) then
971
972	STC(I,L,J)=amax1(STC(I,L,J),STC(II,L,JJ))
973	cur_smc=SMC(I,L,J)
974
975	if ( SMC(II,L,JJ) .gt. 0.005 .and. SMC(II,L,JJ) .lt. 1.0) then
976	aposs_smc=SMC(II,L,JJ)
977
978	if ( cur_smc .eq. 0 ) then
979	cur_smc=aposs_smc
980	SMC(I,L,J)=cur_smc
981	else
982	cur_smc=amin1(cur_smc,aposs_smc)
983	cur_smc=amin1(cur_smc,aposs_smc)
984	SMC(I,L,J)=cur_smc
985	endif
986	endif
987
988	endif ! bounds check
989
990	enddo
991	enddo
992	enddo
993	write(0,*) 'STC, SMC(1) now: ', stc(I,1,J),smc(I,1,J)
994	endif
995
996	if (STC(I,1,J) .eq. 0) then
997	write(0,*) 'STILL troublesome STC value: ', I,J, stc(I,1,J),smc(I,1,J)
998	call wrf_error_fatal("quitting due to significant STC trouble")
999
1000	! do JJ=J-3,J+3
1001	! do L=1, grid%num_soil_layers
1002	! do II=I-3,I+3
1003	! STC(I,L,J)=amax1(STC(I,L,J),STC(II,L,JJ))
1004	! SMC(I,L,J)=amin1(SMC(I,L,J),SMC(II,L,JJ))
1005	! enddo
1006	! enddo
1007	! enddo
1008
1009	endif
1010
1011	ENDDO
1012	ENDDO
1013
1014	!hardwire soil stuff for time being
1015
1016	RTDPTH=0.
1017	RTDPTH(1)=0.1
1018	RTDPTH(2)=0.3
1019	RTDPTH(3)=0.6
1020
1021	! SLDPTH=0.
1022	! SLDPTH(1)=0.1
1023	! SLDPTH(2)=0.3
1024	! SLDPTH(3)=0.6
1025	! SLDPTH(4)=1.0
1026
1027	! write(0,*) 'SLDPTH: ', SLDPTH(1:4)
1028	! write(0,*) 'RTDPTH: ', RTDPTH(1:4)
1029
1030	!!! main body of nmm_specific starts here
1031	!
1032	do J=jts,min(jte,jde-1)
1033	do I=its,min(ite,ide-1)
1034	LMH(I,J)= kme-1 !1
1035	LMV(I,J)= kme-1 !1
1036	RES(I,J)=1.
1037	enddo
1038	enddo
1039
1040	!! HBM2
1041
1042	HBM2=0.
1043
1044	do J=jts,min(jte,jde-1)
1045	do I=its,min(ite,ide-1)
1046
1047	IF ( (J .ge. 3 .and. J .le. (jde-1)-2) .AND. &
1048	(I .ge. 2 .and. I .le. (ide-1)-2+mod(J,2)) ) THEN
1049	HBM2(I,J)=1.
1050	ENDIF
1051	enddo
1052	enddo
1053
1054	!! HBM3
1055	HBM3=0.
1056
1057	!! LOOP OVER LOCAL DIMENSIONS
1058
1059	do J=jts,min(jte,jde-1)
1060	IHWG(J)=mod(J+1,2)-1
1061	IF (J .ge. 4 .and. J .le. (jde-1)-3) THEN
1062	IHL=(ids+1)-IHWG(J)
1063	IHH=(ide-1)-2
1064	do I=its,min(ite,ide-1)
1065	IF (I .ge. IHL .and. I .le. IHH) HBM3(I,J)=1.
1066	enddo
1067	ENDIF
1068	enddo
1069
1070	!! VBM2
1071
1072	VBM2=0.
1073
1074	do J=jts,min(jte,jde-1)
1075	do I=its,min(ite,ide-1)
1076
1077	IF ( (J .ge. 3 .and. J .le. (jde-1)-2) .AND. &
1078	(I .ge. 2 .and. I .le. (ide-1)-1-mod(J,2)) ) THEN
1079
1080	VBM2(I,J)=1.
1081
1082	ENDIF
1083
1084	enddo
1085	enddo
1086
1087	!! VBM3
1088
1089	VBM3=0.
1090
1091	do J=jts,min(jte,jde-1)
1092	do I=its,min(ite,ide-1)
1093
1094	IF ( (J .ge. 4 .and. J .le. (jde-1)-3) .AND. &
1095	(I .ge. 3-mod(J,2) .and. I .le. (ide-1)-2) ) THEN
1096	VBM3(I,J)=1.
1097	ENDIF
1098
1099	enddo
1100	enddo
1101
1102	DTAD=1.0
1103	! IDTCF=DTCF, IDTCF=4
1104	DTCF=4.0 ! used?
1105
1106	DY_NMM=ERAD*DPH
1107	CPGFV=-GRID%DT/(48.*DY_NMM)
1108	EN= GRID%DT/( 4.DY_NMM)DTAD
1109	ENT=GRID%DT/(16.DY_NMM)DTAD
1110
1111	DO J=jts,nnyp
1112	KHL2(J)=(IDE-1)*(J-1)-(J-1)/2+2
1113	KVL2(J)=(IDE-1)*(J-1)-J/2+2
1114	KHH2(J)=(IDE-1)*J-J/2-1
1115	KVH2(J)=(IDE-1)*J-(J+1)/2-1
1116	ENDDO
1117
1118	TPH=SB-DPH
1119
1120	DO J=jts,min(jte,jde-1)
1121	TPH=SB+float(J-1)*DPH
1122	DXP=ERADDLMCOS(TPH)
1123	DXJ(J)=DXP
1124	WPDARJ(J)=-W_NMM * &
1125	((ERADDLMAMIN1(COS(ANBI),COS(SBI)))2+DY_NMM2)/ &
1126	(GRID%DT32.DXP*DY_NMM)
1127
1128	CPGFUJ(J)=-GRID%DT/(48.*DXP)
1129	CURVJ(J)=.5GRID%DTTAN(TPH)/ERAD
1130	FCPJ(J)=GRID%DT/(CP192.DXP*DY_NMM)
1131	FDIVJ(J)=1./(12.DXPDY_NMM)
1132	! EMJ(J)= GRID%DT/( 4.DXP)DTAD
1133	! EMTJ(J)=GRID%DT/(16.DXP)DTAD
1134	FADJ(J)=-GRID%DT/(48.DXPDY_NMM)*DTAD
1135	ACDT=GRID%DTSQRT((ERADDLMAMIN1(COS(ANBI),COS(SBI)))2+DY_NMM*2)
1136	CDDAMP=CODAMP*ACDT
1137	HDACJ(J)=COACACDT/(4.DXP*DY_NMM)
1138	DDMPUJ(J)=CDDAMP/DXP
1139	DDMPVJ(J)=CDDAMP/DY_NMM
1140	ENDDO
1141
1142	!!! wrf_dm_on_monitor block was here, but was causing problems for unknown reasons
1143
1144	DO J=JTS,min(JTE,JDE-1)
1145	TLM=WB-TDLM+MOD(J,2)*DLM
1146	TPH=SB+float(J-1)*DPH
1147	STPH=SIN(TPH)
1148	CTPH=COS(TPH)
1149	DO I=ITS,MIN(ITE,IDE-1)
1150
1151	if (I .eq. ITS) THEN
1152	TLM=TLM+TDLM*ITS
1153	else
1154	TLM=TLM+TDLM
1155	endif
1156
1157	FP=TWOM(CTPH0STPH+STPH0CTPHCOS(TLM))
1158	F(I,J)=0.5GRID%DTFP
1159
1160	ENDDO
1161	ENDDO
1162
1163	! --------------DERIVED VERTICAL GRID CONSTANTS--------------------------
1164
1165	EF4T=.5*GRID%DT/CP
1166	F4Q = -GRID%DT*DTAD
1167	F4D =-.5GRID%DTDTAD
1168
1169	DO L=KDS,KDE-1
1170	RDETA(L)=1./DETA(L)
1171	F4Q2(L)=-.25GRID%DTDTAD/DETA(L)
1172	ENDDO
1173
1174	! shouldnt need to be done globally, right?
1175	DO J=JTS,min(JTE,JDE-1)
1176	DO I=ITS,min(ITE,IDE-1)
1177	DX_NMM(I,J)=DXJ(J)
1178	WPDAR(I,J)=WPDARJ(J)*HBM2(I,J)
1179	CPGFU(I,J)=CPGFUJ(J)*VBM2(I,J)
1180	CURV(I,J)=CURVJ(J)*VBM2(I,J)
1181	FCP(I,J)=FCPJ(J)*HBM2(I,J)
1182	FDIV(I,J)=FDIVJ(J)*HBM2(I,J)
1183	FAD(I,J)=FADJ(J)
1184	HDACV(I,J)=HDACJ(J)*VBM2(I,J)
1185	HDAC(I,J)=HDACJ(J)1.25HBM2(I,J)
1186	ENDDO
1187	ENDDO
1188
1189	! DO J=3,(JDE-1)-2
1190	DO J=JTS, MIN(JDE-1,JTE)
1191
1192	IF (J.LE.5.OR.J.GE.(JDE-1)-4) THEN
1193
1194	KHH=(IDE-1)-2+MOD(J,2) ! KHH is global...loop over I that have
1195	DO I=ITS,MIN(IDE-1,ITE)
1196	IF (I .ge. 2 .and. I .le. KHH) THEN
1197	HDAC(I,J)=HDAC(I,J)* DFC
1198	ENDIF
1199	ENDDO
1200
1201	ELSE
1202
1203	KHH=2+MOD(J,2)
1204	DO I=ITS,MIN(IDE-1,ITE)
1205	IF (I .ge. 2 .and. I .le. KHH) THEN
1206	HDAC(I,J)=HDAC(I,J)* DFC
1207	ENDIF
1208	ENDDO
1209
1210	KHH=(IDE-1)-2+MOD(J,2)
1211
1212	! DO I=(IDE-1)-2,KHH
1213	DO I=ITS,MIN(IDE-1,ITE)
1214	IF (I .ge. (IDE-1)-2 .and. I .le. KHH) THEN
1215	HDAC(I,J)=HDAC(I,J)* DFC
1216	ENDIF
1217	ENDDO
1218	ENDIF
1219	ENDDO
1220
1221	DO J=JTS,min(JTE,JDE-1)
1222	DO I=ITS,min(ITE,IDE-1)
1223	DDMPU(I,J)=DDMPUJ(J)*VBM2(I,J)
1224	DDMPV(I,J)=DDMPVJ(J)*VBM2(I,J)
1225	HDACV(I,J)=HDACV(I,J)*VBM2(I,J)
1226	ENDDO
1227	ENDDO
1228	! --------------INCREASING DIFFUSION ALONG THE BOUNDARIES----------------
1229
1230	! DO J=3,JDE-1-2
1231
1232	DO J=JTS,MIN(JDE-1,JTE)
1233	IF (J.LE.5.OR.J.GE.JDE-1-4) THEN
1234	KVH=(IDE-1)-1-MOD(J,2)
1235	DO I=ITS,min(IDE-1,ITE)
1236	IF (I .ge. 2 .and. I .le. KVH) THEN
1237	DDMPU(I,J)=DDMPU(I,J)*DDFC
1238	DDMPV(I,J)=DDMPV(I,J)*DDFC
1239	HDACV(I,J)=HDACV(I,J)* DFC
1240	ENDIF
1241	ENDDO
1242	ELSE
1243	KVH=3-MOD(J,2)
1244	DO I=ITS,min(IDE-1,ITE)
1245	IF (I .ge. 2 .and. I .le. KVH) THEN
1246	DDMPU(I,J)=DDMPU(I,J)*DDFC
1247	DDMPV(I,J)=DDMPV(I,J)*DDFC
1248	HDACV(I,J)=HDACV(I,J)* DFC
1249	ENDIF
1250	ENDDO
1251	KVH=(IDE-1)-1-MOD(J,2)
1252	DO I=ITS,min(IDE-1,ITE)
1253	IF (I .ge. IDE-1-2 .and. I .le. KVH) THEN
1254	DDMPU(I,J)=DDMPU(I,J)*DDFC
1255	DDMPV(I,J)=DDMPV(I,J)*DDFC
1256	HDACV(I,J)=HDACV(I,J)* DFC
1257	ENDIF
1258	ENDDO
1259	ENDIF
1260	ENDDO
1261
1262	write(0,*) ' grid%num_soil_layers = ', grid%num_soil_layers
1263
1264	write(0,*) 'STC(1)'
1265	do J=min(jde-1,jte),jts,-(jte-jts)/15
1266	write(0,635) (stc(I,1,J),I=its,min(ite,ide-1),(ite-its)/12)
1267	enddo
1268
1269	write(0,*) 'SMC(1)'
1270	do J=min(jde-1,jte),jts,-(jte-jts)/15
1271	write(0,635) (smc(I,1,J),I=its,min(ite,ide-1),(ite-its)/12)
1272	enddo
1273
1274	! write(0,*) 'SM'
1275	! do J=min(jde-1,jte),jts,-(jte-jts)/15
1276	! write(0,635) (smc(I,1,J),I=its,min(ite,ide-1),(ite-its)/12)
1277	! enddo
1278
1279	DO j = jts, MIN(jde-1,jte)
1280	DO i= ITS, MIN(IDE-1,ITE)
1281
1282	if (SM(I,J) .eq. 0 .and. SMC(I,1,J) .gt. 0.5 .and. SICE(I,J) .eq. 0) then
1283	write(0,*) 'wet on land point: ', I,J,SMC(I,1,J),SICE(I,J)
1284	endif
1285
1286	enddo
1287	enddo
1288
1289	!!!! MOVE MONITOR BLOCK HERE
1290
1291	!!! compute EMT, EM on global domain, and only on task 0.
1292
1293	!!! this block also inhibits running as a serial job, it would seem.
1294
1295	IF (wrf_dm_on_monitor()) THEN !!!! NECESSARY TO LIMIT THIS TO TASK ZERO?
1296
1297	ALLOCATE(EMJ(JDS:JDE-1),EMTJ(JDS:JDE-1))
1298
1299	! write(0,*) 'FIGURING OUT EMJ, EMTJ ', JDS, JDE-1
1300	DO J=JDS,JDE-1
1301	TPH=SB+float(J-1)*DPH
1302	DXP=ERADDLMCOS(TPH)
1303	EMJ(J)= GRID%DT/( 4.DXP)DTAD
1304	EMTJ(J)=GRID%DT/(16.DXP)DTAD
1305	! write(0,*) 'J, EMTJ(J): ', J, EMTJ(J)
1306	ENDDO
1307
1308	JA=0
1309	DO 161 J=3,5
1310	JA=JA+1
1311	KHLA(JA)=2
1312	KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
1313	161 EMT(JA)=EMTJ(J)
1314	DO 162 J=(JDE-1)-4,(JDE-1)-2
1315	JA=JA+1
1316	KHLA(JA)=2
1317	KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
1318	162 EMT(JA)=EMTJ(J)
1319	DO 163 J=6,(JDE-1)-5
1320	JA=JA+1
1321	KHLA(JA)=2
1322	KHHA(JA)=2+MOD(J,2)
1323	163 EMT(JA)=EMTJ(J)
1324	DO 164 J=6,(JDE-1)-5
1325	JA=JA+1
1326	KHLA(JA)=(IDE-1)-2
1327	KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
1328	164 EMT(JA)=EMTJ(J)
1329
1330	! --------------SPREADING OF UPSTREAM VELOCITY-POINT ADVECTION FACTOR----
1331
1332	JA=0
1333	DO 171 J=3,5
1334	JA=JA+1
1335	KVLA(JA)=2
1336	KVHA(JA)=(IDE-1)-1-MOD(J,2)
1337	171 EM(JA)=EMJ(J)
1338	DO 172 J=(JDE-1)-4,(JDE-1)-2
1339	JA=JA+1
1340	KVLA(JA)=2
1341	KVHA(JA)=(IDE-1)-1-MOD(J,2)
1342	172 EM(JA)=EMJ(J)
1343	DO 173 J=6,(JDE-1)-5
1344	JA=JA+1
1345	KVLA(JA)=2
1346	KVHA(JA)=2+MOD(J+1,2)
1347	173 EM(JA)=EMJ(J)
1348	DO 174 J=6,(JDE-1)-5
1349	JA=JA+1
1350	KVLA(JA)=(IDE-1)-2
1351	KVHA(JA)=(IDE-1)-1-MOD(J,2)
1352	174 EM(JA)=EMJ(J)
1353
1354	696 continue
1355	ENDIF ! wrf_dm_on_monitor
1356
1357
1358	if(model_config_rec%sf_surface_physics(grid%id).EQ.99) then
1359	call NMM_SH2O(IMS,IME,JMS,JME,ITS,NNXP,JTS,NNYP,grid%num_soil_layers,ISLTYP, &
1360	SM,SICE,STC,SMC,SH2O)
1361	endif
1362
1363	!! must be a better place to put this, but will eliminate "phantom"
1364	!! wind points here (no wind point on eastern boundary of odd numbered rows)
1365
1366	if ( abs(IDE-1-ITE) .eq. 1 ) THEN ! along eastern boundary
1367	write(0,*) 'zero phantom winds'
1368	do K=1,KDE-1
1369	! do J=JTS,JDE-1,2
1370	DO J=JDS,JDE-1,2
1371	if (J .ge. JTS .and. J .le. JTE) THEN
1372	u(IDE-1,K,J)=0.
1373	v(IDE-1,K,J)=0.
1374	endif
1375	enddo
1376	enddo
1377	endif
1378
1379	969 continue
1380
1381	! write(0,*) 'NMM_TSK leaving init_domain_nmm'
1382	! do J=min(jte,jde-1),jts,-(jte-jts)/15
1383	! write(0,635) (NMM_TSK(I,J),I=its,min(ide-1,ite),(ite-its)/12)
1384	! enddo
1385
1386	DO j = jms, jme
1387	DO i = ims, ime
1388
1389	fisx=max(fis(i,j),0.)
1390	Z0(I,J) =SM(I,J)Z0SEA+(1.-SM(I,J)) &
1391	& (Z0(I,J)Z0MAX+FISx FCM+Z0LAND)
1392
1393	ENDDO
1394	ENDDO
1395
1396	write(0,*) 'Z0 over memory, leaving module_initialize_real'
1397	do J=JME,JMS,-(JME-JMS)/20
1398	write(0,635) (Z0(I,J),I=IMS,IME,(IME-IMS)/14)
1399	enddo
1400
1401
1402
1403	write(0,*) 'leaving init_domain_nmm'
1404	!
1405	! Gopal's doing's : following lines moved to namelist_input4 in Registry
1406	!
1407	! write(0,*) 'hardwire'
1408	! sigma=.true.
1409	! grid%IDTAD=2
1410	! grid%NSOIL=4
1411	! grid%NPHS=4
1412	! grid%NCNVC=4
1413	write(message,*)'STUFF MOVED TO REGISTRY:',grid%IDTAD, &
1414	& grid%NSOIL,grid%NRADL,grid%NRADS,grid%NPHS,grid%NCNVC,grid%sigma
1415	CALL wrf_message( TRIM(message) )
1416	!==================================================================================
1417
1418	#define COPY_OUT
1419	#include <nmm_scalar_derefs.inc>
1420	RETURN
1421
1422	END SUBROUTINE init_domain_nmm
1423
1424	!--------------------------------------------------------------------
1425	SUBROUTINE NMM_SH2O(IMS,IME,JMS,JME,ISTART,IM,JSTART,JM,&
1426	NSOIL,ISLTPK, &
1427	SM,SICE,STC,SMC,SH2O)
1428
1429	!! INTEGER, PARAMETER:: NSOTYP=9
1430	! INTEGER, PARAMETER:: NSOTYP=16
1431	INTEGER, PARAMETER:: NSOTYP=19 !!!!!!!!MAYBE???
1432
1433	REAL :: PSIS(NSOTYP),BETA(NSOTYP),SMCMAX(NSOTYP)
1434	REAL :: STC(IMS:IME,NSOIL,JMS:JME), &
1435	SMC(IMS:IME,NSOIL,JMS:JME)
1436	REAL :: SH2O(IMS:IME,NSOIL,JMS:JME),SICE(IMS:IME,JMS:JME),&
1437	SM(IMS:IME,JMS:JME)
1438	REAL :: HLICE,GRAV,T0,BLIM
1439	INTEGER :: ISLTPK(IMS:IME,JMS:JME)
1440
1441	! Constants used in cold start SH2O initialization
1442	DATA HLICE/3.335E5/,GRAV/9.81/,T0/273.15/
1443	DATA BLIM/5.5/
1444	! DATA PSIS /0.04,0.62,0.47,0.14,0.10,0.26,0.14,0.36,0.04/
1445	! DATA BETA /4.26,8.72,11.55,4.74,10.73,8.17,6.77,5.25,4.26/
1446	! DATA SMCMAX /0.421,0.464,0.468,0.434,0.406, &
1447	! 0.465,0.404,0.439,0.421/
1448
1449
1450	!!! NOT SURE...PSIS=SATPSI, BETA=BB??
1451
1452	DATA PSIS /0.069, 0.036, 0.141, 0.759, 0.759, 0.355, &
1453	0.135, 0.617, 0.263, 0.098, 0.324, 0.468, &
1454	0.355, 0.000, 0.069, 0.036, 0.468, 0.069, 0.069 /
1455
1456	DATA BETA/2.79, 4.26, 4.74, 5.33, 5.33, 5.25, &
1457	6.66, 8.72, 8.17, 10.73, 10.39, 11.55, &
1458	5.25, 0.00, 2.79, 4.26, 11.55, 2.79, 2.79 /
1459
1460	DATA SMCMAX/0.339, 0.421, 0.434, 0.476, 0.476, 0.439, &
1461	0.404, 0.464, 0.465, 0.406, 0.468, 0.468, &
1462	0.439, 1.000, 0.200, 0.421, 0.468, 0.200, 0.339/
1463
1464	!tgs - the following table values are from REDPRM_USGS in NMM LSM
1465	! DATA PSIS /0.0350, 0.0363, 0.1413, 0.7586, 0.7586, 0.3548, &
1466	! 0.1349, 0.6166, 0.2630, 0.0977, 0.3236, 0.4677, &
1467	! 0.3548, 0.0000, 0.0350, 0.0363, 0.4677, 0.0350, &
1468	! 0.0350 /
1469
1470	! DATA BETA /4.05, 4.26, 4.74, 5.33, 5.33, 5.25, &
1471	! 6.77, 8.72, 8.17, 10.73, 10.39, 11.55, &
1472	! 5.25, 0.00, 4.05, 4.26, 11.55, 4.05, &
1473	! 4.05 /
1474
1475	! DATA SMCMAX/0.395, 0.421, 0.434, 0.476, 0.476, 0.439, &
1476	! 0.404, 0.464, 0.465, 0.406, 0.468, 0.457, &
1477	! 0.464, 0.000, 0.200, 0.421, 0.457, 0.200, &
1478	! 0.395 /
1479
1480
1481	write(0,*) 'define SH2O over IM,JM: ', IM,JM
1482	DO K=1,NSOIL
1483	DO J=JSTART,JM
1484	DO I=ISTART,IM
1485	if(i==169.and.j==111.and.k==1)then
1486	write(0,*)' enter NMM_SH2O k=',k,' smc=',smc(i,k,j),' sh2o=',sh2o(i,k,j)
1487	endif
1488	!tst
1489	IF (SMC(I,K,J) .gt. SMCMAX(ISLTPK(I,J))) then
1490	! if (K .eq. 1) then
1491	! write(0,*) 'I,J,reducing SMC from ' ,I,J,SMC(I,K,J), 'to ', SMCMAX(ISLTPK(I,J))
1492	! endif
1493	SMC(I,K,J)=SMCMAX(ISLTPK(I,J))
1494	ENDIF
1495	!tst
1496
1497	if(i==056.and.j==265.and.k==1)then
1498	write(0,*)' NMM_SH2O point 2 k=',k,' smc=',smc(i,k,j),' sh2o=',sh2o(i,k,j)
1499	endif
1500	IF ( (SM(I,J) .lt. 0.5) .and. (SICE(I,J) .lt. 0.5) ) THEN
1501
1502	IF (ISLTPK(I,J) .gt. 19) THEN
1503	WRITE(6,*) 'FORCING ISLTPK at : ', I,J
1504	ISLTPK(I,J)=9
1505	ELSEIF (ISLTPK(I,J) .le. 0) then
1506	WRITE(6,*) 'FORCING ISLTPK at : ', I,J
1507	ISLTPK(I,J)=1
1508	ENDIF
1509
1510
1511	! cold start: determine liquid soil water content (SH2O)
1512	! SH2O <= SMC for T < 273.149K (-0.001C)
1513
1514	IF (STC(I,K,J) .LT. 273.149) THEN
1515
1516	! first guess following explicit solution for Flerchinger Eqn from Koren
1517	! et al, JGR, 1999, Eqn 17 (KCOUNT=0 in FUNCTION FRH2O).
1518
1519	BX = BETA(ISLTPK(I,J))
1520	IF ( BETA(ISLTPK(I,J)) .GT. BLIM ) BX = BLIM
1521
1522	if ( GRAV*(-PSIS(ISLTPK(I,J))) .eq. 0 ) then
1523	write(0,*) 'TROUBLE'
1524	write(0,*) 'I,J: ', i,J
1525	write(0,*) 'grav, isltpk, psis(isltpk): ', grav,isltpk(I,J),&
1526	psis(isltpk(I,J))
1527	endif
1528
1529	if (BX .eq. 0 .or. STC(I,K,J) .eq. 0) then
1530	write(0,*) 'I,J,BX, STC: ', I,J,BX,STC(I,K,J)
1531	endif
1532	FK = (((HLICE/(GRAV(-PSIS(ISLTPK(I,J))))) &
1533	((STC(I,K,J)-T0)/STC(I,K,J)))** &
1534	(-1/BX))*SMCMAX(ISLTPK(I,J))
1535	IF (FK .LT. 0.02) FK = 0.02
1536	SH2O(I,K,J) = MIN ( FK, SMC(I,K,J) )
1537	! ----------------------------------------------------------------------
1538	! now use iterative solution for liquid soil water content using
1539	! FUNCTION FRH2O (from the Eta "NOAH" land-surface model) with the
1540	! initial guess for SH2O from above explicit first guess.
1541
1542	SH2O(I,K,J)=FRH2O_init(STC(I,K,J),SMC(I,K,J),SH2O(I,K,J), &
1543	SMCMAX(ISLTPK(I,J)),BETA(ISLTPK(I,J)), &
1544	PSIS(ISLTPK(I,J)))
1545
1546	ELSE ! above freezing
1547	SH2O(I,K,J)=SMC(I,K,J)
1548	ENDIF
1549
1550
1551	ELSE ! water point
1552	SH2O(I,K,J)=SMC(I,K,J)
1553
1554	ENDIF ! test on land/ice/sea
1555	if (SH2O(I,K,J) .gt. SMCMAX(ISLTPK(I,J))) then
1556	write(0,*) 'SH2O > THAN SMCMAX ', I,J,SH2O(I,K,J),SMCMAX(ISLTPK(I,J)),SMC(I,K,J)
1557	endif
1558
1559	if(i==169.and.j==111.and.k==1)then
1560	write(0,*)' exit NMM_SH2O k=',k,' smc=',smc(i,k,j),' sh2o=',sh2o(i,k,j)
1561	endif
1562
1563	ENDDO
1564	ENDDO
1565	ENDDO
1566
1567	END SUBROUTINE NMM_SH2O
1568
1569	!-------------------------------------------------------------------
1570
1571	FUNCTION FRH2O_init(TKELV,SMC,SH2O,SMCMAX,B,PSIS)
1572
1573	IMPLICIT NONE
1574
1575	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1576	! PURPOSE: CALCULATE AMOUNT OF SUPERCOOLED LIQUID SOIL WATER CONTENT
1577	! IF TEMPERATURE IS BELOW 273.15K (T0). REQUIRES NEWTON-TYPE ITERATION
1578	! TO SOLVE THE NONLINEAR IMPLICIT EQUATION GIVEN IN EQN 17 OF
1579	! KOREN ET AL. (1999, JGR, VOL 104(D16), 19569-19585).
1580	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1581	!
1582	! New version (JUNE 2001): much faster and more accurate newton iteration
1583	! achieved by first taking log of eqn cited above -- less than 4
1584	! (typically 1 or 2) iterations achieves convergence. Also, explicit
1585	! 1-step solution option for special case of parameter Ck=0, which reduces
1586	! the original implicit equation to a simpler explicit form, known as the
1587	! ""Flerchinger Eqn". Improved handling of solution in the limit of
1588	! freezing point temperature T0.
1589	!
1590	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1591	!
1592	! INPUT:
1593	!
1594	! TKELV.........Temperature (Kelvin)
1595	! SMC...........Total soil moisture content (volumetric)
1596	! SH2O..........Liquid soil moisture content (volumetric)
1597	! SMCMAX........Saturation soil moisture content (from REDPRM)
1598	! B.............Soil type "B" parameter (from REDPRM)
1599	! PSIS..........Saturated soil matric potential (from REDPRM)
1600	!
1601	! OUTPUT:
1602	! FRH2O.........supercooled liquid water content.
1603	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1604
1605	REAL B
1606	REAL BLIM
1607	REAL BX
1608	REAL CK
1609	REAL DENOM
1610	REAL DF
1611	REAL DH2O
1612	REAL DICE
1613	REAL DSWL
1614	REAL ERROR
1615	REAL FK
1616	REAL FRH2O_init
1617	REAL GS
1618	REAL HLICE
1619	REAL PSIS
1620	REAL SH2O
1621	REAL SMC
1622	REAL SMCMAX
1623	REAL SWL
1624	REAL SWLK
1625	REAL TKELV
1626	REAL T0
1627
1628	INTEGER NLOG
1629	INTEGER KCOUNT
1630	PARAMETER (CK=8.0)
1631	! PARAMETER (CK=0.0)
1632	PARAMETER (BLIM=5.5)
1633	! PARAMETER (BLIM=7.0)
1634	PARAMETER (ERROR=0.005)
1635
1636	PARAMETER (HLICE=3.335E5)
1637	PARAMETER (GS = 9.81)
1638	PARAMETER (DICE=920.0)
1639	PARAMETER (DH2O=1000.0)
1640	PARAMETER (T0=273.15)
1641
1642	! ### LIMITS ON PARAMETER B: B < 5.5 (use parameter BLIM) ####
1643	! ### SIMULATIONS SHOWED IF B > 5.5 UNFROZEN WATER CONTENT ####
1644	! ### IS NON-REALISTICALLY HIGH AT VERY LOW TEMPERATURES ####
1645	! ################################################################
1646	!
1647	BX = B
1648	IF ( B .GT. BLIM ) BX = BLIM
1649	! ------------------------------------------------------------------
1650
1651	! INITIALIZING ITERATIONS COUNTER AND ITERATIVE SOLUTION FLAG.
1652	NLOG=0
1653	KCOUNT=0
1654
1655	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1656	! C IF TEMPERATURE NOT SIGNIFICANTLY BELOW FREEZING (T0), SH2O = SMC
1657	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1658
1659
1660	IF (TKELV .GT. (T0 - 1.E-3)) THEN
1661
1662	FRH2O_init=SMC
1663
1664	ELSE
1665
1666	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1667	IF (CK .NE. 0.0) THEN
1668
1669	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1670	! CCCCCCCCC OPTION 1: ITERATED SOLUTION FOR NONZERO CK CCCCCCCCCCC
1671	! CCCCCCCCCCCC IN KOREN ET AL, JGR, 1999, EQN 17 CCCCCCCCCCCCCCCCC
1672
1673	! INITIAL GUESS FOR SWL (frozen content)
1674	SWL = SMC-SH2O
1675	! KEEP WITHIN BOUNDS.
1676	IF (SWL .GT. (SMC-0.02)) SWL=SMC-0.02
1677	IF(SWL .LT. 0.) SWL=0.
1678	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1679	! C START OF ITERATIONS
1680	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1681	DO WHILE (NLOG .LT. 10 .AND. KCOUNT .EQ. 0)
1682	NLOG = NLOG+1
1683	DF = ALOG(( PSISGS/HLICE ) ( ( 1.+CKSWL )2. ) &
1684	( SMCMAX/(SMC-SWL) )**BX) - ALOG(-(TKELV-T0)/TKELV)
1685	DENOM = 2. * CK / ( 1.+CK*SWL ) + BX / ( SMC - SWL )
1686	SWLK = SWL - DF/DENOM
1687	! BOUNDS USEFUL FOR MATHEMATICAL SOLUTION.
1688	IF (SWLK .GT. (SMC-0.02)) SWLK = SMC - 0.02
1689	IF(SWLK .LT. 0.) SWLK = 0.
1690	! MATHEMATICAL SOLUTION BOUNDS APPLIED.
1691	DSWL=ABS(SWLK-SWL)
1692	SWL=SWLK
1693	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1694	! CC IF MORE THAN 10 ITERATIONS, USE EXPLICIT METHOD (CK=0 APPROX.)
1695	! CC WHEN DSWL LESS OR EQ. ERROR, NO MORE ITERATIONS REQUIRED.
1696	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1697	IF ( DSWL .LE. ERROR ) THEN
1698	KCOUNT=KCOUNT+1
1699	END IF
1700	END DO
1701	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1702	! C END OF ITERATIONS
1703	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1704	! BOUNDS APPLIED WITHIN DO-BLOCK ARE VALID FOR PHYSICAL SOLUTION.
1705	FRH2O_init = SMC - SWL
1706
1707	! CCCCCCCCCCCCCCCCCCCCCCCC END OPTION 1 CCCCCCCCCCCCCCCCCCCCCCCCCCC
1708
1709	ENDIF
1710
1711	IF (KCOUNT .EQ. 0) THEN
1712	! Print*,'Flerchinger used in NEW version. Iterations=',NLOG
1713
1714	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1715	! CCCCC OPTION 2: EXPLICIT SOLUTION FOR FLERCHINGER EQ. i.e. CK=0 CCCCCCCC
1716	! CCCCCCCCCCCCC IN KOREN ET AL., JGR, 1999, EQN 17 CCCCCCCCCCCCCCC
1717
1718	FK=(((HLICE/(GS(-PSIS)))((TKELV-T0)/TKELV))*(-1/BX))SMCMAX
1719	! APPLY PHYSICAL BOUNDS TO FLERCHINGER SOLUTION
1720	IF (FK .LT. 0.02) FK = 0.02
1721	FRH2O_init = MIN ( FK, SMC )
1722
1723	! CCCCCCCCCCCCCCCCCCCCCCCCC END OPTION 2 CCCCCCCCCCCCCCCCCCCCCCCCCC
1724
1725	ENDIF
1726
1727	ENDIF
1728
1729	RETURN
1730
1731	END FUNCTION FRH2O_init
1732
1733
1734	!--------------------------------------------------------------------
1735
1736	SUBROUTINE init_module_initialize
1737	END SUBROUTINE init_module_initialize
1738
1739	!---------------------------------------------------------------------
1740
1741	END MODULE module_initialize

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

Download in other formats:

Original Format