CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
module_BL_MYJPBL.F90
1!-----------------------------------------------------------------------
2!
4!
5!-----------------------------------------------------------------------
6!
7!*** THE MYJ PBL SCHEME
8!
9!-----------------------------------------------------------------------
10!
11! USE MODULE_INCLUDE
12!
13! USE MODULE_CONSTANTS,ONLY : A2,A3,A4,CP,ELIV,ELWV,ELIWV &
14! ,EP_1,EPSQ &
15! ,G,P608,PI,PQ0,R_D,R_V,RHOWATER &
16! ,STBOLT,CAPPA
17
18 USE machine, only: kfpt => kind_phys, &
19 kint => kind_integer, &
20 klog => kind_logical
21
22!-----------------------------------------------------------------------
23!
24 IMPLICIT NONE
25!
26!-----------------------------------------------------------------------
27! integer,parameter :: isingle=selected_int_kind(r=9)
28! integer,parameter :: idouble=selected_int_kind(r=18)
29! integer,parameter :: single=selected_real_kind(p=6,r=37)
30! integer,parameter :: double=selected_real_kind(p=13,r=200)
31
32! integer,parameter:: &
33! klog=4 &
34! ,kint=isingle &
35! ,kdin=idouble &
36! ,kfpt=single &
37! ,kdbl=double
38
39! real (kind=kfpt),parameter :: r4_in=x'ffbfffff'
40! real (kind=kdbl),parameter :: r8_in=x'fff7ffffffffffff'
41! integer(kind=kint),parameter :: i4_in=-999 ! -huge(1)
42
43 ! integer,parameter:: &
44 ! klog=4 & ! logical variables
45 ! ,kint=4 & ! integer variables
46 ! !,kfpt=4 & ! floating point variables
47 ! ,kfpt=8 & ! floating point variables
48 ! ,kdbl=8 ! double precision
49
50 REAL(kind=kfpt),PARAMETER :: a2=17.2693882,a3=273.15,a4=35.86,cp=1004.6 &
51 ,eliv=2.850e6,elwv=2.501e6,r_v=461.6 &
52! ,EPSQ=1.e-12,EPSQ2=0.02,G=9.8060226 &
53 ,epsq=1.e-12,g=9.8060226 &
54 ,pq0=379.90516,r_d=287.04,ep_1=r_v/r_d-1. &
55 ,p608=r_v/r_d-1.,pi=3.141592653589793 &
56 ,rhowater=1000.,stbolt=5.67051e-8,cappa=r_d/cp
57 REAL(kind=kfpt),PARAMETER :: eliwv=2.683e6
58!
59 REAL(kind=kfpt),PARAMETER :: conw=1./g,cont=cp/g,conq=elwv/g
60
61!-----------------------------------------------------------------------
62!
63 PRIVATE
64!
65 PUBLIC:: myjpbl_init, myjpbl
66!
67!-----------------------------------------------------------------------
68!-----------------------------------------------------------------------
69!*** FOR MYJ TURBULENCE
70!-----------------------------------------------------------------------
71!-----------------------------------------------------------------------
72!
73 REAL(kind=kfpt),PARAMETER:: &
74 elevfc=0.6
75!
76 REAL(kind=kfpt),PARAMETER:: &
77 vkarman=0.4 &
78!
79 ,xls=eliv,xlv=elwv &
80 ,rlivwv=xls/xlv,elocp=2.72e6/cp &
81!
82 ,eps1=1.e-12,eps2=0. &
83 ,epsru=1.e-7,epsrs=1.e-7 &
84 ,epstrb=1.e-24 &
85 ,fh=1.10 &
86!
87 ,alph=0.30,beta=1./273.,el0max=1000.,el0min=1. &
88! ,ELFC=0.5,GAM1=0.2222222222222222222 &
89! ,ELFC=0.23*0.25,GAM1=0.2222222222222222222 &
90 ,elfc=1.,gam1=0.2222222222222222222 &
91!
92 ,a1=0.659888514560862645 &
93 ,a2x=0.6574209922667784586 &
94 ,b1=11.87799326209552761 &
95 ,b2=7.226971804046074028 &
96 ,c1=0.000830955950095854396 &
97 ,elz0=0.,esq=5.0 &
98!
99 ,seafc=0.98,pq0sea=pq0*seafc &
100!
101 ,btg=beta*g &
102 ,esqhf=0.5*5.0 &
103 ,rb1=1./b1
104!
105 REAL(kind=kfpt),PARAMETER:: &
106 adnh= 9.*a1*a2x*a2x*(12.*a1+3.*b2)*btg*btg &
107 ,adnm=18.*a1*a1*a2x*(b2-3.*a2x)*btg &
108 ,anmh=-9.*a1*a2x*a2x*btg*btg &
109 ,anmm=-3.*a1*a2x*(3.*a2x+3.*b2*c1+18.*a1*c1-b2)*btg &
110 ,bdnh= 3.*a2x*(7.*a1+b2)*btg &
111 ,bdnm= 6.*a1*a1 &
112 ,beqh= a2x*b1*btg+3.*a2x*(7.*a1+b2)*btg &
113 ,beqm=-a1*b1*(1.-3.*c1)+6.*a1*a1 &
114 ,bnmh=-a2x*btg &
115 ,bnmm=a1*(1.-3.*c1) &
116 ,bshh=9.*a1*a2x*a2x*btg &
117 ,bshm=18.*a1*a1*a2x*c1 &
118 ,bsmh=-3.*a1*a2x*(3.*a2x+3.*b2*c1+12.*a1*c1-b2)*btg &
119 ,cesh=a2x &
120 ,cesm=a1*(1.-3.*c1) &
121 ,cnv=ep_1*g/btg
122!
123!-----------------------------------------------------------------------
124!*** FREE TERM IN THE EQUILIBRIUM EQUATION FOR (L/Q)**2
125!-----------------------------------------------------------------------
126!
127 REAL(kind=kfpt),PARAMETER:: &
128 aeqh=9.*a1*a2x*a2x*b1*btg*btg &
129 +9.*a1*a2x*a2x*(12.*a1+3.*b2)*btg*btg &
130 ,aeqm=3.*a1*a2x*b1*(3.*a2x+3.*b2*c1+18.*a1*c1-b2) &
131 *btg+18.*a1*a1*a2x*(b2-3.*a2x)*btg
132!
133!-----------------------------------------------------------------------
134!*** FORBIDDEN TURBULENCE AREA
135!-----------------------------------------------------------------------
136!
137 REAL(kind=kfpt),PARAMETER:: &
138 requ=-aeqh/aeqm &
139 ,epsgh=1.e-9,epsgm=requ*epsgh
140!
141!-----------------------------------------------------------------------
142!*** NEAR ISOTROPY FOR SHEAR TURBULENCE, WW/Q2 LOWER LIMIT
143!-----------------------------------------------------------------------
144!
145 REAL(kind=kfpt),PARAMETER:: &
146 ubryl=(18.*requ*a1*a1*a2x*b2*c1*btg &
147 +9.*a1*a2x*a2x*b2*btg*btg) &
148 /(requ*adnm+adnh) &
149 ,ubry=(1.+epsrs)*ubryl,ubry3=3.*ubry
150!
151 REAL(kind=kfpt),PARAMETER:: &
152 aubh=27.*a1*a2x*a2x*b2*btg*btg-adnh*ubry3 &
153 ,aubm=54.*a1*a1*a2x*b2*c1*btg -adnm*ubry3 &
154 ,bubh=(9.*a1*a2x+3.*a2x*b2)*btg-bdnh*ubry3 &
155 ,bubm=18.*a1*a1*c1 -bdnm*ubry3 &
156 ,cubr=1. - ubry3 &
157 ,rcubr=1./cubr
158!
159!-----------------------------------------------------------------------
160!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
161!---LOOK-UP TABLES------------------------------------------------------
162INTEGER(KIND=KINT),PARAMETER:: &
163 itbl=401 & ! CONVECTION TABLES, DIMENSION 1
164,jtbl=1201 & ! CONVECTION TABLES, DIMENSION 2
165,kerfm=301 & ! SIZE OF ERF HALF TABLE
166,kerfm2=kerfm-2 ! INTERNAL POINTS OF ERF HALF TABLE
167
168REAL(kind=kfpt),PARAMETER:: &
169 pl=2500. & ! LOWER BOUND OF PRESSURE RANGE
170,ph=105000. & ! UPPER BOUND OF PRESSURE RANGE
171,thl=210. & ! LOWER BOUND OF POTENTIAL TEMPERATURE RANGE
172,thh=365. & ! UPPER BOUND OF POTENTIAL TEMPERATURE RANGE
173,xemin=0. & ! LOWER BOUND OF ERF HALF TABLE
174,xemax=3. ! UPPER BOUND OF ERF HALF TABLE
175
176REAL(kind=kfpt),PRIVATE,SAVE:: &
177 rdp & ! SCALING FACTOR FOR PRESSURE
178,rdq & ! SCALING FACTOR FOR HUMIDITY
179,rdth & ! SCALING FACTOR FOR POTENTIAL TEMPERATURE
180,rdthe & ! SCALING FACTOR FOR EQUIVALENT POT. TEMPERATURE
181,rdxe ! ERF HALF TABLE SCALING FACTOR
182
183REAL(kind=kfpt),DIMENSION(1:ITBL),PRIVATE,SAVE:: &
184 sthe & ! RANGE FOR EQUIVALENT POTENTIAL TEMPERATURE
185,the0 ! BASE FOR EQUIVALENT POTENTIAL TEMPERATURE
186
187REAL(kind=kfpt),DIMENSION(1:JTBL),PRIVATE,SAVE:: &
188 qs0 & ! BASE FOR SATURATION SPECIFIC HUMIDITY
189,sqs ! RANGE FOR SATURATION SPECIFIC HUMIDITY
190
191REAL(kind=kfpt),DIMENSION(1:KERFM),PRIVATE,SAVE:: &
192 herff ! HALF ERF TABLE
193
194REAL(kind=kfpt),DIMENSION(1:ITBL,1:JTBL),PRIVATE,SAVE:: &
195 ptbl ! SATURATION PRESSURE TABLE
196
197REAL(kind=kfpt),DIMENSION(1:JTBL,1:ITBL),PRIVATE,SAVE:: &
198 ttbl ! TEMPERATURE TABLE
199!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
200!-----------------------------------------------------------------------
201!
202 CONTAINS
203!
204!-----------------------------------------------------------------------
205!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
206!-----------------------------------------------------------------------
207!
208! REFERENCES: JANJIC (2001), NCEP OFFICE NOTE 437
209!
210! ABSTRACT:
211! MYJ UPDATES THE TURBULENT KINETIC ENERGY WITH THE PRODUCTION/
212! DISSIPATION TERM AND THE VERTICAL DIFFUSION TERM
213! (USING AN IMPLICIT FORMULATION) FROM MELLOR-YAMADA
214! LEVEL 2.5 AS EXTENDED BY JANJIC. EXCHANGE COEFFICIENTS FOR
215! THE SURFACE LAYER ARE COMPUTED FROM THE MONIN-OBUKHOV THEORY.
216! THE TURBULENT VERTICAL EXCHANGE IS THEN EXECUTED.
217!
218!-----------------------------------------------------------------------
219 SUBROUTINE myjpbl(NTSD,ME,DT_PHS,EPSL,EPSQ2,HT,STDH,DZ,DEL &
220 ,PMID,PINH,TH,T,EXNER,Q,CWM,U,V &
221 ,TSK,QSFC,CHKLOWQ,THZ0,QZ0,UZ0,VZ0 &
222 ,XLAND,SICE,SNOW &
223 ,Q2,EXCH_H,USTAR,Z0,EL_MYJ,PBLH,KPBL,CT &
224 ,AKHS,AKMS,ELFLX,MIXHT,THLM,QLM &
225 ,RUBLTEN,RVBLTEN,RTHBLTEN,RQBLTEN,RQCBLTEN &
226 ,DUSFC,DVSFC,DTSFC,DQSFC,xkzo,xkzmo,ICT &
227 ,IDS,IDE,JDS,JDE &
228 ,IMS,IME,JMS,JME &
229 ,ITS,ITE,JTS,JTE,LM)
230
231! SUBROUTINE MYJPBL(DT,NPHS,EPSL,EPSQ2,HT,STDH,DZ &
232! ,PMID,PINH,TH,T,EXNER,Q,CWM,U,V &
233! ,TSK,QSFC,CHKLOWQ,THZ0,QZ0,UZ0,VZ0 &
234! ,XLAND,SICE,SNOW &
235! ,Q2,EXCH_H,USTAR,Z0,EL_MYJ,PBLH,KPBL,CT &
236! ,AKHS,AKMS,ELFLX,MIXHT &
237! ,RUBLTEN,RVBLTEN,RTHBLTEN,RQBLTEN,RQCBLTEN &
238! ,IDS,IDE,JDS,JDE &
239! ,IMS,IME,JMS,JME &
240! ,ITS,ITE,JTS,JTE,LM)
241
242!----------------------------------------------------------------------
243!
244 IMPLICIT NONE
245!
246 logical(kind=klog),save:: &
247 reinit
248!----------------------------------------------------------------------
249 INTEGER(KIND=KINT),INTENT(IN):: &
250 ids,ide,jds,jde &
251 ,ims,ime,jms,jme &
252 ,its,ite,jts,jte,lm
253!
254 INTEGER,INTENT(IN) :: ict,me,ntsd
255
256! INTEGER(KIND=KINT),INTENT(IN):: &
257! NPHS
258!
259 INTEGER(KIND=KINT),DIMENSION(IMS:IME,JMS:JME),INTENT(OUT):: &
260 kpbl
261!
262 REAL(kind=kfpt),INTENT(IN):: &
263 dt_phs
264! DT
265!
266 real(kind=kfpt),dimension(1:lm-1),intent(inout):: epsl
267 real(kind=kfpt),dimension(1:lm),intent(in):: epsq2
268!
269 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME),INTENT(IN):: &
270 ht,sice,snow,stdh &
271 ,tsk,xland &
272 ,chklowq,elflx,thlm,qlm
273!
274 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME,1:LM),INTENT(IN):: &
275 dz,exner,pmid,q,cwm,u,v,t,th,del,xkzo,xkzmo
276!
277 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME,1:LM+1),INTENT(IN):: &
278 pinh
279!
280 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME),INTENT(OUT):: &
281 mixht &
282 ,pblh
283!
284 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME,1:LM),INTENT(OUT):: &
285 el_myj
286!
287 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME,1:LM),INTENT(OUT):: &
288 rqcblten &
289 ,rublten,rvblten &
290 ,rthblten,rqblten
291!
292 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: &
293 dusfc,dvsfc &
294 ,dtsfc,dqsfc
295!
296 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT):: &
297 akhs,akms
298!
299 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT):: &
300 ct,qsfc,qz0 &
301 ,thz0,ustar &
302 ,uz0,vz0,z0
303!
304 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME,1:LM),INTENT(INOUT):: &
305 exch_h &
306 ,q2
307!
308!----------------------------------------------------------------------
309!***
310!*** LOCAL VARIABLES
311!***
312 INTEGER(KIND=KINT):: &
313 i,iqtb,ittb,j,k,llow,lmh,lmxl
314!
315 INTEGER(KIND=KINT),DIMENSION(IMS:IME,JMS:JME):: &
316 lpbl
317!
318 REAL(kind=kfpt):: &
319 akhs_dens,akms_dens,bq,bqs00k,bqs10k &
320 ,dcdt,deltaz,dqdt,dtdif,dtdt,dtturbl &
321 ,p00k,p01k,p10k,p11k,pelevfc,pp1,psfc,psp,ptop &
322 ,qbt,qfc1,qlow,qq1,qx &
323 ,rdtturbl,rg,rsqdt,rxners,rxnsfc &
324 ,seamask,sq,sqs00k,sqs10k &
325 ,thbt,thnew,thold,tq,tth &
326 ,ulow,vlow,rstdh,stdfac,zsf,zsx,zsy,zuv
327!
328 REAL(kind=kfpt),DIMENSION(1:LM):: &
329 cwmk,pk,psk,q2k,qk,rhok,rxnerk,thek,thk,thvk,tk,uk,vk
330!
331 REAL(kind=kfpt),DIMENSION(1:LM-1):: &
332 akhk,akmk,dcol,el,gh,gm
333!
334 REAL(kind=kfpt),DIMENSION(1:LM+1):: &
335 zhk
336!
337 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME):: &
338 thsk
339!
340 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME,1:LM):: &
341 rxner,thv
342!
343 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME,1:LM-1):: &
344 akh,akm
345!
346 REAL(kind=kfpt),DIMENSION(IMS:IME,JMS:JME,1:LM+1):: &
347 zint
348!
349!*** Begin debugging
350 REAL(kind=kfpt):: zsl_diag
351 INTEGER(KIND=KINT):: imd,jmd,print_diag
352!*** End debugging
353!-----------------------------------------------------------------------
354!***********************************************************************
355 data reinit/.false./
356!-----------------------------------------------------------------------
357! if(reinit) then
358! call MYJPBL_INIT( &
359! 1,IDE,1,1,LM, &
360! 1,IDE,1,1, &
361! 1,IDE,1,1)
362! reinit=.false.
363! endif
364!
365!----------------------------------------------------------------------
366!**********************************************************************
367!----------------------------------------------------------------------
368!
369!*** Begin debugging
370 imd=(ims+ime)/2
371 jmd=(jms+jme)/2
372!*** End debugging
373!
374!*** MAKE PREPARATIONS
375!
376!----------------------------------------------------------------------
377 stdfac=1.
378!----------------------------------------------------------------------
379! DTTURBL=DT*NPHS
380 dtturbl=dt_phs
381 rdtturbl=1./dtturbl
382 rsqdt=sqrt(rdtturbl)
383 dtdif=dtturbl
384 rg=1./g
385!
386 DO k=1,lm-1
387 DO j=jts,jte
388 DO i=its,ite
389 akm(i,j,k)=0.
390 ENDDO
391 ENDDO
392 ENDDO
393!
394 DO k=1,lm+1
395 DO j=jts,jte
396 DO i=its,ite
397 zint(i,j,k)=0.
398 ENDDO
399 ENDDO
400 ENDDO
401!
402 DO j=jts,jte
403 DO i=its,ite
404 zint(i,j,lm+1)=ht(i,j) ! Z AT BOTTOM OF LOWEST SIGMA LAYER
405 ENDDO
406 ENDDO
407!
408 DO k=lm,1,-1
409 DO j=jts,jte
410 DO i=its,ite
411 zint(i,j,k)=zint(i,j,k+1)+dz(i,j,k)
412 rxner(i,j,k)=1./exner(i,j,k)
413 thv(i,j,k)=(q(i,j,k)*0.608+(1.-cwm(i,j,k)))*th(i,j,k)
414 ENDDO
415 ENDDO
416 ENDDO
417!
418 DO j=jts,jte
419 DO i=its,ite
420 el_myj(i,j,lm)=0.
421 ENDDO
422 ENDDO
423 DO j=jts,jte
424 DO i=its,ite
425 dusfc(i,j)=0.
426 dvsfc(i,j)=0.
427 dtsfc(i,j)=0.
428 dqsfc(i,j)=0.
429 ENDDO
430 ENDDO
431
432!
433!----------------------------------------------------------------------
434!.......................................................................
435!ZJ$OMP PARALLEL DO &
436!ZJ$OMP PRIVATE(J,I,LMH,PTOP,PSFC,SEAMASK,K,TK,THVK,QK,Q2K,RXNERK, &
437!ZJ$OMP PK,UK,VK,Q2K,ZHK,LMXL,GM,GH,EL,AKMK,AKHK,DELTAZ), &
438!ZJ$OMP SCHEDULE(DYNAMIC)
439!.......................................................................
440!----------------------------------------------------------------------
441 setup_integration: DO j=jts,jte
442!----------------------------------------------------------------------
443!
444 DO i=its,ite
445!
446 lmh=lm
447!
448 ptop=pinh(i,j,1)
449 psfc=pinh(i,j,lmh+1)
450!
451!*** CONVERT LAND MASK (1 FOR SEA; 0 FOR LAND)
452!
453 seamask=xland(i,j)-1.
454!
455!*** FILL 1-D VERTICAL ARRAYS
456!
457 DO k=lm,1,-1
458 pk(k)=pmid(i,j,k)
459 tk(k)=t(i,j,k)
460 qk(k)=q(i,j,k)
461 thvk(k)=thv(i,j,k)
462 rxnerk(k)=rxner(i,j,k)
463 uk(k)=u(i,j,k)
464 vk(k)=v(i,j,k)
465 q2k(k)=q2(i,j,k)
466!
467!*** COMPUTE THE HEIGHTS OF THE LAYER INTERFACES
468!
469 zhk(k)=zint(i,j,k)
470!
471 ENDDO
472 zhk(lm+1)=ht(i,j) ! Z AT BOTTOM OF LOWEST SIGMA LAYER
473!
474!*** POTENTIAL INSTABILITY
475!
476 pelevfc=pmid(i,j,lmh)*elevfc
477!
478 DO k=lmh,1,-1
479!-----------------------------------------------------------------------
480 IF(k==lmh .OR. pmid(i,j,k)>pelevfc) THEN
481!---PREPARATION FOR SEARCH FOR MAX CAPE---------------------------------
482 qbt=qk(k)
483 thbt=th(i,j,k)
484 tth=(thbt-thl)*rdth
485 qq1=tth-aint(tth)
486 ittb=int(tth)+1
487!---KEEPING INDICES WITHIN THE TABLE------------------------------------
488 IF(ittb.LT.1)THEN
489 ittb=1
490 qq1=0.
491 ELSE IF(ittb.GE.jtbl)THEN
492 ittb=jtbl-1
493 qq1=0.
494 ENDIF
495!---BASE AND SCALING FACTOR FOR SPEC. HUMIDITY--------------------------
496 bqs00k=qs0(ittb)
497 sqs00k=sqs(ittb)
498 bqs10k=qs0(ittb+1)
499 sqs10k=sqs(ittb+1)
500!--------------SCALING SPEC. HUMIDITY & TABLE INDEX---------------------
501 bq=(bqs10k-bqs00k)*qq1+bqs00k
502 sq=(sqs10k-sqs00k)*qq1+sqs00k
503 tq=(qbt-bq)/sq*rdq
504 pp1=tq-aint(tq)
505 iqtb=int(tq)+1
506!----------------KEEPING INDICES WITHIN THE TABLE-----------------------
507 IF(iqtb.LT.1)THEN
508 iqtb=1
509 pp1=0.
510 ELSEIF(iqtb.GE.itbl)THEN
511 iqtb=itbl-1
512 pp1=0.
513 ENDIF
514!--------------SATURATION PRESSURE AT FOUR SURROUNDING TABLE PTS.-------
515 p00k=ptbl(iqtb ,ittb )
516 p10k=ptbl(iqtb+1,ittb )
517 p01k=ptbl(iqtb ,ittb+1)
518 p11k=ptbl(iqtb+1,ittb+1)
519!--------------SATURATION POINT VARIABLES AT THE BOTTOM-----------------
520 psp=p00k+(p10k-p00k)*pp1+(p01k-p00k)*qq1 &
521 +(p00k-p10k-p01k+p11k)*pp1*qq1
522 rxners=(1.e5/psp)**cappa
523 thek(k)=thbt*exp(elocp*qbt*rxners/thbt)
524 psk(k)=psp
525!-----------------------------------------------------------------------
526 ELSE
527!-----------------------------------------------------------------------
528 thek(k)=thek(k+1)
529 psk(k)=pinh(i,j,1)
530!-----------------------------------------------------------------------
531 ENDIF
532!-----------------------------------------------------------------------
533 ENDDO
534!
535!*** Begin debugging
536! IF(I==IMD.AND.J==JMD)THEN
537! PRINT_DIAG=1
538! ELSE
539! PRINT_DIAG=0
540! ENDIF
541! IF(I==227.AND.J==363)PRINT_DIAG=2
542!*** End debugging
543!
544!----------------------------------------------------------------------
545!***
546!*** FIND THE MIXING LENGTH
547!***
548 CALL mixlen(lmh,rsqdt,uk,vk,thvk,thek &
549 ,q2k,epsl,epsq2,zhk,pk,psk,rxnerk,gm,gh,el &
550 ,pblh(i,j),lpbl(i,j),lmxl,ct(i,j),mixht(i,j) &
551 ,i,j,lm)
552!
553!----------------------------------------------------------------------
554!***
555!*** SOLVE FOR THE PRODUCTION/DISSIPATION OF
556!*** THE TURBULENT KINETIC ENERGY
557!***
558!
559 CALL prodq2(ntsd,me,lmh,dtturbl,ustar(i,j),gm,gh,el,q2k &
560 ,epsl,epsq2,i,j,lm)
561
562! if(i.eq.4)print*,'11ql test Q2(LMH)=',Q2K(LMH),B1,USTAR(I,J)
563!
564!----------------------------------------------------------------------
565!*** THE MODEL LAYER (COUNTING UPWARD) CONTAINING THE TOP OF THE PBL
566!----------------------------------------------------------------------
567!
568 kpbl(i,j)=lpbl(i,j)
569!
570!----------------------------------------------------------------------
571!***
572!*** FIND THE EXCHANGE COEFFICIENTS IN THE FREE ATMOSPHERE
573!***
574 CALL difcof(ntsd,me,lmh,lmxl,gm,gh,el,tk,q2k,zhk,akmk,akhk,i,j,lm &
575 ,print_diag,kpbl(i,j))
576!
577!*** COUNTING DOWNWARD FROM THE TOP, THE EXCHANGE COEFFICIENTS AKH
578!*** ARE DEFINED ON THE BOTTOMS OF THE LAYERS 1 TO LM-1. COUNTING
579!*** COUNTING UPWARD FROM THE BOTTOM, THOSE SAME COEFFICIENTS EXCH_H
580!*** ARE DEFINED ON THE TOPS OF THE LAYERS 1 TO LM-1.
581!
582 DO k=1,lm-1
583
584 deltaz=0.5*(zhk(k)-zhk(k+2))
585 akhk(k)=max(akhk(k),xkzo(i,j,k)/deltaz) ! add minimum background diffusion
586 akmk(k)=max(akmk(k),xkzmo(i,j,k)/deltaz)
587 if((thvk(lm)-thvk(k)).GT.0.) then
588 akhk(k)=max(akhk(k),3./deltaz) ! add minimum background diffusion
589 akmk(k)=max(akmk(k),3./deltaz)
590 end if
591 akh(i,j,k)=akhk(k)
592 akm(i,j,k)=akmk(k)
593 exch_h(i,j,k)=akhk(k)*deltaz
594 ENDDO
595!
596!----------------------------------------------------------------------
597!***
598!*** CARRY OUT THE VERTICAL DIFFUSION OF
599!*** TURBULENT KINETIC ENERGY
600!***
601!
602 CALL vdifq(lmh,dtdif,q2k,el,zhk,i,j,lm)
603!
604!*** SAVE THE NEW Q2 AND MIXING LENGTH.
605!
606 DO k=1,lm
607 q2(i,j,k)=max(q2k(k),epsq2(k))
608 IF(k<lm)el_myj(i,j,k)=el(k) ! EL IS NOT DEFINED AT LM
609 ENDDO
610!
611 ENDDO
612!
613!----------------------------------------------------------------------
614!
615 ENDDO setup_integration
616!
617!.......................................................................
618!ZJ$OMP END PARALLEL DO
619!.......................................................................
620!----------------------------------------------------------------------
621!
622!*** CONVERT SURFACE SENSIBLE TEMPERATURE TO POTENTIAL TEMPERATURE.
623!
624 DO j=jts,jte
625 DO i=its,ite
626 psfc=pinh(i,j,lm+1)
627 thsk(i,j)=tsk(i,j)*(1.e5/psfc)**cappa
628 ENDDO
629 ENDDO
630!
631!----------------------------------------------------------------------
632!
633!----------------------------------------------------------------------
634!.......................................................................
635!ZJ$OMP PARALLEL DO PRIVATE(I,J,K &
636!ZJ$OMP & ,THK,QK,CWMK,ZHK,RHOK &
637!ZJ$OMP & ,AKHK,SEAMASK,LLOW,AKHS_DENS,QFC1,QLOW,PSFC,RXNSFC,LMH &
638!ZJ$OMP & ,THOLD,THNEW,DTDT,DQDT,DCDT,ZSL_DIAG,AKMK,AKMS_DENS &
639!ZJ$OMP & ,UK,VK &
640!ZJ$OMP & ),SCHEDULE(DYNAMIC)
641!.......................................................................
642!----------------------------------------------------------------------
643 main_integration: DO j=jts,jte
644!----------------------------------------------------------------------
645!
646 DO i=its,ite
647!
648!*** FILL 1-D VERTICAL ARRAYS
649!
650 DO k=lm,1,-1
651 thk(k)=th(i,j,k)
652 qk(k)=q(i,j,k)
653 cwmk(k)=cwm(i,j,k)
654 zhk(k)=zint(i,j,k)
655 rhok(k)=pmid(i,j,k)/(r_d*t(i,j,k)*(1.+p608*qk(k)-cwmk(k)))
656 ENDDO
657!
658!*** COUNTING DOWNWARD FROM THE TOP, THE EXCHANGE COEFFICIENTS AKH
659!*** ARE DEFINED ON THE BOTTOMS OF THE LAYERS 1 TO LM-1. THESE COEFFICIENTS
660!*** ARE ALSO MULTIPLIED BY THE DENSITY AT THE BOTTOM INTERFACE LEVEL.
661!
662 DO k=1,lm-1
663 akhk(k)=akh(i,j,k)*0.5*(rhok(k)+rhok(k+1))
664 ENDDO
665!
666 zhk(lm+1)=zint(i,j,lm+1)
667!
668 seamask=xland(i,j)-1.
669 thz0(i,j)=(1.-seamask)*thsk(i,j)+seamask*thz0(i,j)
670!
671 llow=lm
672 akhs_dens=akhs(i,j)*rhok(lm)
673!
674 IF(seamask<0.5)THEN
675 qfc1=xlv*chklowq(i,j)*akhs_dens
676!
677 IF(snow(i,j)>0..OR.sice(i,j)>0.5)THEN
678 qfc1=qfc1*rlivwv
679 ENDIF
680!
681 IF(qfc1>0.)THEN
682 qlow=qk(lm)
683!ql QSFC(I,J)=QLOW+ELFLX(I,J)/QFC1
684 ENDIF
685!
686 ELSE
687 psfc=pinh(i,j,lm+1)
688 rxnsfc=(1.e5/psfc)**cappa
689
690!ql QSFC(I,J)=PQ0SEA/PSFC &
691!ql & *EXP(A2*(THSK(I,J)-A3*RXNSFC)/(THSK(I,J)-A4*RXNSFC))
692 ENDIF
693!
694 qz0(i,j)=(1.-seamask)*qsfc(i,j)+seamask*qz0(i,j)
695!
696 lmh=lm
697!
698!----------------------------------------------------------------------
699!*** CARRY OUT THE VERTICAL DIFFUSION OF
700!*** TEMPERATURE AND WATER VAPOR
701!----------------------------------------------------------------------
702!
703 CALL vdifh(dtdif,lmh,thz0(i,j),qz0(i,j) &
704 ,akhs_dens,chklowq(i,j),ct(i,j) &
705 ,thk,qk,cwmk,akhk,zhk,rhok,i,j,lm)
706!----------------------------------------------------------------------
707!***
708! QL set lower bondary
709! THK(LM)=THLM(I,J)
710! QK(LM)=QLM(I,J)
711!*** COMPUTE PRIMARY VARIABLE TENDENCIES
712!***
713 DO k=1,lm
714 rthblten(i,j,k)=(thk(k)-th(i,j,k))*rdtturbl
715 rqblten(i,j,k)=(qk(k)-q(i,j,k))*rdtturbl
716 rqcblten(i,j,k)=(cwmk(k)-cwm(i,j,k))*rdtturbl
717 dtsfc(i,j)=dtsfc(i,j)+cont*del(i,j,k)*rthblten(i,j,k)*exner(i,j,k)
718 dqsfc(i,j)=dqsfc(i,j)+conq*del(i,j,k)*rqblten(i,j,k)
719 ENDDO
720!
721!*** Begin debugging
722! IF(I==IMD.AND.J==JMD)THEN
723! PRINT_DIAG=0
724! ELSE
725! PRINT_DIAG=0
726! ENDIF
727! IF(I==227.AND.J==363)PRINT_DIAG=0
728!*** End debugging
729!
730 psfc=.01*pinh(i,j,lm+1)
731 zsl_diag=0.5*dz(i,j,lm)
732!
733!*** Begin debugging
734! IF(PRINT_DIAG==1)THEN
735!
736! WRITE(6,"(A, 2I5, 2I3, 2F8.2, F6.2, 2F8.2)") &
737! '{TURB4 I,J, KPBL, KMXL, PSFC, ZSFC, ZSL, ZPBL, ZMXL = ' &
738! , I, J, KPBL(I,J), LM-LMXL+1, PSFC, ZHK(LMH+1), ZSL_DIAG &
739! , PBLH(I,J), ZHK(LMXL)-ZHK(LMH+1)
740! WRITE(6,"(A, 2F7.2, F7.3, 3E11.4)") &
741! '{TURB4 TSK, THSK, QZ0, Q**2_0, AKHS, EXCH_0 = ' &
742! , TSK(I,J)-273.15, THSK(I,J), 1000.*QZ0(I,J) &
743! , Q2(I,1,J), AKHS(I,J), AKHS(I,J)*ZSL_DIAG
744! WRITE(6,"(A)") &
745! '{TURB5 K, PMID, PINH_1, TC, TH, DTH, GH, GM, EL, Q**2, AKH, EXCH_H, DZ, DP'
746! DO K=1,LM/2
747! WRITE(6,"(A,I3, 2F8.2, 2F8.3, 3E12.4, 4E11.4, F7.2, F6.2)") &
748! '{TURB5 ', K, .01*PMID(I,K,J),.01*PINH(I,K,J), T(I,K,J)-273.15 &
749! , TH(I,K,J), DTTURBL*RTHBLTEN(I,K,J), GH(K), GM(K) &
750! , EL_MYJ(I,K,J), Q2(I,K+1,J), AKH(I,K,J) &
751! , EXCH_H(I,K,J), DZ(I,K,J), .01*(PINH(I,K,J)-PINH(I,K+1,J))
752! ENDDO
753!
754! ELSEIF(PRINT_DIAG==2)THEN
755!
756! WRITE(6,"(A, 2I5, 2I3, 2F8.2, F6.2, 2F8.2)") &
757! '}TURB4 I,J, KPBL, KMXL, PSFC, ZSFC, ZSL, ZPBL, ZMXL = ' &
758! , I, J, KPBL(I,J), LM-LMXL+1, PSFC, ZHK(LMH+1), ZSL_DIAG &
759! , PBLH(I,J), ZHK(LMXL)-ZHK(LMH+1)
760! WRITE(6,"(A, 2F7.2, F7.3, 3E11.4)") &
761! '}TURB4 TSK, THSK, QZ0, Q**2_0, AKHS, EXCH_0 = ' &
762! , TSK(I,J)-273.15, THSK(I,J), 1000.*QZ0(I,J) &
763! , Q2(I,1,J), AKHS(I,J), AKHS(I,J)*ZSL_DIAG
764! WRITE(6,"(A)") &
765! '}TURB5 K, PMID, PINH_1, TC, TH, DTH, GH, GM, EL, Q**2, AKH, EXCH_H, DZ, DP'
766! DO K=1,LM/2
767! WRITE(6,"(A,I3, 2F8.2, 2F8.3, 3E12.4, 4E11.4, F7.2, F6.2)") &
768! '}TURB5 ', K, .01*PMID(I,K,J),.01*PINH(I,K,J), T(I,K,J)-273.15 &
769! , TH(I,K,J), DTTURBL*RTHBLTEN(I,K,J), GH(K), GM(K) &
770! , EL_MYJ(I,K,J), Q2(I,K+1,J), AKH(I,K,J) &
771! , EXCH_H(I,K,J), DZ(I,K,J), .01*(PINH(I,K,J)-PINH(I,K+1,J))
772! ENDDO
773! ENDIF
774!*** End debugging
775!
776!----------------------------------------------------------------------
777!
778 seamask=xland(i,j)-1.
779!
780 IF(seamask.LT.0.5.AND.stdh(i,j).GT.1.) THEN
781 rstdh=1./stdh(i,j)
782 ELSE
783 rstdh=0.
784 ENDIF
785 zhk(lm+1)=zint(i,j,lm+1)
786 zsf=stdh(i,j)*stdfac+zhk(lm+1)
787!
788!----------------------------------------------------------------------
789!
790!*** FILL 1-D VERTICAL ARRAYS
791!
792 DO k=1,lm-1
793 akmk(k)=akm(i,j,k)
794 akmk(k)=akmk(k)*(rhok(k)+rhok(k+1))*0.5
795 ENDDO
796!
797 akms_dens=akms(i,j)*rhok(lm)
798!
799 DO k=lm,1,-1
800 uk(k)=u(i,j,k)
801 vk(k)=v(i,j,k)
802 zhk(k)=zint(i,j,k)
803 ENDDO
804 zhk(lm+1)=zint(i,j,lm+1)
805!
806!----------------------------------------------------------------------
807!
808 DO k=1,lm-1
809!jun23 IF(SEAMASK.GT.0.5) THEN
810!jun23 DCOL(K)=0.
811!jun23 ELSE
812!jun23 ZUV=(ZHK(K)+ZHK(K+1))*0.5
813!jun23 IF(ZUV.GT.ZSF) THEN
814!jun23 DCOL(K)=0.
815!jun23 ELSE
816!jun23 DCOL(K)=HERF((((ZUV-ZHK(LM+1))*RSTDH)**2)*0.5)
817!jun23 ENDIF
818!jun23 ENDIF
819!WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
820 dcol(k)=0. !ZJ
821!MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
822 ENDDO
823!
824!----------------------------------------------------------------------
825!*** CARRY OUT THE VERTICAL DIFFUSION OF
826!*** VELOCITY COMPONENTS
827!----------------------------------------------------------------------
828!
829 CALL vdifv(lmh,dtdif,uz0(i,j),vz0(i,j) &
830 & ,akms_dens,dcol,uk,vk,akmk,zhk,rhok,i,j,lm)
831!
832!----------------------------------------------------------------------
833!***
834!*** COMPUTE PRIMARY VARIABLE TENDENCIES
835!***
836 DO k=1,lm
837 rublten(i,j,k)=(uk(k)-u(i,j,k))*rdtturbl
838 rvblten(i,j,k)=(vk(k)-v(i,j,k))*rdtturbl
839 dusfc(i,j)=dusfc(i,j)+conw*del(i,j,k)*rublten(i,j,k)
840 dvsfc(i,j)=dvsfc(i,j)+conw*del(i,j,k)*rvblten(i,j,k)
841 ENDDO
842!
843 ENDDO
844!----------------------------------------------------------------------
845!
846 ENDDO main_integration
847!JAA!ZJ$OMP END PARALLEL DO
848!
849!----------------------------------------------------------------------
850!
851 END SUBROUTINE myjpbl
852!
853!----------------------------------------------------------------------
854!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
855!----------------------------------------------------------------------
856 SUBROUTINE mixlen &
857!----------------------------------------------------------------------
858! ******************************************************************
859! * *
860! * LEVEL 2.5 MIXING LENGTH *
861! * *
862! ******************************************************************
863!
864 (lmh,rsqdt,u,v,thv,the,q2,epsl,epsq2,z,p,ps,rxner &
865 ,gm,gh,el,pblh,lpbl,lmxl,ct,mixht,i,j,lm)
866!
867!----------------------------------------------------------------------
868!
869 IMPLICIT NONE
870!
871!----------------------------------------------------------------------
872 INTEGER(KIND=KINT),INTENT(IN):: &
873 LMH,I,J,LM
874!
875 REAL(KIND=kfpt),INTENT(IN):: &
876 rsqdt
877!
878 INTEGER(KIND=KINT),INTENT(OUT):: &
879 LMXL,LPBL
880!
881 real(kind=kfpt),dimension(1:lm-1),intent(inout):: epsl
882 REAL(KIND=kfpt),DIMENSION(1:LM),INTENT(IN):: &
883 p,ps,epsq2,rxner,the,thv,u,v
884! P,PS,Q2,EPSQ2,RXNER,THE,THV,U,V
885!
886 REAL(KIND=kfpt),DIMENSION(1:LM),INTENT(INOUT):: q2
887!
888 REAL(KIND=kfpt),DIMENSION(1:LM+1),INTENT(IN):: &
889 z
890!
891 REAL(KIND=kfpt),INTENT(OUT):: &
892 mixht &
893 ,pblh
894!
895 REAL(KIND=kfpt),DIMENSION(1:LM-1),INTENT(OUT):: &
896 el,gh,gm
897!
898 REAL(KIND=kfpt),INTENT(INOUT):: ct
899!----------------------------------------------------------------------
900!***
901!*** LOCAL VARIABLES
902!***
903 INTEGER(KIND=KINT):: &
904 K,LPBLM
905!
906 REAL(KIND=kfpt):: &
907 aden,bden,aubr,bubr,blmx,cubry,dthv,dz &
908 ,el0,eloq2x,ghl,gml &
909 ,qol2st,qol2un,qdzl &
910 ,rdz,sq,srel,szq,vkrmz,wcon
911!
912 REAL(KIND=kfpt),DIMENSION(1:LM):: &
913 q1
914!
915 REAL(KIND=kfpt),DIMENSION(1:LM-1):: &
916 elm,rel
917!
918!----------------------------------------------------------------------
919!***********************************************************************
920!--------1---------2---------3---------4---------5---------6---------7--
921 cubry=ubry*1.5 !*2.
922!--------------FIND THE HEIGHT OF THE PBL-------------------------------
923 lpbl=lmh
924! LPBL=LMH-1
925 DO k=lmh-1,1,-1
926! EPSL(K)=1.
927 if((thv(lmh)-thv(k)).GT.0.) then
928 q2(k)=max(q2(k),1.0)
929 epsl(k)=10.
930 ENDIF
931 ENDDO
932!
933 DO k=lmh-1,1,-1
934 if(q2(k)-epsq2(k)+epsq2(lm).le.epsq2(lm)*fh) then
935 lpbl=k
936 GO TO 110
937 ENDIF
938 ENDDO
939!
940 lpbl=1
941!
942!--------------THE HEIGHT OF THE PBL------------------------------------
943!
944 110 pblh=z(lpbl+1)-z(lmh+1)
945!
946!-----------------------------------------------------------------------
947 DO k=1,lmh
948 q1(k)=0.
949 ENDDO
950!-----------------------------------------------------------------------
951 DO k=1,lmh-1
952 dz=(z(k)-z(k+2))*0.5
953 rdz=1./dz
954 gml=((u(k)-u(k+1))**2+(v(k)-v(k+1))**2)*rdz*rdz
955 gm(k)=max(gml,epsgm)
956!
957 dthv=thv(k)-thv(k+1)
958!----------------------------------------------------------------------
959 IF(dthv.GT.0.) THEN
960 IF(the(k+1).GT.the(k)) THEN
961 IF(ps(k+1).GT.p(k)) THEN
962!
963 wcon=(p(k+1)-ps(k+1))/(p(k+1)-p(k))
964!
965 if( &
966 (q2(k).gt.epsq2(k)) .and. &
967 (q2(k)*cubry.gt.(dz*wcon*rsqdt)**2) &
968 ) then
969!
970 dthv=(the(k)-the(k+1))+dthv
971!
972 ENDIF
973 ENDIF
974 ENDIF
975 ENDIF
976!--------------------------------------------------------------------------
977!
978 ghl=dthv*rdz
979 IF(abs(ghl)<=epsgh)ghl=epsgh
980 gh(k)=ghl
981 ENDDO
982!
983 ct=0.
984!
985!----------------------------------------------------------------------
986!*** FIND MAXIMUM MIXING LENGTHS AND THE LEVEL OF THE PBL TOP
987!----------------------------------------------------------------------
988!
989 lmxl=lmh
990!
991 DO k=1,lmh-1
992 gml=gm(k)
993 ghl=gh(k)
994!
995 IF(ghl>=epsgh)THEN
996 IF(gml/ghl<=requ)THEN
997 elm(k)=epsl(k)
998 lmxl=k+1
999 ELSE
1000 aubr=(aubm*gml+aubh*ghl)*ghl
1001 bubr= bubm*gml+bubh*ghl
1002 qol2st=(-0.5*bubr+sqrt(bubr*bubr*0.25-aubr*cubr))*rcubr
1003 eloq2x=1./max(epsgh, qol2st)
1004 elm(k)=max(sqrt(eloq2x*q2(k)),epsl(k))
1005 ENDIF
1006 ELSE
1007 aden=(adnm*gml+adnh*ghl)*ghl
1008 bden= bdnm*gml+bdnh*ghl
1009 qol2un=-0.5*bden+sqrt(bden*bden*0.25-aden)
1010 eloq2x=1./(qol2un+epsru) ! REPSR1/QOL2UN
1011 elm(k)=max(sqrt(eloq2x*q2(k)),epsl(k))
1012 ENDIF
1013 ENDDO
1014!
1015 IF(elm(lmh-1)==epsl(lmh-1))lmxl=lmh
1016!
1017!----------------------------------------------------------------------
1018!*** THE HEIGHT OF THE MIXED LAYER
1019!----------------------------------------------------------------------
1020!
1021 blmx=z(lmxl)-z(lmh+1)
1022 mixht=blmx
1023!
1024!----------------------------------------------------------------------
1025 DO k=lpbl,lmh
1026 q1(k)=sqrt(q2(k))
1027 ENDDO
1028!----------------------------------------------------------------------
1029 szq=0.
1030 sq =0.
1031!
1032 DO k=1,lmh-1
1033 qdzl=(q1(k)+q1(k+1))*(z(k+1)-z(k+2))
1034 szq=(z(k+1)+z(k+2)-z(lmh+1)-z(lmh+1))*qdzl+szq
1035 sq=qdzl+sq
1036 ENDDO
1037!
1038!----------------------------------------------------------------------
1039!*** COMPUTATION OF ASYMPTOTIC L IN BLACKADAR FORMULA
1040!----------------------------------------------------------------------
1041!
1042 el0=min(alph*szq*0.5/sq,el0max)
1043 el0=max(el0 ,el0min)
1044!
1045!----------------------------------------------------------------------
1046!*** ABOVE THE PBL TOP
1047!----------------------------------------------------------------------
1048!
1049 lpblm=max(lpbl-1,1)
1050!
1051 DO k=1,lpblm
1052 el(k)=min((z(k)-z(k+2))*elfc,elm(k))
1053 rel(k)=el(k)/elm(k)
1054 ENDDO
1055!
1056!----------------------------------------------------------------------
1057!*** INSIDE THE PBL
1058!----------------------------------------------------------------------
1059!
1060 IF(lpbl<lmh)THEN
1061 DO k=lpbl,lmh-1
1062 vkrmz=(z(k+1)-z(lmh+1))*vkarman
1063 el(k)=min(vkrmz/(vkrmz/el0+1.),elm(k))
1064 rel(k)=el(k)/elm(k)
1065 ENDDO
1066 ENDIF
1067!
1068 DO k=lpbl+1,lmh-2
1069 srel=min(((rel(k-1)+rel(k+1))*0.5+rel(k))*0.5,rel(k))
1070 el(k)=max(srel*elm(k),epsl(k))
1071 ENDDO
1072!
1073!----------------------------------------------------------------------
1074 END SUBROUTINE mixlen
1075!----------------------------------------------------------------------
1076!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1077!----------------------------------------------------------------------
1078 SUBROUTINE prodq2 &
1079!----------------------------------------------------------------------
1080! ******************************************************************
1081! * *
1082! * LEVEL 2.5 Q2 PRODUCTION/DISSIPATION *
1083! * *
1084! ******************************************************************
1085!
1086 (ntsd,me,lmh,dtturbl,ustar,gm,gh,el,q2,epsl,epsq2,i,j,lm)
1087!----------------------------------------------------------------------
1088!
1089 IMPLICIT NONE
1090!
1091!----------------------------------------------------------------------
1092 INTEGER(KIND=KINT),INTENT(IN):: &
1093 LMH,I,J,LM,ME,NTSD
1094!
1095 REAL(KIND=kfpt),INTENT(IN):: &
1096 dtturbl,ustar
1097!
1098 REAL(KIND=kfpt),DIMENSION(1:LM-1),INTENT(IN):: &
1099 gh,gm
1100!
1101 real(kind=kfpt),dimension(1:lm-1),intent(in):: epsl
1102 real(kind=kfpt),dimension(1:lm),intent(in):: epsq2
1103!
1104 REAL(KIND=kfpt),DIMENSION(1:LM-1),INTENT(INOUT):: &
1105 el
1106!
1107 REAL(KIND=kfpt),DIMENSION(1:LM),INTENT(INOUT):: &
1108 q2
1109!----------------------------------------------------------------------
1110!***
1111!*** LOCAL VARIABLES
1112!***
1113 INTEGER(KIND=KINT):: &
1114 K
1115!,IFLAG
1116!
1117 REAL(KIND=kfpt):: &
1118 aden,aequ,anum,arhs,bden,bequ,bnum,brhs,cden,crhs &
1119 ,dloq1,eloq11,eloq12,eloq13,eloq21,eloq22,eloq31,eloq32 &
1120 ,eloq41,eloq42,eloq51,eloq52,eloqn,eqol2,ghl,gml &
1121 ,rden1,rden2,rhs2,rhsp1,rhsp2,rhst2
1122!
1123!----------------------------------------------------------------------
1124!**********************************************************************
1125!----------------------------------------------------------------------
1126!
1127
1128! IFLAG=1
1129
1130 main_integration: DO k=1,lmh-1
1131 gml=gm(k)
1132 ghl=gh(k)
1133!
1134!----------------------------------------------------------------------
1135!*** COEFFICIENTS OF THE EQUILIBRIUM EQUATION
1136!----------------------------------------------------------------------
1137!
1138 aequ=(aeqm*gml+aeqh*ghl)*ghl
1139 bequ= beqm*gml+beqh*ghl
1140!
1141!----------------------------------------------------------------------
1142!*** EQUILIBRIUM SOLUTION FOR L/Q
1143!----------------------------------------------------------------------
1144!
1145 eqol2=-0.5*bequ+sqrt(bequ*bequ*0.25-aequ)
1146!
1147!----------------------------------------------------------------------
1148!*** IS THERE PRODUCTION/DISSIPATION ?
1149!----------------------------------------------------------------------
1150!
1151 IF(((gml+ghl*ghl<=epstrb) &
1152 & .OR.(ghl>=epsgh.AND.gml/ghl<=requ) &
1153 & .OR.(eqol2<=eps2)))THEN
1154! & .OR.(EQOL2<=EPS2)).and.IFLAG.EQ.1)THEN
1155!
1156! if(ntsd.eq.23.and.me.eq.76.and.I.eq.32)then
1157! print*,'no turb=',K,GML,GHL,EPSTRB,EPSGH,REQU,EQOL2,EPS2,GML/GHL
1158! end if
1159!----------------------------------------------------------------------
1160!*** NO TURBULENCE
1161!----------------------------------------------------------------------
1162!
1163 q2(k)=epsq2(k)
1164 el(k)=epsl(k)
1165! IFLAG=2
1166!----------------------------------------------------------------------
1167!
1168 ELSE
1169!
1170!----------------------------------------------------------------------
1171!*** TURBULENCE
1172!----------------------------------------------------------------------
1173!----------------------------------------------------------------------
1174!*** COEFFICIENTS OF THE TERMS IN THE NUMERATOR
1175!----------------------------------------------------------------------
1176!
1177 anum=(anmm*gml+anmh*ghl)*ghl
1178 bnum= bnmm*gml+bnmh*ghl
1179!
1180!----------------------------------------------------------------------
1181!*** COEFFICIENTS OF THE TERMS IN THE DENOMINATOR
1182!----------------------------------------------------------------------
1183!
1184 aden=(adnm*gml+adnh*ghl)*ghl
1185 bden= bdnm*gml+bdnh*ghl
1186 cden= 1.
1187!
1188!----------------------------------------------------------------------
1189!*** COEFFICIENTS OF THE NUMERATOR OF THE LINEARIZED EQ.
1190!----------------------------------------------------------------------
1191!
1192 arhs=-(anum*bden-bnum*aden)*2.
1193 brhs=- anum*4.
1194 crhs=- bnum*2.
1195!
1196!----------------------------------------------------------------------
1197!*** INITIAL VALUE OF L/Q
1198!----------------------------------------------------------------------
1199!
1200 dloq1=el(k)/sqrt(q2(k))
1201!
1202!----------------------------------------------------------------------
1203!*** FIRST ITERATION FOR L/Q, RHS=0
1204!----------------------------------------------------------------------
1205!
1206 eloq21=1./eqol2
1207 eloq11=sqrt(eloq21)
1208 eloq31=eloq21*eloq11
1209 eloq41=eloq21*eloq21
1210 eloq51=eloq21*eloq31
1211!
1212!----------------------------------------------------------------------
1213!*** 1./DENOMINATOR
1214!----------------------------------------------------------------------
1215!
1216 rden1=1./(aden*eloq41+bden*eloq21+cden)
1217!
1218!----------------------------------------------------------------------
1219!*** D(RHS)/D(L/Q)
1220!----------------------------------------------------------------------
1221!
1222 rhsp1=(arhs*eloq51+brhs*eloq31+crhs*eloq11)*rden1*rden1
1223!
1224!----------------------------------------------------------------------
1225!*** FIRST-GUESS SOLUTION
1226!----------------------------------------------------------------------
1227!
1228 eloq12=eloq11+(dloq1-eloq11)*exp(rhsp1*dtturbl)
1229 eloq12=max(eloq12,eps1)
1230!
1231!----------------------------------------------------------------------
1232!*** SECOND ITERATION FOR L/Q
1233!----------------------------------------------------------------------
1234!
1235 eloq22=eloq12*eloq12
1236 eloq32=eloq22*eloq12
1237 eloq42=eloq22*eloq22
1238 eloq52=eloq22*eloq32
1239!
1240!----------------------------------------------------------------------
1241!*** 1./DENOMINATOR
1242!----------------------------------------------------------------------
1243!
1244 rden2=1./(aden*eloq42+bden*eloq22+cden)
1245 rhs2 =-(anum*eloq42+bnum*eloq22)*rden2+rb1
1246 rhsp2= (arhs*eloq52+brhs*eloq32+crhs*eloq12)*rden2*rden2
1247 rhst2=rhs2/rhsp2
1248!
1249!----------------------------------------------------------------------
1250!*** CORRECTED SOLUTION
1251!----------------------------------------------------------------------
1252!
1253 eloq13=eloq12-rhst2+(rhst2+dloq1-eloq12)*exp(rhsp2*dtturbl)
1254 eloq13=amax1(eloq13,eps1)
1255!
1256!----------------------------------------------------------------------
1257!*** TWO ITERATIONS IS ENOUGH IN MOST CASES ...
1258!----------------------------------------------------------------------
1259!
1260 eloqn=eloq13
1261!
1262 IF(eloqn>eps1)THEN
1263 q2(k)=el(k)*el(k)/(eloqn*eloqn)
1264 q2(k)=amax1(q2(k),epsq2(k))
1265 IF(q2(k)==epsq2(k))THEN
1266 el(k)=epsl(k)
1267 ENDIF
1268 ELSE
1269 q2(k)=epsq2(k)
1270 el(k)=epsl(k)
1271 ENDIF
1272!
1273!----------------------------------------------------------------------
1274!*** END OF TURBULENT BRANCH
1275!----------------------------------------------------------------------
1276!
1277 ENDIF
1278!----------------------------------------------------------------------
1279!*** END OF PRODUCTION/DISSIPATION LOOP
1280!----------------------------------------------------------------------
1281!
1282 ENDDO main_integration
1283!
1284!----------------------------------------------------------------------
1285!*** LOWER BOUNDARY CONDITION FOR Q2
1286!----------------------------------------------------------------------
1287!
1288 q2(lmh)=amax1(b1**(2./3.)*ustar*ustar,epsq2(lmh))
1289! if(I.eq.4)print*,'12ql test Q2(LMH)=',LMH,Q2(LMH),B1,USTAR
1290
1291!----------------------------------------------------------------------
1292!
1293 END SUBROUTINE prodq2
1294!
1295!----------------------------------------------------------------------
1296!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1297!----------------------------------------------------------------------
1298 SUBROUTINE difcof &
1299! ******************************************************************
1300! * *
1301! * LEVEL 2.5 DIFFUSION COEFFICIENTS *
1302! * *
1303! ******************************************************************
1304 (ntsd,me,lmh,lmxl,gm,gh,el,t,q2,z,akm,akh,i,j,lm,print_diag,kpbl)
1305!----------------------------------------------------------------------
1306!
1307 IMPLICIT NONE
1308!
1309!----------------------------------------------------------------------
1310 INTEGER(KIND=KINT),INTENT(IN):: &
1311 LMH,LMXL,I,J,LM,ME,NTSD,KPBL
1312!
1313 REAL(KIND=kfpt),DIMENSION(1:LM),INTENT(IN):: &
1314 q2,t
1315!
1316 REAL(KIND=kfpt),DIMENSION(1:LM-1),INTENT(IN):: &
1317 el,gh,gm
1318!
1319 REAL(KIND=kfpt),DIMENSION(1:LM+1),INTENT(IN):: &
1320 z
1321!
1322 REAL(KIND=kfpt),DIMENSION(1:LM-1),INTENT(OUT):: &
1323 akh,akm
1324!----------------------------------------------------------------------
1325!***
1326!*** LOCAL VARIABLES
1327!***
1328 INTEGER(KIND=KINT):: &
1329 K,KINV
1330!
1331 REAL(KIND=kfpt):: &
1332 aden,akmin,bden,besh,besm,cden,d2t,ell,eloq2,eloq4,elqdz &
1333 ,esh,esm,ghl,gml,q1l,rden,rdz
1334!
1335!*** Begin debugging
1336 INTEGER(KIND=KINT),INTENT(IN):: PRINT_DIAG
1337! REAL(KIND=KFPT):: D2TMIN
1338!*** End debugging
1339!
1340!----------------------------------------------------------------------
1341!**********************************************************************
1342!----------------------------------------------------------------------
1343!
1344 DO k=1,lmh-1
1345 ell=el(k)
1346!
1347 eloq2=ell*ell/q2(k)
1348 eloq4=eloq2*eloq2
1349!
1350 gml=gm(k)
1351 ghl=gh(k)
1352!
1353!----------------------------------------------------------------------
1354!*** COEFFICIENTS OF THE TERMS IN THE DENOMINATOR
1355!----------------------------------------------------------------------
1356!
1357 aden=(adnm*gml+adnh*ghl)*ghl
1358 bden= bdnm*gml+bdnh*ghl
1359 cden= 1.
1360!
1361!----------------------------------------------------------------------
1362!*** COEFFICIENTS FOR THE SM DETERMINANT
1363!----------------------------------------------------------------------
1364!
1365 besm=bsmh*ghl
1366!
1367!----------------------------------------------------------------------
1368!*** COEFFICIENTS FOR THE SH DETERMINANT
1369!----------------------------------------------------------------------
1370!
1371 besh=bshm*gml+bshh*ghl
1372!
1373!----------------------------------------------------------------------
1374!*** 1./DENOMINATOR
1375!----------------------------------------------------------------------
1376!
1377 rden=1./(aden*eloq4+bden*eloq2+cden)
1378!
1379!----------------------------------------------------------------------
1380!*** SM AND SH
1381!----------------------------------------------------------------------
1382!
1383 esm=(besm*eloq2+cesm)*rden
1384 esh=(besh*eloq2+cesh)*rden
1385!
1386!----------------------------------------------------------------------
1387!*** DIFFUSION COEFFICIENTS
1388!----------------------------------------------------------------------
1389!
1390 rdz=2./(z(k)-z(k+2))
1391 q1l=sqrt(q2(k))
1392 elqdz=ell*q1l*rdz
1393 akm(k)=elqdz*esm
1394 akh(k)=elqdz*esh
1395! if(NTSD.gt.22.and.me.eq.76.and.I.eq.32)then
1396! if(AKM(K).lt.RDZ*3.)then
1397! print*,'1K,ELQDZ,ESH,ELL,Q1L,RDZ,Q2=',K,ELQDZ,ESH &
1398! ,ELL,Q1L,RDZ,Q2(K),BESH,ELOQ2,CESH,RDEN &
1399! ,ADEN,ELOQ4,BDEN,CDEN,BSHM,GML,BSHH,GHL,BSMH &
1400! ,BDNM,BDNH,ADNM,ADNH
1401! else
1402! print*,'2K,ELQDZ,ESH,ELL,Q1L,RDZ,Q2=',K,ELQDZ,ESH &
1403! ,ELL,Q1L,RDZ,Q2(K),BESH,ELOQ2,CESH,RDEN &
1404! ,ADEN,ELOQ4,BDEN,CDEN,BSHM,GML,BSHH,GHL,BSMH &
1405! ,BDNM,BDNH,ADNM,ADNH
1406! end if
1407! if(K.eq.(LMH-1))stop
1408! end if
1409!WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1410! if(K.gt.KPBL)then
1411! AKM(K)=MAX(AKM(K),RDZ*3.)
1412! AKH(K)=MAX(AKH(K),RDZ*3.)
1413! end if
1414! AKM(K)=MAX(AKM(K),RDZ*3.)
1415! AKH(K)=MAX(AKH(K),RDZ*3.)
1416! AKM(K)=MAX(AKM(K),RDZ)
1417! AKH(K)=MAX(AKH(K),RDZ)
1418!MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
1419!----------------------------------------------------------------------
1420 ENDDO
1421! qingfu test
1422! K=LM-1
1423! RDZ=2./(Z(K)-Z(K+2))
1424! AKH(LM-1)=AKH(LM-1)*10.
1425! AKM(K)=MAX(AKM(K),RDZ*3.)*10.
1426! AKH(K)=MAX(AKH(K),RDZ*3.)*10.
1427!----------------------------------------------------------------------
1428!
1429!----------------------------------------------------------------------
1430!*** INVERSIONS
1431!----------------------------------------------------------------------
1432!
1433! IF(LMXL==LMH)THEN
1434! KINV=LMH
1435! D2TMIN=0.
1436!
1437! DO K=LMH/2,LMH-1
1438! D2T=T(K-1)-2.*T(K)+T(K+1)
1439! IF(D2T<D2TMIN)THEN
1440! D2TMIN=D2T
1441! IF(D2T<0)KINV=K
1442! ENDIF
1443! ENDDO
1444!
1445! IF(KINV<LMH)THEN
1446! DO K=KINV-1,LMH-1
1447! RDZ=2./(Z(K)-Z(K+2))
1448! AKMIN=0.5*RDZ
1449! AKM(K)=MAX(AKM(K),AKMIN)
1450! AKH(K)=MAX(AKH(K),AKMIN)
1451! ENDDO
1452!
1453!*** Begin debugging
1454! IF(PRINT_DIAG>0)THEN
1455! WRITE(6,"(A,3I3)") '{TURB1 LMXL,LMH,KINV=',LMXL,LMH,KINV
1456! WRITE(6,"(A,3I3)") '}TURB1 LMXL,LMH,KINV=',LMXL,LMH,KINV
1457! IF(PRINT_DIAG==1)THEN
1458! WRITE(6,"(A)") &
1459! '{TURB3 K, T, D2T, RDZ, Z(K), Z(K+2), AKMIN, AKH '
1460! ELSE
1461! WRITE(6,"(A)") &
1462! '}TURB3 K, T, D2T, RDZ, Z(K), Z(K+2), AKMIN, AKH '
1463! ENDIF
1464! DO K=LMH-1,KINV-1,-1
1465! D2T=T(K-1)-2.*T(K)+T(K+1)
1466! RDZ=2./(Z(K)-Z(K+2))
1467! AKMIN=0.5*RDZ
1468! IF(PRINT_DIAG==1)THEN
1469! WRITE(6,"(A,I3,F8.3,2E12.5,2F9.2,2E12.5)") '{TURB3 ' &
1470! ,K,T(K)-273.15,D2T,RDZ,Z(K),Z(K+2),AKMIN,AKH(K)
1471! ELSE
1472! WRITE(6,"(A,I3,F8.3,2E12.5,2F9.2,2E12.5)") '}TURB3 ' &
1473! ,K,T(K)-273.15,D2T,RDZ,Z(K),Z(K+2),AKMIN,AKH(K)
1474! ENDIF
1475! ENDDO
1476! ENDIF !- IF (PRINT_DIAG > 0) THEN
1477! ENDIF !- IF(KINV<LMH)THEN
1478!*** End debugging
1479!
1480! ENDIF
1481!----------------------------------------------------------------------
1482!
1483 END SUBROUTINE difcof
1484!
1485!----------------------------------------------------------------------
1486!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1487!----------------------------------------------------------------------
1488 SUBROUTINE vdifq &
1489! ******************************************************************
1490! * *
1491! * VERTICAL DIFFUSION OF Q2 (TKE) *
1492! * *
1493! ******************************************************************
1494 (lmh,dtdif,q2,el,z,i,j,lm)
1495!----------------------------------------------------------------------
1496!
1497 IMPLICIT NONE
1498!
1499!----------------------------------------------------------------------
1500 INTEGER(KIND=KINT),INTENT(IN):: &
1501 LMH,I,J,LM
1502!
1503 REAL(KIND=kfpt),INTENT(IN):: &
1504 dtdif
1505!
1506 REAL(KIND=kfpt),DIMENSION(1:LM-1),INTENT(IN):: &
1507 el
1508!
1509 REAL(KIND=kfpt),DIMENSION(1:LM+1),INTENT(IN):: &
1510 z
1511!
1512 REAL(KIND=kfpt),DIMENSION(1:LM),INTENT(INOUT):: &
1513 q2
1514!----------------------------------------------------------------------
1515!***
1516!*** LOCAL VARIABLES
1517!***
1518 INTEGER(KIND=KINT):: &
1519 K
1520!
1521 REAL(KIND=kfpt):: &
1522 aden,akqs,bden,besh,besm,cden,cf,dtozs,ell,eloq2,eloq4 &
1523 ,elqdz,esh,esm,esqhf,ghl,gml,q1l,rden,rdz
1524!
1525 REAL(KIND=kfpt),DIMENSION(1:LM-2):: &
1526 akq,cm,cr,dtoz,rsq2
1527!----------------------------------------------------------------------
1528!**********************************************************************
1529!----------------------------------------------------------------------
1530!***
1531!*** VERTICAL TURBULENT DIFFUSION
1532!***
1533!----------------------------------------------------------------------
1534 esqhf=0.5*esq
1535!
1536 DO k=1,lmh-2
1537 dtoz(k)=(dtdif+dtdif)/(z(k)-z(k+2))
1538 akq(k)=sqrt((q2(k)+q2(k+1))*0.5)*(el(k)+el(k+1))*esqhf &
1539 & /(z(k+1)-z(k+2))
1540 cr(k)=-dtoz(k)*akq(k)
1541 ENDDO
1542!
1543 cm(1)=dtoz(1)*akq(1)+1.
1544 rsq2(1)=q2(1)
1545!
1546 DO k=1+1,lmh-2
1547 cf=-dtoz(k)*akq(k-1)/cm(k-1)
1548 cm(k)=-cr(k-1)*cf+(akq(k-1)+akq(k))*dtoz(k)+1.
1549 rsq2(k)=-rsq2(k-1)*cf+q2(k)
1550 ENDDO
1551!
1552 dtozs=(dtdif+dtdif)/(z(lmh-1)-z(lmh+1))
1553 akqs=sqrt((q2(lmh-1)+q2(lmh))*0.5)*(el(lmh-1)+elz0)*esqhf &
1554 & /(z(lmh)-z(lmh+1))
1555!
1556 cf=-dtozs*akq(lmh-2)/cm(lmh-2)
1557!
1558 q2(lmh-1)=(dtozs*akqs*q2(lmh)-rsq2(lmh-2)*cf+q2(lmh-1)) &
1559 & /((akq(lmh-2)+akqs)*dtozs-cr(lmh-2)*cf+1.)
1560!
1561 DO k=lmh-2,1,-1
1562 q2(k)=(-cr(k)*q2(k+1)+rsq2(k))/cm(k)
1563 ENDDO
1564!----------------------------------------------------------------------
1565!
1566 END SUBROUTINE vdifq
1567!
1568!----------------------------------------------------------------------
1569!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1570!---------------------------------------------------------------------
1571 SUBROUTINE vdifh(DTDIF,LMH,THZ0,QZ0,RKHS,CHKLOWQ,CT &
1572 ,TH,Q,CWM,RKH,Z,RHO,I,J,LM)
1573! ***************************************************************
1574! * *
1575! * VERTICAL DIFFUSION OF MASS VARIABLES *
1576! * *
1577! ***************************************************************
1578!---------------------------------------------------------------------
1579!
1580 IMPLICIT NONE
1581!
1582!---------------------------------------------------------------------
1583 INTEGER(KIND=KINT),INTENT(IN):: &
1584 LMH,I,J,LM
1585!
1586 REAL(KIND=kfpt),INTENT(IN):: &
1587 chklowq,ct,dtdif,qz0,rkhs,thz0
1588!
1589 REAL(KIND=kfpt),DIMENSION(1:LM-1),INTENT(IN):: &
1590 rkh
1591!
1592 REAL(KIND=kfpt),DIMENSION(1:LM),INTENT(IN):: &
1593 rho
1594!
1595 REAL(KIND=kfpt),DIMENSION(1:LM+1),INTENT(IN):: &
1596 z
1597!
1598 REAL(KIND=kfpt),DIMENSION(1:LM),INTENT(INOUT):: &
1599 cwm,q,th
1600!
1601!----------------------------------------------------------------------
1602!***
1603!*** LOCAL VARIABLES
1604!***
1605 INTEGER(KIND=KINT):: &
1606 K
1607!
1608 REAL(KIND=kfpt):: &
1609 cf,cmb,cmcb,cmqb,cmtb,cthf,dtozl,dtozs &
1610 ,rcml,rkhh,rkqs,rscb,rsqb,rstb &
1611 ,tem1
1612!
1613 REAL(KIND=kfpt),DIMENSION(1:LM-1):: &
1614 cm,cr,dtoz,rkct,rsc,rsq,rst
1615!
1616!----------------------------------------------------------------------
1617!**********************************************************************
1618!----------------------------------------------------------------------
1619 cthf=0.5*ct
1620!
1621 DO k=1,lmh-1
1622 dtoz(k)=dtdif/(z(k)-z(k+1))
1623 cr(k)=-dtoz(k)*rkh(k)
1624 rkct(k)=rkh(k)*(z(k)-z(k+2))*cthf
1625 ENDDO
1626!
1627 cm(1)=dtoz(1)*rkh(1)+rho(1)
1628!----------------------------------------------------------------------
1629 rst(1)=th(1)*rho(1)-rkct(1)*dtoz(1)
1630 rsq(1)=q(1)*rho(1)
1631 rsc(1)=cwm(1)*rho(1)
1632!----------------------------------------------------------------------
1633 DO k=1+1,lmh-1
1634 dtozl=dtoz(k)
1635 cf=-dtozl*rkh(k-1)/cm(k-1)
1636 cm(k)=-cr(k-1)*cf+(rkh(k-1)+rkh(k))*dtozl+rho(k)
1637 rst(k)=-rst(k-1)*cf+(rkct(k-1)-rkct(k))*dtozl+th(k)*rho(k)
1638 rsq(k)=-rsq(k-1)*cf+q(k) *rho(k)
1639 rsc(k)=-rsc(k-1)*cf+cwm(k)*rho(k)
1640 ENDDO
1641!
1642 dtozs=dtdif/(z(lmh)-z(lmh+1))
1643 rkhh=rkh(lmh-1)
1644!
1645 cf=-dtozs*rkhh/cm(lmh-1)
1646 rkqs=rkhs*chklowq
1647!
1648 cmb=cr(lmh-1)*cf
1649 cmtb=-cmb+(rkhh+rkhs)*dtozs+rho(lmh)
1650 cmqb=-cmb+(rkhh+rkqs)*dtozs+rho(lmh)
1651 cmcb=-cmb+(rkhh )*dtozs+rho(lmh)
1652!
1653 rstb=-rst(lmh-1)*cf+rkct(lmh-1)*dtozs+th(lmh)*rho(lmh)
1654 rsqb=-rsq(lmh-1)*cf+q(lmh) *rho(lmh)
1655 rscb=-rsc(lmh-1)*cf+cwm(lmh)*rho(lmh)
1656!----------------------------------------------------------------------
1657 tem1=th(lmh)
1658 th(lmh) =(dtozs*rkhs*thz0+rstb)/cmtb
1659 q(lmh) =(dtozs*rkqs*qz0 +rsqb)/cmqb
1660 cwm(lmh)=( rscb)/cmcb
1661 if(abs(tem1-th(lmh)).gt.8.)then
1662! print*,'TH(LMH),DTOZS,RKHS,THZ0,RSTB,CMTB=', &
1663! DTOZS*RKHS*THZ0,RSTB,CMTB, &
1664! TH(LMH),tem1,thz0,DTOZS,RKHS, &
1665! CMB,RKHH,RKHS,RHO(LMH),CR(LMH-1),CF,RKH(LMH-1), &
1666! DTOZS*RKQS*QZ0,RSQB,CMQB
1667 end if
1668!----------------------------------------------------------------------
1669 DO k=lmh-1,1,-1
1670 rcml=1./cm(k)
1671 th(k) =(-cr(k)* th(k+1)+rst(k))*rcml
1672 q(k) =(-cr(k)* q(k+1)+rsq(k))*rcml
1673 cwm(k)=(-cr(k)*cwm(k+1)+rsc(k))*rcml
1674 ENDDO
1675!----------------------------------------------------------------------
1676!
1677 END SUBROUTINE vdifh
1678!
1679!---------------------------------------------------------------------
1680!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1681!---------------------------------------------------------------------
1682 SUBROUTINE vdifv(LMH,DTDIF,UZ0,VZ0,RKMS,D,U,V,RKM,Z,RHO,I,J,LM)
1683! ***************************************************************
1684! * *
1685! * VERTICAL DIFFUSION OF VELOCITY COMPONENTS *
1686! * *
1687! ***************************************************************
1688!---------------------------------------------------------------------
1689!
1690 IMPLICIT NONE
1691!
1692!---------------------------------------------------------------------
1693 INTEGER(KIND=KINT),INTENT(IN):: &
1694 LMH,I,J,LM
1695!
1696 REAL(KIND=kfpt),INTENT(IN):: &
1697 rkms,dtdif,uz0,vz0
1698!
1699 REAL(KIND=kfpt),DIMENSION(1:LM-1),INTENT(IN):: &
1700 d,rkm
1701!
1702 REAL(KIND=kfpt),DIMENSION(1:LM),INTENT(IN):: &
1703 rho
1704!
1705 REAL(KIND=kfpt),DIMENSION(1:LM+1),INTENT(IN):: &
1706 z
1707!
1708 REAL(KIND=kfpt),DIMENSION(1:LM),INTENT(INOUT):: &
1709 u,v
1710!----------------------------------------------------------------------
1711!***
1712!*** LOCAL VARIABLES
1713!***
1714 INTEGER(KIND=KINT):: &
1715 K
1716!
1717 REAL(KIND=kfpt):: &
1718 cf,dtozak,dtozl,dtozs,rcml,rcmvb,rhok,rkmh
1719!
1720 REAL(KIND=kfpt),DIMENSION(1:LM-1):: &
1721 cm,cr,dtoz,rsu,rsv
1722!----------------------------------------------------------------------
1723!**********************************************************************
1724!----------------------------------------------------------------------
1725 DO k=1,lmh-1
1726 dtoz(k)=dtdif/(z(k)-z(k+1))
1727 cr(k)=-dtoz(k)*rkm(k)
1728 ENDDO
1729!
1730 rhok=rho(1)
1731 cm(1)=dtoz(1)*rkm(1)+rhok
1732 rsu(1)=u(1)*rhok
1733 rsv(1)=v(1)*rhok
1734!----------------------------------------------------------------------
1735 DO k=2,lmh-1
1736 dtozl=dtoz(k)
1737 cf=-dtozl*rkm(k-1)/cm(k-1)
1738 rhok=rho(k)
1739 cm(k)=-cr(k-1)*cf+(rkm(k-1)+rkm(k)+d(k)*rkms)*dtozl+rhok
1740 rsu(k)=-rsu(k-1)*cf+u(k)*rhok
1741 rsv(k)=-rsv(k-1)*cf+v(k)*rhok
1742 ENDDO
1743!----------------------------------------------------------------------
1744 dtozs=dtdif/(z(lmh)-z(lmh+1))
1745 rkmh=rkm(lmh-1)
1746!
1747 cf=-dtozs*rkmh/cm(lmh-1)
1748 rhok=rho(lmh)
1749 rcmvb=1./((rkmh+rkms)*dtozs-cr(lmh-1)*cf+rhok)
1750 dtozak=dtozs*rkms
1751!----------------------------------------------------------------------
1752 u(lmh)=(dtozak*uz0-rsu(lmh-1)*cf+u(lmh)*rhok)*rcmvb
1753 v(lmh)=(dtozak*vz0-rsv(lmh-1)*cf+v(lmh)*rhok)*rcmvb
1754!----------------------------------------------------------------------
1755 DO k=lmh-1,1,-1
1756 rcml=1./cm(k)
1757 u(k)=(-cr(k)*u(k+1)+rsu(k))*rcml
1758 v(k)=(-cr(k)*v(k+1)+rsv(k))*rcml
1759 ENDDO
1760!----------------------------------------------------------------------
1761!
1762 END SUBROUTINE vdifv
1763!
1764!-----------------------------------------------------------------------
1765!
1766!=======================================================================
1767! SUBROUTINE MYJPBL_INIT(EXCH_H &
1768 SUBROUTINE myjpbl_init(IDS,IDE,JDS,JDE,LM &
1769 & ,IMS,IME,JMS,JME &
1770 & ,ITS,ITE,JTS,JTE )
1771!-----------------------------------------------------------------------
1772 IMPLICIT NONE
1773!-----------------------------------------------------------------------
1774! LOGICAL(KIND=KLOG),INTENT(IN):: &
1775! RESTART
1776
1777 INTEGER(KIND=KINT),INTENT(IN):: &
1778 ids,ide,jds,jde,lm &
1779 ,ims,ime,jms,jme &
1780 ,its,ite,jts,jte
1781
1782! REAL(KIND=KFPT),DIMENSION(IMS:IME,JMS:JME,1:LM),INTENT(OUT):: &
1783! EXCH_H
1784!
1785 INTEGER(KIND=KINT):: &
1786 i,j,k,itf,jtf,ktf
1787!-----------------------------------------------------------------------
1788!-----------------------------------------------------------------------
1789!
1790 jtf=min0(jte,jde-1)
1791 ktf=lm
1792 itf=min0(ite,ide-1)
1793!
1794! IF(.NOT.RESTART)THEN
1795! DO K=1,KTF
1796! DO J=JTS,JTF
1797! DO I=ITS,ITF
1798! EXCH_H(I,J,K)=0.
1799! ENDDO
1800! ENDDO
1801! ENDDO
1802! ENDIF
1803!
1804!-----------------------------------------------------------------------
1805 CALL herftbl
1806!-----------------------------------------------------------------------
1807 CALL tablept
1808 CALL tablett
1809!-----------------------------------------------------------------------
1810!
1811 END SUBROUTINE myjpbl_init
1812!
1813!-----------------------------------------------------------------------
1814!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1815 SUBROUTINE herftbl
1816! ******************************************************************
1817! * *
1818! * POSITIVE HALF OF ERF FUNCTION TABLE *
1819! * RESPONSIBLE PERSON: Z.JANJIC *
1820! * *
1821! ******************************************************************
1822!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1823IMPLICIT NONE
1824
1825!--LOCAL VARIABLES------------------------------------------------------
1826INTEGER(KIND=KINT):: &
1827 K ! INDEX
1828
1829REAL(KIND=kfpt):: &
1830 dxe & ! ARGUMENT INCREMENT
1831,dxh & !
1832,rden & !
1833,x & ! ARGUMENT
1834,xgmin !
1835
1836REAL(KIND=kfpt),DIMENSION(1:KERFM+1):: &
1837 gauss !
1838
1839!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1840!-----------------------------------------------------------------------
1841 dxe=(xemax-xemin)/float(kerfm-1)
1842 dxh=dxe*0.5
1843 xgmin=xemin-dxh
1844!
1845 rdxe=1./dxe
1846!
1847 DO k=1,kerfm+1
1848 x=(k-1)*dxe+xgmin
1849 gauss(k)=exp(-x*x*0.5)
1850 ENDDO
1851!
1852 rden=1./sqrt(pi+pi)
1853 herff(1)=0.
1854 DO k=2,kerfm
1855 herff(k)=((gauss(k)+gauss(k+1))*dxh)*rden+herff(k-1)
1856 ENDDO
1857!
1858 DO k=1,kerfm
1859 herff(k)=0.5-herff(k)
1860 ENDDO
1861!-----------------------------------------------------------------------
1862 ENDSUBROUTINE herftbl
1863!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1864 FUNCTION herf(X)
1865! ******************************************************************
1866! * *
1867! * ERF FUNCTION HALF-TABLE *
1868! * RESPONSIBLE PERSON: Z.JANJIC *
1869! * *
1870! ******************************************************************
1871!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1872IMPLICIT NONE
1873
1874!-----------------------------------------------------------------------
1875REAL(kind=kfpt):: &
1876 herf
1877
1878!--LOCAL VARIABLES------------------------------------------------------
1879INTEGER(KIND=KINT):: &
1880 k ! INDEX
1881
1882REAL(kind=kfpt):: &
1883 ak & ! POSITION IN TABLE
1884,x ! ARGUMENT
1885!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1886!-----------------------------------------------------------------------
1887 ak=(x-xemin)*rdxe
1888 k=int(ak)
1889 k=max(k,0)
1890 k=min(k,kerfm2)
1891!
1892 herf=(herff(k+2)-herff(k+1))*(ak-real(k))+herff(k+1)
1893!-----------------------------------------------------------------------
1894 ENDFUNCTION herf
1895!-----------------------------------------------------------------------
1896!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1897 SUBROUTINE tablept
1898! ******************************************************************
1899! * *
1900! * GENERATES THE TABLE FOR FINDING PRESSURE FROM *
1901! * SATURATION SPECIFIC HUMIDITY AND POTENTIAL TEMPERATURE *
1902! * *
1903! ******************************************************************
1904!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1905IMPLICIT NONE
1906
1907!-----------------------------------------------------------------------
1908REAL(KIND=kfpt),PARAMETER:: &
1909 eps=1.e-10
1910!--LOCAL VARIABLES------------------------------------------------------
1911INTEGER(KIND=KINT):: &
1912 KTH & ! INDEX
1913,KP ! INDEX
1914
1915REAL(KIND=kfpt):: &
1916 rxner & ! 1./EXNER FUNCTION
1917,dth & ! POTENTIAL TEMPERATURE STEP
1918,dp & ! PRESSURE STEP
1919,dqs & ! SATURATION SPECIFIC HUMIDITY STEP
1920,p & ! PRESSURE
1921,qs0k & ! BASE VALUE FOR SATURATION HUMIDITY
1922,sqsk & ! SATURATION SPEC. HUMIDITY RANGE
1923,th ! POTENTIAL TEMPERATURE
1924
1925REAL(KIND=kfpt),DIMENSION(1:ITBL):: &
1926 app & ! TEMPORARY
1927,aqp & ! TEMPORARY
1928,pnew & ! NEW PRESSURES
1929,pold & ! OLD PRESSURE
1930,qsnew & ! NEW SATURATION SPEC. HUMIDITY
1931,qsold & ! OLD SATURATION SPEC. HUMIDITY
1932,y2p ! TEMPORARY
1933!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1934!-----------------------------------------------------------------------
1935 dth=(thh-thl)/real(jtbl-1)
1936 dp=(ph-pl)/real(itbl-1)
1937 rdth=1./dth
1938!-----------------------------------------------------------------------
1939 th=thl-dth
1940 DO kth=1,jtbl
1941 th=th+dth
1942 p=pl-dp
1943 DO kp=1,itbl
1944 p=p+dp
1945 rxner=(100000./p)**cappa
1946 qsold(kp)=pq0/p*exp(a2*(th-a3*rxner)/(th-a4*rxner))
1947 pold(kp)=p
1948 ENDDO
1949!
1950 qs0k=qsold(1)
1951 sqsk=qsold(itbl)-qsold(1)
1952 qsold(1)=0.
1953 qsold(itbl)=1.
1954!
1955 DO kp=2,itbl-1
1956 qsold(kp)=(qsold(kp)-qs0k)/sqsk
1957!WWWWWWWWWWWWWW FIX DUE TO 32 BIT PRECISION LIMITATION WWWWWWWWWWWWWWWWW
1958 IF((qsold(kp)-qsold(kp-1)).LT.eps) THEN
1959 qsold(kp)=qsold(kp-1)+eps
1960 ENDIF
1961!MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
1962 ENDDO
1963!
1964 qs0(kth)=qs0k
1965 sqs(kth)=sqsk
1966!
1967 qsnew(1)=0.
1968 qsnew(itbl)=1.
1969 dqs=1./real(itbl-1)
1970 rdq=1./dqs
1971!
1972 DO kp=2,itbl-1
1973 qsnew(kp)=qsnew(kp-1)+dqs
1974 ENDDO
1975!
1976 y2p(1)=0.
1977 y2p(itbl)=0.
1978!
1979 CALL spline(jtbl,itbl,qsold,pold,y2p,itbl,qsnew,pnew,app,aqp)
1980!
1981 DO kp=1,itbl
1982 ptbl(kp,kth)=pnew(kp)
1983 ENDDO
1984!-----------------------------------------------------------------------
1985 ENDDO
1986!-----------------------------------------------------------------------
1987 ENDSUBROUTINE tablept
1988!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1989 SUBROUTINE tablett
1990! ******************************************************************
1991! * *
1992! * GENERATES THE TABLE FOR FINDING TEMPERATURE FROM *
1993! * PRESSURE AND EQUIVALENT POTENTIAL TEMPERATURE *
1994! * *
1995! ******************************************************************
1996!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1997IMPLICIT NONE
1998
1999!-----------------------------------------------------------------------
2000REAL(KIND=kfpt),PARAMETER:: &
2001 eps=1.e-10
2002!--LOCAL VARIABLES------------------------------------------------------
2003INTEGER(KIND=KINT):: &
2004 KTH & ! INDEX
2005,KP ! INDEX
2006
2007REAL(KIND=kfpt):: &
2008 rxner & ! 1./EXNER FUNCTION
2009,dth & ! POTENTIAL TEMPERATURE STEP
2010,dp & ! PRESSURE STEP
2011,dthe & ! EQUIVALENT POT. TEMPERATURE STEP
2012,p & ! PRESSURE
2013,qs & ! SATURATION SPECIFIC HUMIDITY
2014,the0k & ! BASE VALUE FOR EQUIVALENT POT. TEMPERATURE
2015,sthek & ! EQUIVALENT POT. TEMPERATURE RANGE
2016,th ! POTENTIAL TEMPERATURE
2017
2018REAL(KIND=kfpt),DIMENSION(1:JTBL):: &
2019 apt & ! TEMPORARY
2020,aqt & ! TEMPORARY
2021,tnew & ! NEW TEMPERATURE
2022,told & ! OLD TEMPERATURE
2023,thenew & ! NEW EQUIVALENT POTENTIAL TEMPERATURE
2024,theold & ! OLD EQUIVALENT POTENTIAL TEMPERATURE
2025,y2t ! TEMPORARY
2026!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2027!-----------------------------------------------------------------------
2028 dth=(thh-thl)/real(jtbl-1)
2029 dp=(ph-pl)/real(itbl-1)
2030 rdp=1./dp
2031!-----------------------------------------------------------------------
2032 p=pl-dp
2033 DO kp=1,itbl
2034 p=p+dp
2035 th=thl-dth
2036 DO kth=1,jtbl
2037 th=th+dth
2038 rxner=(100000./p)**cappa
2039 qs=pq0/p*exp(a2*(th-a3*rxner)/(th-a4*rxner))
2040 told(kth)=th/rxner
2041 theold(kth)=th*exp(eliwv*qs/(cp*told(kth)))
2042 ENDDO
2043!
2044 the0k=theold(1)
2045 sthek=theold(jtbl)-theold(1)
2046 theold(1)=0.
2047 theold(jtbl)=1.
2048!
2049 DO kth=2,jtbl-1
2050 theold(kth)=(theold(kth)-the0k)/sthek
2051!WWWWWWWWWWWWWW FIX DUE TO 32 BIT PRECISION LIMITATION WWWWWWWWWWWWWWWWW
2052 IF((theold(kth)-theold(kth-1)).LT.eps) THEN
2053 theold(kth)=theold(kth-1)+eps
2054 ENDIF
2055!MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
2056 ENDDO
2057!
2058 the0(kp)=the0k
2059 sthe(kp)=sthek
2060!
2061 thenew(1)=0.
2062 thenew(jtbl)=1.
2063 dthe=1./real(jtbl-1)
2064 rdthe=1./dthe
2065!
2066 DO kth=2,jtbl-1
2067 thenew(kth)=thenew(kth-1)+dthe
2068 ENDDO
2069!
2070 y2t(1)=0.
2071 y2t(jtbl)=0.
2072!
2073 CALL spline(jtbl,jtbl,theold,told,y2t,jtbl,thenew,tnew,apt,aqt)
2074!
2075 DO kth=1,jtbl
2076 ttbl(kth,kp)=tnew(kth)
2077 ENDDO
2078!-----------------------------------------------------------------------
2079 ENDDO
2080!-----------------------------------------------------------------------
2081 ENDSUBROUTINE tablett
2082!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2083 SUBROUTINE spline(JTBL,NOLD,XOLD,YOLD,Y2,NNEW,XNEW,YNEW,P,Q)
2084! ******************************************************************
2085! * *
2086! * THIS IS A ONE-DIMENSIONAL CUBIC SPLINE FITTING ROUTINE *
2087! * PROGRAMED FOR A SMALL SCALAR MACHINE. *
2088! * *
2089! * PROGRAMER: Z. JANJIC, YUGOSLAV FED. HYDROMET. INST., BEOGRAD *
2090! * *
2091! * *
2092! * *
2093! * NOLD - NUMBER OF GIVEN VALUES OF THE FUNCTION. MUST BE GE 3. *
2094! * XOLD - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
2095! * FUNCTION ARE GIVEN. MUST BE IN ASCENDING ORDER. *
2096! * YOLD - THE GIVEN VALUES OF THE FUNCTION AT THE POINTS XOLD. *
2097! * Y2 - THE SECOND DERIVATIVES AT THE POINTS XOLD. IF NATURAL *
2098! * SPLINE IS FITTED Y2(1)=0. AND Y2(NOLD)=0. MUST BE *
2099! * SPECIFIED. *
2100! * NNEW - NUMBER OF VALUES OF THE FUNCTION TO BE CALCULATED. *
2101! * XNEW - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
2102! * FUNCTION ARE CALCULATED. XNEW(K) MUST BE GE XOLD(1) *
2103! * AND LE XOLD(NOLD). *
2104! * YNEW - THE VALUES OF THE FUNCTION TO BE CALCULATED. *
2105! * P, Q - AUXILIARY VECTORS OF THE LENGTH NOLD-2. *
2106! * *
2107! ******************************************************************
2108 IMPLICIT REAL(KIND=KFPT)(a-h,o-z),INTEGER(KIND=KINT)(I-N)
2109!-----------------------------------------------------------------------
2110 dimension &
2111 xold(jtbl),yold(jtbl),y2(jtbl),p(jtbl),q(jtbl) &
2112 ,xnew(jtbl),ynew(jtbl)
2113!-----------------------------------------------------------------------
2114 noldm1=nold-1
2115!
2116 dxl=xold(2)-xold(1)
2117 dxr=xold(3)-xold(2)
2118 dydxl=(yold(2)-yold(1))/dxl
2119 dydxr=(yold(3)-yold(2))/dxr
2120 rtdxc=.5/(dxl+dxr)
2121!
2122 p(1)= rtdxc*(6.*(dydxr-dydxl)-dxl*y2(1))
2123 q(1)=-rtdxc*dxr
2124!
2125 IF(nold.EQ.3) GO TO 700
2126!-----------------------------------------------------------------------
2127 k=3
2128!
2129 100 dxl=dxr
2130 dydxl=dydxr
2131 dxr=xold(k+1)-xold(k)
2132 dydxr=(yold(k+1)-yold(k))/dxr
2133 dxc=dxl+dxr
2134 den=1./(dxl*q(k-2)+dxc+dxc)
2135!
2136 p(k-1)= den*(6.*(dydxr-dydxl)-dxl*p(k-2))
2137 q(k-1)=-den*dxr
2138!
2139 k=k+1
2140 IF(k.LT.nold) GO TO 100
2141!-----------------------------------------------------------------------
2142 700 k=noldm1
2143!
2144 200 y2(k)=p(k-1)+q(k-1)*y2(k+1)
2145!
2146 k=k-1
2147 IF(k.GT.1) GO TO 200
2148!-----------------------------------------------------------------------
2149 k1=1
2150!
2151 300 xk=xnew(k1)
2152!
2153 DO 400 k2=2,nold
2154 IF(xold(k2).LE.xk) GO TO 400
2155 kold=k2-1
2156 GO TO 450
2157 400 CONTINUE
2158 ynew(k1)=yold(nold)
2159 GO TO 600
2160!
2161 450 IF(k1.EQ.1) GO TO 500
2162 IF(k.EQ.kold) GO TO 550
2163!
2164 500 k=kold
2165!
2166 y2k=y2(k)
2167 y2kp1=y2(k+1)
2168 dx=xold(k+1)-xold(k)
2169 rdx=1./dx
2170!
2171 ak=.1666667*rdx*(y2kp1-y2k)
2172 bk=.5*y2k
2173 ck=rdx*(yold(k+1)-yold(k))-.1666667*dx*(y2kp1+y2k+y2k)
2174!
2175 550 x=xk-xold(k)
2176 xsq=x*x
2177!
2178 ynew(k1)=ak*xsq*x+bk*xsq+ck*x+yold(k)
2179!
2180 600 k1=k1+1
2181 IF(k1.LE.nnew) GO TO 300
2182!-----------------------------------------------------------------------
2183 ENDSUBROUTINE spline
2184!-----------------------------------------------------------------------
2185!
2186 END MODULE module_bl_myjpbl
2187!
2188!-----------------------------------------------------------------------