UPP  001
 All Data Structures Files Functions Pages
SLP_new.f
1  SUBROUTINE memslp(TPRES,QPRES,FIPRES)
2 !
3 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
4 ! . . .
5 ! SUBROUTINE: MEMSLP MEMBRANE SLP REDUCTION
6 !
7 ! ABSTRACT: THIS ROUTINE COMPUTES THE SEA LEVEL PRESSURE
8 ! REDUCTION USING THE MESINGER RELAXATION
9 ! METHOD FOR SIGMA COORDINATES.
10 ! A BY-PRODUCT IS THE
11 ! SET OF VALUES FOR THE UNDERGROUND TEMPERATURES
12 ! ON THE SPECIFIED PRESSURE LEVELS
13 !
14 ! PROGRAM HISTORY LOG:
15 ! 99-09-23 T BLACK - REWRITTEN FROM ROUTINE SLP (ETA
16 ! COORDINATES)
17 ! 02-07-26 H CHUANG - PARALLIZE AND MODIFIED FOR WRF A/C GRIDS
18 ! ALSO REDUCE S.O.R. COEFF FROM 1.75 to 1.25
19 ! BECAUSE THERE WAS NUMERICAL INSTABILITY
20 ! 02-08-21 H CHUANG - MODIFIED TO ALWAYS USE OLD TTV FOR RELAXATION
21 ! SO THAT THERE WAS BIT REPRODUCIBILITY BETWEEN
22 ! USING ONE AND MULTIPLE TASKS
23 ! 11-04-29 H CHUANG - FIX GFS GIBSING BY USING LM-1 STATE VARIABLES
24 ! TO DERIVE SLP HYDROSTATICALLY
25 ! 13-12-06 H CHUANG - REMOVE EXTRA SMOOTHING OF SLP ITSELF
26 ! CHANGES TO AVOID RELAXATION FOR ABOVE G GIBSING
27 ! ARE COMMENTED OUT FOR NOW
28 ! 19-10-30 Bo CUI - REMOVE "GOTO" STATEMENT
29 ! 21-07-26 W Meng - Restrict computation from undefined grids
30 ! 21-09-25 W Meng - Further modification for restricting computation
31 ! from undefined grids.
32 !
33 ! USAGE: CALL SLPSIG FROM SUBROUITNE ETA2P
34 !
35 ! INPUT ARGUMENT LIST:
36 ! PD - SFC PRESSURE MINUS PTOP
37 ! FIS - SURFACE GEOPOTENTIAL
38 ! T - TEMPERATURE
39 ! Q - SPECIFIC HUMIDITY
40 ! FI - GEOPOTENTIAL
41 ! PT - TOP PRESSURE OF DOMAIN
42 !
43 ! OUTPUT ARGUMENT LIST:
44 ! PSLP - THE FINAL REDUCED SEA LEVEL PRESSURE ARRAY
45 !
46 ! SUBPROGRAMS CALLED:
47 ! UNIQUE:
48 ! NONE
49 !
50 !-----------------------------------------------------------------------
51  use vrbls3d, only: pint, zint, t, q
52  use vrbls2d, only: pslp, fis
53  use masks, only: lmh
54  use params_mod, only: overrc, ad05, cft0, g, rd, d608, h1, kslpd
55  use ctlblk_mod, only: jend, jsta, spval, spl, num_procs, mpi_comm_comp, lsmp1, &
56  jsta_m, jend_m, lm, im, jsta_2l, jend_2u, lsm, jm,&
57  im_jm
58 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59  implicit none
60 !
61  include "mpif.h"
62 !-----------------------------------------------------------------------
63  integer,PARAMETER :: nfill=0,nrlx1=500,nrlx2=100
64  real,parameter:: def_of_mountain=2.0
65 !-----------------------------------------------------------------------
66  real,dimension(IM,JSTA_2L:JEND_2U,LSM),intent(in) :: qpres
67  real,dimension(IM,JSTA_2L:JEND_2U,LSM),intent(inout) :: tpres,fipres
68  REAL :: ttv(im,jsta_2l:jend_2u),tnew(im,jsta_2l:jend_2u) &
69  , p1(im,jsta_2l:jend_2u),htm2d(im,jsta_2l:jend_2u)
70  REAL :: htmo(im,jsta_2l:jend_2u,lsm)
71  real :: p2,tlyr,gz1,gz2,spll,psfc,pchk,slope,tvrtc,dis,tvrt,tem
72 !-----------------------------------------------------------------------
73 !-----------------------------------------------------------------------
74  INTEGER :: kmntm(lsm),imnt(im_jm,lsm),jmnt(im_jm,lsm) &
75  , lmho(im,jsta_2l:jend_2u)
76  INTEGER :: ihe(jm),ihw(jm),ive(jm),ivw(jm),ihs(jm),ihn(jm)
77  integer ii,jj,i,j,l,n,llmh,km,ks,ihh2,kount,kmn,nrlx,lhmnt, &
78  lmhij,lmap1,kmm,lp,lxxx,ierr
79 ! dong
80  real a1,a2,a3,a4,a5,a6,a7,a8
81 !-----------------------------------------------------------------------
82  LOGICAL :: done(im,jsta_2l:jend_2u)
83 !-----------------------------------------------------------------------
84 !***
85 !*** CALCULATE THE I-INDEX EAST-WEST INCREMENTS
86 !***
87 !
88  ii = im/2
89  jj = (jend-jsta)/2
90  DO j=1,jm
91  ihe(j) = 1
92  ihw(j) = -1
93  ihs(j) = -1
94  ihn(j) = 1
95  ive(j) = mod(j,2)
96  ivw(j) = ive(j)-1
97  ENDDO
98 ! print*,'relaxation coeff= ',OVERRC
99 !-----------------------------------------------------------------------
100 !***
101 !*** INITIALIZE ARRAYS. LOAD SLP ARRAY WITH SURFACE PRESSURE.
102 !***
103 !$omp parallel do private(i,j,llmh)
104  DO j=jsta,jend
105  DO i=1,im
106  llmh = nint(lmh(i,j))
107  pslp(i,j) = pint(i,j,llmh+1)
108 ! dong
109 ! TTV(I,J) = 0.
110  ttv(i,j) = spval
111  tnew(i,j) = spval
112 
113 
114  lmho(i,j) = lsm
115  done(i,j) = .false.
116  ENDDO
117  ENDDO
118 !
119 !--------------------------------------------------------------------
120 !***
121 !*** CREATE A 3-D "HEIGHT MASK" FOR THE SPECIFIED PRESSURE LEVELS
122 !*** (1 => ABOVE GROUND) AND A 2-D INDICATOR ARRAY THAT SAYS
123 !*** WHICH PRESSURE LEVEL IS THE LOWEST ONE ABOVE THE GROUND
124 !***
125  DO l=1,lsm
126  spll = spl(l)
127 !
128 !$omp parallel do private(j,i,psfc,pchk)
129  DO j=jsta,jend
130  DO i=1,im
131 
132  htmo(i,j,l)=1.
133  if(pslp(i,j)<spval) then
134 
135  psfc = pslp(i,j)
136  pchk = psfc
137  IF(nfill > 0) THEN
138  pchk = pint(i,j,nint(lmh(i,j))+1-nfill)
139  ENDIF
140  IF(fis(i,j) < 1.) pchk = psfc
141 !
142  IF(spll < pchk) THEN
143  htmo(i,j,l) = 1.
144  ELSE
145  htmo(i,j,l) = 0.
146  IF(l > 1 .AND. htmo(i,j,l-1) > 0.5) lmho(i,j) = l-1
147  ENDIF
148  IF(l == lsm .AND. htmo(i,j,l) > 0.5) lmho(i,j) = lsm
149 !
150 ! test new idea of filtering above-ground pressure levels for Gibsing
151 ! IF(L==LSM.AND.HTMO(I,J,L)>0.5)THEN
152 ! IF(FIS(I,J)>0.)THEN
153 ! LMHO(I,J)=LSM
154 ! ELSE
155 ! LMHO(I,J)=LSM-2
156 ! HTMO(I,J,LSM)=0.
157 ! HTMO(I,J,LSM-1)=0.
158 ! END IF
159 ! END IF
160 ! if(i==ii.and.j==jj)print*,'Debug: HTMO= ',HTMO(I,J,L)
161 
162  endif !if pslp
163 
164  ENDDO
165  ENDDO
166 !
167  ENDDO
168 ! if(jj>=jsta.and.jj<=jend) print*,'Debug: LMHO=',LMHO(ii,jj)
169 !--------------------------------------------------------------------
170 !***
171 !*** WE REACH THIS LINE IF WE WANT THE MESINGER ETA SLP REDUCTION
172 !*** BASED ON RELAXATION TEMPERATURES. THE FIRST STEP IS TO
173 !*** FIND THE HIGHEST LAYER CONTAINING MOUNTAINS.
174 !***
175  lhmnt = lsm
176  loop210: DO l=lsm,1,-1
177 
178  DO j=jsta,jend
179  DO i=1,im
180  if(pslp(i,j)<spval) then
181  IF(htmo(i,j,l) < 0.5) cycle loop210
182  endif
183  ENDDO
184  ENDDO
185  lhmnt = l+1
186  EXIT loop210
187  ENDDO loop210
188  !210 continue
189  !220 continue
190 
191 ! print*,'Debug in SLP: LHMNT=',LHMNT
192 
193  if ( num_procs > 1 ) then
194  CALL mpi_allreduce &
195  (lhmnt,lxxx,1,mpi_integer,mpi_min,mpi_comm_comp,ierr)
196  lhmnt = lxxx
197  end if
198 
199  IF(lhmnt == lsmp1) go to 325
200 
201 ! print*,'Debug in SLP: LHMNT A ALLREDUCE=',LHMNT
202 !***
203 !*** NOW GATHER THE ADDRESSES OF ALL THE UNDERGROUND POINTS.
204 !***
205 !!$omp parallel do private(kmn,kount)
206  DO l=lhmnt,lsm
207  kmn = 0
208  kmntm(l) = 0
209  kount = 0
210 ! DO 240 J=JSTA_M2,JEND_M2
211  DO j=jsta_m,jend_m
212  DO i=2,im-1
213  if(pslp(i,j)<spval) then
214  kount = kount + 1
215  imnt(kount,l) = 0
216  jmnt(kount,l) = 0
217  IF(htmo(i,j,l) > 0.5) cycle
218  kmn = kmn + 1
219  imnt(kmn,l) = i
220  jmnt(kmn,l) = j
221  endif
222  enddo
223  enddo
224  kmntm(l) = kmn
225  enddo
226 !
227 !
228 !*** CREATE A TEMPORARY TV ARRAY, AND FOLLOW BY SEQUENTIAL
229 !*** OVERRELAXATION, DOING NRLX PASSES.
230 !
231 ! IF(NTSD==1)THEN
232  nrlx = nrlx2
233 ! ELSE
234 ! NRLX=NRLX2
235 ! ENDIF
236 !
237 !!$omp parallel do private(i,j,ttv,tem,kmma (Can this loop be threaded?))
238  DO l=lhmnt,lsm
239 !
240 !$omp parallel do private(i,j)
241  DO j=jsta,jend
242  DO i=1,im
243 ! dong
244 ! if (QPRES(I,J,LSM) < spval) then
245  !if(PSLP(I,J)<spval) then
246  ttv(i,j) = tpres(i,j,l)
247  htm2d(i,j) = htmo(i,j,l)
248  !end if ! spval if
249 ! end if ! spval if
250 ! IF(TTV(I,J)<150. .and. TTV(I,J)>325.0)print* &
251 ! ,'abnormal IC for T relaxation',i,j,TTV(I,J)
252  enddo
253  enddo
254 !
255 !*** FOR GRID BOXES NEXT TO MOUNTAINS, COMPUTE TV TO USE AS
256 !*** BOUNDARY CONDITIONS FOR THE RELAXATION UNDERGROUND
257 !
258  CALL exch(htm2d(1,jsta_2l)) !ONLY NEED TO EXCHANGE ONE ROW FOR A/C GRID
259 ! DO J=JSTA_M2,JEND_M2
260 !$omp parallel do private(i,j,tem)
261  DO j=jsta_m,jend_m
262  DO i=2,im-1
263 
264  if(pslp(i,j)<spval) then
265 
266 !HC IF(HTM2D(I,J,L)>0.5.AND.
267 !HC 1 HTM2D(I+IHW(J),J-1,L)*HTM2D(I+IHE(J),J-1,L)
268 !HC 2 *HTM2D(I+IHW(J),J+1,L)*HTM2D(I+IHE(J),J+1,L)
269 !HC 3 *HTM2D(I-1 ,J ,L)*HTM2D(I+1 ,J ,L)
270 !HC 4 *HTM2D(I ,J-2,L)*HTM2D(I ,J+2,L)<0.5)THEN
271 !HC MODIFICATION FOR C AND A GRIDS
272 
273  tem = htm2d(i-1,j)*htm2d(i+1,j)*htm2d(i,j-1)*htm2d(i,j+1) &
274  * htm2d(i-1,j-1)*htm2d(i+1,j-1)*htm2d(i-1,j+1)*htm2d(i+1,j+1)
275  IF(htm2d(i,j) > 0.5 .AND. tem < 0.5) then
276  ttv(i,j) = tpres(i,j,l)*(1.+0.608*qpres(i,j,l))
277  ENDIF
278 ! if(i==ii.and.j==jj)print*,'Debug:L,TTV B SMOO= ',l,TTV(I,J)
279  end if ! spval
280  ENDDO
281  ENDDO
282 !
283  kmm = kmntm(l)
284 
285 ! print*,'Debug:L,KMM=',L,KMM
286 !
287  DO n=1,nrlx
288  CALL exch(ttv(1,jsta_2l))
289 !!$omp parallel do private(i,j,km,a1,a2,a3,a4,a5,a6,a7,a8)
290  DO km=1,kmm
291  i = imnt(km,l)
292  j = jmnt(km,l)
293 
294  if(pslp(i,j)<spval) then
295 
296 !HC TTV(I,J)=AD05*(4.*(TTV(I+IHW(J),J-1)+TTV(I+IHE(J),J-1)
297 !HC 1 +TTV(I+IHW(J),J+1)+TTV(I+IHE(J),J+1))
298 !HC 2 +TTV(I-1,J) +TTV(I+1,J)
299 !HC 3 +TTV(I,J-2) +TTV(I,J+2))
300 !HC 4 -CFT0*TTV(I,J)
301 !HC MODIFICATION FOR C AND A GRIDS
302 ! eight point relaxation using updated TTV to the lower and left
303 ! TTV(I,J)=AD05*(4.*(TTV(I-1,J)+TTV(I+1,J)
304 ! 1 +TTV(I,J-1)+TTV(I,J+1))
305 ! 2 +TTV(I-1,J-1)+TTV(I+1,J-1)
306 ! 3 +TTV(I-1,J+1)+TTV(I+1,J+1))
307 ! 4 -CFT0*TTV(I,J)
308 ! eight point relaxation using old TTV
309  a1=ttv(i-1,j)
310  a2=ttv(i+1,j)
311  a3=ttv(i,j-1)
312  a4=ttv(i,j+1)
313  a5=ttv(i-1,j-1)
314  a6=ttv(i+1,j-1)
315  a7=ttv(i-1,j+1)
316  a8=ttv(i+1,j+1)
317 ! if ((a1-spval) <= 1e-10) a1=TTV(I,J)
318 ! if ((a2-spval) <= 1e-10) a2=TTV(I,J)
319 ! if ((a3-spval) <= 1e-10) a3=TTV(I,J)
320 ! if ((a4-spval) <= 1e-10) a4=TTV(I,J)
321 ! if ((a5-spval) <= 1e-10) a5=TTV(I,J)
322 ! if ((a6-spval) <= 1e-10) a6=TTV(I,J)
323 ! if ((a7-spval) <= 1e-10) a7=TTV(I,J)
324 ! if ((a8-spval) <= 1e-10) a8=TTV(I,J)
325 
326  if ((a1 < spval) .and. &
327  (a2 < spval) .and. &
328  (a3 < spval) .and. &
329  (a4 < spval) .and. &
330  (a5 < spval) .and. &
331  (a6 < spval) .and. &
332  (a7 < spval) .and. &
333  (a8 < spval) .and. (ttv(i,j) < spval)) then
334 
335 ! TNEW(I,J) = AD05*(4.*(a1 +a2 +a3 &
336 ! +a4) +a5 +a6 &
337 ! +a7+a8)-TTV(I,J)*CFT0
338 
339  tnew(i,j) = ad05*(4.*(ttv(i-1,j) +ttv(i+1,j) +ttv(i,j-1) &
340  +ttv(i,j+1)) +ttv(i-1,j-1) +ttv(i+1,j-1) &
341  +ttv(i-1,j+1)+ttv(i+1,j+1))-ttv(i,j)*cft0
342  else
343  tnew(i,j) = ttv(i,j)
344  end if ! spval
345 
346 ! four point relaxation using old TTV
347 ! TNEW(I,J)=TTV(I,J)+1.0*((TTV(I-1,J)+TTV(I+1,J)
348 ! 1 +TTV(I,J-1)+TTV(I,J+1)-4.0*TTV(I,J))/4.0)
349 ! four point relaxation using updated TTV to the lower and left
350 ! TTV(I,J)=TTV(I,J)+1.0*((TTV(I-1,J)+TTV(I+1,J)
351 ! 1 +TTV(I,J-1)+TTV(I,J+1)-4.0*TTV(I,J))/4.0)
352 !
353 ! if(i==ii.and.j==jj)print*,'Debug: L,TTV A S'
354 ! 1,l,TTV(I,J),N
355 ! 1,l,TNEW(I,J),N
356 
357  end if ! spval
358 
359  enddo
360 !
361 !!$omp parallel do private(i,j,km)
362  DO km=1,kmm
363  i = imnt(km,l)
364  j = jmnt(km,l)
365  if(pslp(i,j)<spval .and. tnew(i,j)< spval/100.) then
366  ttv(i,j) = tnew(i,j)
367  end if ! spval
368  END DO
369  END DO ! NRLX loop
370 !
371 !!$omp parallel do private(i,j,km)
372  DO km=1,kmm
373  i = imnt(km,l)
374  j = jmnt(km,l)
375 
376  if(pslp(i,j)<spval) then
377 
378 ! dong try to fix missing value for hgtprs at 1000 mb
379  tpres(i,j,l) = ttv(i,j)
380  end if ! spval
381 
382 ! if (QPRES(I,J,L) < 1000) TPRES(I,J,L) = TTV(I,J)
383 ! if (QPRES(I,J,L) < 1000) TPRES(I,J,L) = 1
384 
385  END DO
386  enddo ! end of l loop
387 !----------------------------------------------------------------
388 !***
389 !*** CALCULATE THE SEA LEVEL PRESSURE AS PER THE NEW SCHEME.
390 !*** INTEGRATE THE HYDROSTATIC EQUATION DOWNWARD FROM THE
391 !*** GROUND THROUGH EACH OUTPUT PRESSURE LEVEL (WHERE TV
392 !*** IS NOW KNOWN) TO FIND GZ AT THE NEXT MIDPOINT BETWEEN
393 !*** PRESSURE LEVELS. WHEN GZ=0 IS REACHED, SOLVE FOR THE
394 !*** PRESSURE.
395 !***
396 !
397 !*** BEFORE APPLYING RELAXATION FOR UNDERGROUND POINTS,
398 !*** FIRST FIND GRID POINTS AT/NEAR/BELOW SEA LEVEL AND DERIVE
399 !*** SEA LEVEL PRESSURE TO AVOID MEMBRANE RELAXATION
400 !*** AT THESE GRID POINTS. E.G. HURRICANE CENTER NEAR COAST
401 !
402  kount = 0
403  DO j=jsta,jend
404  DO i=1,im
405 
406  if(pslp(i,j)<spval) then
407 ! P1(I,J)=SPL(NINT(LMH(I,J)))
408 ! DONE(I,J)=.FALSE.
409 
410  IF(abs(fis(i,j)) < 1.) THEN
411  pslp(i,j) = pint(i,j,nint(lmh(i,j))+1)
412  done(i,j) = .true.
413  kount = kount + 1
414 ! if(i==ii.and.j==jj)print*,'Debug:DONE,PSLP A S1=' &
415 ! ,done(i,j),PSLP(I,J)
416  ELSE IF(fis(i,j) < -1.0) THEN
417  DO l=lm,1,-1
418  IF(zint(i,j,l) > 0.)THEN
419 ! PSLP(I,J)=PINT(I,J,L)/EXP(-ZINT(I,J,L)*G &
420 ! /(RD*T(I,J,L)*(Q(I,J,L)*D608+1.0)))
421 
422  tem = 0.5*(t(i,j,l)+t(i,j,l-1))*(1.0+0.5*d608*(q(i,j,l)+q(i,j,l-1)))
423  pslp(i,j) = pint(i,j,l-1)/exp(-zint(i,j,l-1)*g/(rd*tem))
424  done(i,j) = .true.
425 ! if(i==ii.and.j==jj)print* &
426 ! ,'Debug:DONE,PINT,PSLP A S1=' &
427 ! ,done(i,j),PINT(I,J,L),PSLP(I,J)
428  exit
429  END IF
430  END DO
431  ENDIF
432 
433  end if ! spval
434 
435  ENDDO
436  ENDDO
437 !
438  kmm = kmntm(lsm)
439 !!$omp parallel do private(gz1,gz2,i,j,lmap1,p1,p2),shared(pslp)
440 loop320:DO km=1,kmm
441  i = imnt(km,lsm)
442  j = jmnt(km,lsm)
443 
444  if(pslp(i,j)<spval) then
445 
446  IF(done(i,j)) cycle
447  lmhij = lmho(i,j)
448  gz1 = fipres(i,j,lmhij)
449  p1(i,j) = spl(lmhij)
450 !
451  lmap1 = lmhij+1
452  DO l=lmap1,lsm
453  p2 = spl(l)
454  tlyr = 0.5*(tpres(i,j,l)+tpres(i,j,l-1))
455  gz2 = gz1 + rd*tlyr*log(p1(i,j)/p2)
456  fipres(i,j,l) = gz2
457 ! if(i==ii.and.j==jj)print*,'Debug:L,FI A S2=',L,GZ2
458  IF(gz2 <= 0.)THEN
459  pslp(i,j) = p1(i,j)/exp(-gz1/(rd*tpres(i,j,l-1)))
460 ! if(i==ii.and.j==jj)print*,'Debug:PSLP A S2=',PSLP(I,J)
461  done(i,j) = .true.
462  kount = kount + 1
463  cycle loop320
464  ENDIF
465  p1(i,j) = p2
466  gz1 = gz2
467  ENDDO
468 !HC EXPERIMENT
469  lp = lsm
470  slope = -6.6e-4
471  tlyr = tpres(i,j,lp)-0.5*fipres(i,j,lp)*slope
472  pslp(i,j) = spl(lp)/exp(-fipres(i,j,lp)/(rd*tlyr))
473  done(i,j) = .true.
474 ! if(i==ii.and.j==jj)print*,'Debug:spl,FI,TLYR,PSLPA3=' &
475 ! ,spl(lp),FIPRES(I,J,LP),TLYR,PSLP(I,J)
476 !HC EXPERIMENT
477  end if ! spval
478 
479 ENDDO loop320
480  320 CONTINUE
481 !
482 !*** WHEN SEA LEVEL IS BELOW THE LOWEST OUTPUT PRESSURE LEVEL,
483 !*** SOLVE THE HYDROSTATIC EQUATION BY CHOOSING A TEMPERATURE
484 !*** AT THE MIDPOINT OF THE LAYER BETWEEN THAT LOWEST PRESSURE
485 !*** LEVEL AND THE GROUND BY EXTRAPOLATING DOWNWARD FROM T ON
486 !*** THE LOWEST PRESSURE LEVEL USING THE DT/DFI BETWEEN THE
487 !*** LOWEST PRESSURE LEVEL AND THE ONE ABOVE IT.
488 !
489 ! TOTAL=(IM-2)*(JM-4)
490 !
491 !HC DO 340 LP=LSM,1,-1
492 ! IF(KOUNT==TOTAL)GO TO 350
493 !HC MODIFICATION FOR SMALL HILL HIGH PRESSURE SITUATION
494 !HC IF SURFACE PRESSURE IS CLOSER TO SEA LEVEL THAN LWOEST
495 !HC OUTPUT PRESSURE LEVEL, USE SURFACE PRESSURE TO DO EXTRAPOLATION
496 
497  325 CONTINUE
498  lp = lsm
499  DO j=jsta,jend
500  DO i=1,im
501 
502  if(pslp(i,j)<spval) then
503 
504 ! if(i==ii.and.j==jj)print*,'Debug: with 330 loop'
505  IF(done(i,j)) cycle
506 
507 ! if(i==ii.and.j==jj)print*,'Debug: still within 330 loop'
508 !HC Comment out the following line for situation with terrain
509 !HC at boundary (ie FIPRES<0)
510 !HC because they were not counted as undergound point for 8 pt
511 !HC relaxation
512 !HC IF(FIPRES(I,J,LP)<0.)GO TO 330
513 ! IF(FIPRES(I,J,LP)<0.)THEN
514 ! DO LP=LSM,1,-1
515 ! IF (FIPRES(I,J) <= 0)
516 
517 ! IF(FIPRES(I,J,LP)<0..OR.DONE(I,J))GO TO 330
518 ! SLOPE=(TPRES(I,J,LP)-TPRES(I,J,LP-1))
519 ! & /(FIPRES(I,J,LP)-FIPRES(I,J,LP-1))
520 
521  slope = -6.6e-4
522  IF(pint(i,j,nint(lmh(i,j))+1) > spl(lp))THEN
523  llmh = nint(lmh(i,j))
524  tvrt = t(i,j,llmh)*(h1+d608*q(i,j,llmh))
525  dis = zint(i,j,llmh+1)-zint(i,j,llmh)+0.5*zint(i,j,llmh+1)
526  tlyr = tvrt-dis*g*slope
527  pslp(i,j) = pint(i,j,llmh+1)*exp(zint(i,j,llmh+1)*g &
528  /(rd*tlyr))
529 ! if(i==ii.and.j==jj)print*,'Debug:PSFC,zsfc,TLYR,PSLPA3='
530 ! 1,PINT(I,J,LLMH+1),ZINT(I,J,LLMH+1),TLYR,PSLP(I,J)
531  ELSE
532  tlyr=tpres(i,j,lp)-0.5*fipres(i,j,lp)*slope
533  pslp(i,j)=spl(lp)/exp(-fipres(i,j,lp)/(rd*tlyr))
534 ! if(i==ii.and.j==jj)print*,'Debug:spl,FI,TLYR,PSLPA3=' &
535 ! ,spl(lp),FIPRES(I,J,LP),TLYR,PSLP(I,J)
536  END IF
537  done(i,j) = .true.
538  kount = kount + 1
539  end if ! spval
540 
541  enddo
542  enddo
543 !HC 340 CONTINUE
544 !
545 ! 350 CONTINUE
546 !--------------------------------------------------------------------
547  RETURN
548  END
masks
Definition: MASKS_mod.f:1
params_mod
Definition: params.F:1