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