UPP (develop)
Loading...
Searching...
No Matches
INITPOST_NEMS.f
Go to the documentation of this file.
1
20!----------------------------------------------------------------------
26!----------------------------------------------------------------------
27 SUBROUTINE initpost_nems(NREC,nfile)
28
29 use vrbls3d, only: t, q, uh, vh, q2, cwm, f_ice, f_rain, f_rimef, cfr, pint,&
30 pint, alpint, pmid, pmidv, zint, zmid, wh, rlwtt, rswtt,&
31 ttnd, tcucn, train, el_pbl, exch_h, omga, qqni, qqnr, qqw, qqi, &
32 qqr, qqs, qqg, ref_10cm, radius_cloud, radius_ice, radius_snow, &
33 extcof55, aextc55
34 use vrbls2d, only: f, pd, fis, pblh, mixht, ustar, z0, ths, qs, twbs, qwbs, prec,&
35 acprec, cuprec,ancprc, lspa, sno, snoavg, psfcavg, t10avg, t10m, akhsavg, akmsavg,&
36 refd_max, w_up_max, w_dn_max, up_heli_max, si, cldefi, th10, q10, pshltr,&
37 tshltr, qshltr, maxtshltr, mintshltr, maxrhshltr, minrhshltr, akhs, akms, albase,&
38 albedo, czen, cfracl, cfracm, islope, cmc, grnflx, pctsno, soiltb, vegfrc,&
39 acfrcv, acfrst, ssroff, bgroff, czmean, mxsnal, radot, sigt4, tg, sr, cfrach,&
40 rlwin, rlwtoa, alwin, alwout, alwtoa, rswin, rswinc, rswout, aswin,aswout,&
41 aswtoa, sfcshx, sfclhx, subshx, snopcx, sfcuvx, potevp, ncfrcv, ncfrst, u10h,&
42 u10, v10h, v10, u10max, v10max, smstav, smstot, sfcevp, ivgtyp, acsnow, acsnom,&
43 sst, thz0, qz0, uz0, vz0, htop, isltyp, sfcexc, hbot, htopd, htops, cuppt, cprate,&
44 hbotd, hbots, prate_max, fprate_max
45 use soil, only: sldpth, sh2o, smc, stc
46 use masks, only: lmv, lmh, htm, vtm, dx, dy, hbm2, gdlat, gdlon, sm, sice
47 use kinds, only: i_llong
48 use wrf_io_flags_mod, only:
49 use params_mod, only: pi, dtr, g, d608, rd,tfrz
50 use lookup_mod, only: thl, plq, ptbl, ttbl, rdq, rdth, rdp, rdthe, pl, qs0, sqs, sthe, the0,&
51 ttblq, rdpq, rdtheq, stheq, the0q
52 use ctlblk_mod, only: me, mpi_comm_comp, global, icnt, idsp, jsta, ihrst, imin, idat, sdat,&
53 ifhr, ifmin, filename, restrt, imp_physics, isf_surface_physics, icu_physics, jend,&
54 dt, spval, gdsdegr, grib, pdtop, pt, tmaxmin, nsoil, lp1, jend_m, nprec, nphs, avrain,&
55 avcnvc, ardlw, ardsw, asrfc, novegtype, spl, lsm, dtq2, tsrfc, trdlw, trdsw, theat, tclod,&
56 tprec, alsl, lm , im, jm, jsta_2l, jend_2u, ivegsrc, pthresh
57 use gridspec_mod, only: dyval, dxval, cenlat, cenlon, maptype, gridtype, latstart, latlast, latnw,&
58 latse, lonstart, lonlast, lonnw, lonse, latstartv, latlastv, cenlatv, lonstartv,&
59 lonlastv, cenlonv
60 use nemsio_module, only: nemsio_gfile, nemsio_getfilehead, nemsio_close, nemsio_getheadvar
61 use upp_math, only: h2u
62!
63! INCLUDE/SET PARAMETERS.
64 implicit none
65!
66 type(nemsio_gfile),intent(inout) :: nfile
67!
68 include "mpif.h"
69! This version of INITPOST shows how to initialize, open, read from, and
70! close a NetCDF dataset. In order to change it to read an internal (binary)
71! dataset, do a global replacement of _ncd_ with _int_.
72
73 character(len=20) :: VarName
74 character(len=20) :: VcoordName
75 integer :: Status
76 character startdate*19,SysDepInfo*80,cgar*1
77 character startdate2(19)*4
78!
79! NOTE: SOME INTEGER VARIABLES ARE READ INTO DUMMY ( A REAL ). THIS IS OK
80! AS LONG AS REALS AND INTEGERS ARE THE SAME SIZE.
81!
82! ALSO, EXTRACT IS CALLED WITH DUMMY ( A REAL ) EVEN WHEN THE NUMBERS ARE
83! INTEGERS - THIS IS OK AS LONG AS INTEGERS AND REALS ARE THE SAME SIZE.
84 logical, parameter :: debugprint = .false.
85! logical, parameter :: debugprint = .true.
86 logical :: convert_rad_to_deg=.false.
87! logical global
88! CHARACTER BLANK*4
89 integer nfhour ! forecast hour from nems io file
90 INTEGER IDATE(8),JDATE(8)
91!
92! DECLARE VARIABLES.
93!
94 REAL FACT,tsph,tstart
95 REAL RINC(5)
96 REAL ETA1(LM+1), ETA2(LM+1)
97 REAL GARB
98 REAL DUM1D (LM+1)
99 REAL DUMMY ( IM, JM )
100 REAL DUMMY2 ( IM, JM )
101 real, allocatable :: fi(:,:,:)
102
103 integer ibuf(im,jsta_2l:jend_2u)
104 real buf(im,jsta_2l:jend_2u)
105 character*8,allocatable:: recname(:)
106 character*16,allocatable :: reclevtyp(:)
107 integer,allocatable:: reclev(:)
108 real, allocatable:: bufy(:)
109 real, allocatable:: glat1d(:),glon1d(:)
110!jw
111 integer ii,jj,js,je,jev,iyear,imn,iday,itmp,ioutcount,istatus, &
112 nsrfc,nrdlw,nrdsw,nheat,nclod, &
113 iunit,nrec,i,j,l, iret,nframe,impf,jmpf,nframed2, &
114 igdout,ll,n,im1,jm1,iim1,item
115!
116! DATA BLANK/' '/
117!
118!***********************************************************************
119! START INIT HERE.
120!
121 if(me==0)WRITE(6,*)'INITPOST: ENTER INITPOST'
122!
123!
124! STEP 1. READ MODEL OUTPUT FILE
125!
126!***
127! LMH always = LM for sigma-type vert coord
128! LMV always = LM for sigma-type vert coord
129
130 do j = jsta_2l, jend_2u
131 do i = 1, im
132 lmv( i, j ) = lm
133 lmh( i, j ) = lm
134 end do
135 end do
136
137! HTM VTM all 1 for sigma-type vert coord
138
139 do l = 1, lm
140 do j = jsta_2l, jend_2u
141 do i = 1, im
142 htm( i, j, l ) = 1.0
143 vtm( i, j, l ) = 1.0
144 end do
145 end do
146 end do
147
148! The end j row is going to be jend_2u for all variables except for V.
149 js=jsta_2l
150 je=jend_2u
151 IF (jend_2u==jm) THEN
152 jev=jend_2u+1
153 ELSE
154 jev=jend_2u
155 ENDIF
156! sample print point
157 ii=(1+im)/2
158 jj=(1+jm)/2
159! get start date
160 idate=0
161 if (me == 0)then
162 print*,'nrec=',nrec
163 allocate(recname(nrec),reclevtyp(nrec),reclev(nrec))
164
165 call nemsio_getfilehead(nfile,iret=iret &
166 ,idate=idate(1:7),nfhour=nfhour,recname=recname &
167 ,reclevtyp=reclevtyp,reclev=reclev,nframe=nframe)
168! if(iret/=0)print*,'error getting idate,fhour, stopping';stop
169 if (debugprint) then
170 print *,'printing an inventory of NEMS file'
171 do i=1,nrec
172 print *,'recname,reclevtyp,reclev=',trim(recname(i)),' ', &
173 trim(reclevtyp(i)),reclev(i)
174 end do
175 endif
176! print *,'reclevtyp=',(trim(reclevtyp(i)),i=1,nrec)
177! print *,'reclev=',(reclev(i),i=1,nrec)
178 deallocate(recname,reclevtyp,reclev)
179 impf=im+nframe*2
180 jmpf=jm+nframe*2
181! nframed2=nframe/2
182 print*,'nframe,impf,jmpf= ',nframe,impf,jmpf
183 allocate(glat1d(impf*jmpf),glon1d(impf*jmpf) )
184 call nemsio_getfilehead(nfile,dx=glat1d &
185 ,dy=glon1d,iret=iret)
186 if(iret/=0)print*,'did not find dx dy'
187! do j=1,impf*jmpf
188! print*,'dx before scatter= ',j,glat1d(j)
189! end do
190!$omp parallel do private(i,j,item)
191 do j=1,jm
192 item = (j-1)*impf + nframe
193 do i=1,im
194 dummy(i,j) = glat1d(item+i)
195 dummy2(i,j) = glon1d(item+i)
196! dummy(i,j)=glat1d(i-nframe,j-nframe)
197! dummy2(i,j)=glon1d(i-nframe,j-nframe)
198 end do
199 end do
200 deallocate(glat1d,glon1d)
201! latstart=nint(dummy(1,1)*1000.)
202! latlast=nint(dummy(im,jm)*1000.)
203! lonstart=nint(dummy2(1,1)*1000.)
204! lonlast=nint(dummy2(im,jm)*1000.)
205! dyval=nint((dummy(1,2)-dummy(1,1))*1000.)
206! dxval=nint((dummy(2,1)-dummy(1,1))*1000.)
207! cenlat=nint(dummy(ii,jj)*1000.)
208! cenlon=nint(dummy2(ii,jj)*1000.)
209 print*,'idate before broadcast = ',(idate(i),i=1,7)
210 end if
211 call mpi_bcast(idate(1),7,mpi_integer,0,mpi_comm_comp,iret)
212 call mpi_bcast(nfhour,1,mpi_integer,0,mpi_comm_comp,iret)
213 call mpi_bcast(nframe,1,mpi_integer,0,mpi_comm_comp,iret)
214! call mpi_bcast(latstart,1,MPI_INTEGER,0,mpi_comm_comp,iret)
215! call mpi_bcast(latlast,1,MPI_INTEGER,0,mpi_comm_comp,iret)
216! call mpi_bcast(lonstart,1,MPI_INTEGER,0,mpi_comm_comp,iret)
217! call mpi_bcast(lonlast,1,MPI_INTEGER,0,mpi_comm_comp,iret)
218! call mpi_bcast(cenlat,1,MPI_INTEGER,0,mpi_comm_comp,iret)
219! call mpi_bcast(cenlon,1,MPI_INTEGER,0,mpi_comm_comp,iret)
220! print*, 'latstart,latlast A calling bcast=',latstart,latlast
221! print*,'lonstart,lonlast A calling bcast=',lonstart,lonlast
222! print*,'cenlat,cenlon A calling bcast=',cenlat,cenlon
223
224! if(me == 0)then
225! call nemsio_getheadvar(nfile,'global',global,iret)
226! if (iret /= 0) then
227! print*,"global not found in file-Assigned false"
228! global=.FALSE.
229! end if
230! end if
231! call mpi_bcast(global,1,MPI_LOGICAL,0,mpi_comm_comp,iret)
232
233! print*,'Is this a global run ',global
234 IF(.not. global)THEN
235! nframe=0 ! Wang added option to read without halos, so specify nframe=0
236 impf=im+nframe*2
237 jmpf=jm+nframe*2
238! nframed2=nframe/2
239 ELSE
240! nframe=1 !
241 impf=im+1 ! post cut im off because it's the same as i=1 but data from model is till im
242 jmpf=jm
243! nframed2=nframe/2
244 END IF
245
246 if (debugprint) then
247 print*,'impf,jmpf,nframe for reading fields = ',impf,jmpf,nframe
248 print*,'idate after broadcast = ',(idate(i),i=1,7)
249 print*,'nfhour = ',nfhour
250 end if
251 call mpi_scatterv(dummy(1,1),icnt,idsp,mpi_real &
252 ,dx(1,jsta),icnt(me),mpi_real,0,mpi_comm_comp,iret)
253 call mpi_scatterv(dummy2(1,1),icnt,idsp,mpi_real &
254 ,dy(1,jsta),icnt(me),mpi_real,0,mpi_comm_comp,iret)
255
256! print *,'before call EXCH,mype=',me,'max(gdlat)=',maxval(gdlat),'max(gdlon)=', &
257! maxval(gdlon)
258! CALL EXCH(gdlat(1,JSTA_2L))
259! print *,'after call EXCH,mype=',me
260
261 iyear = idate(1)
262 imn = idate(2) ! ask Jun
263 iday = idate(3) ! ask Jun
264 ihrst = idate(4)
265 imin = idate(5)
266 jdate = 0
267 idate = 0
268!
269! read(startdate,15)iyear,imn,iday,ihrst,imin
270 15 format(i4,1x,i2,1x,i2,1x,i2,1x,i2)
271 if (me==0) then
272 print*,'start yr mo day hr min =',iyear,imn,iday,ihrst,imin
273 print*,'processing yr mo day hr min=' &
274 ,idat(3),idat(1),idat(2),idat(4),idat(5)
275 endif
276!
277 idate(1) = iyear
278 idate(2) = imn
279 idate(3) = iday
280 idate(5) = ihrst
281 idate(6) = imin
282 sdat(1) = imn
283 sdat(2) = iday
284 sdat(3) = iyear
285 jdate(1) = idat(3)
286 jdate(2) = idat(1)
287 jdate(3) = idat(2)
288 jdate(5) = idat(4)
289 jdate(6) = idat(5)
290!
291
292! CALL W3DIFDAT(JDATE,IDATE,2,RINC)
293! ifhr=nint(rinc(2))
294!
295 CALL w3difdat(jdate,idate,0,rinc)
296!
297 ifhr=nint(rinc(2)+rinc(1)*24.)
298 ifmin=nint(rinc(3))
299! if(ifhr /= nfhour)print*,'find wrong Model input file';stop
300 if (me==0)print*,' in INITPOST ifhr ifmin fileName=',ifhr,ifmin,trim(filename)
301
302! Getting tstart
303 tstart=0.
304
305! Getiing restart
306
307 restrt=.true. ! set RESTRT as default
308! call ext_int_get_dom_ti_integer(DataHandle,'RESTARTBIN',itmp
309! + ,1,ioutcount,istatus)
310
311! IF(itmp < 1)THEN
312! RESTRT=.FALSE.
313! ELSE
314! RESTRT=.TRUE.
315! END IF
316
317! print*,'status for getting RESTARTBIN= ',istatus
318
319! print*,'Is this a restrt run? ',RESTRT
320
321 IF(tstart > 1.0e-2)THEN
322 ifhr=ifhr+nint(tstart)
323 rinc=0
324 idate=0
325 rinc(2)=-1.0*ifhr
326 call w3movdat(rinc,jdate,idate)
327 sdat(1)=idate(2)
328 sdat(2)=idate(3)
329 sdat(3)=idate(1)
330 ihrst=idate(5)
331 print*,'new forecast hours for restrt run= ',ifhr
332 print*,'new start yr mo day hr min =',sdat(3),sdat(1) &
333 ,sdat(2),ihrst,imin
334 END IF
335
336
337!KRF: Initialize extinction coef for aerosol to zero to avoid failure.
338! These are not in NEMS model output, but new CALVIS_GSD methods uses
339! these fields from ARW, and if not initialized here will cause failure.
340 extcof55=0.
341 aextc55=0.
342
343!Chuang: set default to Ferrier's MP scheme. NPS does not write MP option
344!used in model to nemsio output
345 varname='mp_physi'
346 if(me == 0)then
347 call nemsio_getheadvar(nfile,trim(varname),imp_physics,iret)
348 if (iret /= 0) then
349 print*,varname," not found in file- go to 16 character "
350 varname='mp_physics'
351 call nemsio_getheadvar(nfile,trim(varname),imp_physics,iret)
352 if (iret /= 0) then
353 print*,varname," not found in file-Assigned 1000"
354! imp_physics=1000 ! should this line be uncommented? - Moorthi
355 imp_physics=5
356 end if
357 end if
358 end if
359 call mpi_bcast(imp_physics,1,mpi_integer,0,mpi_comm_comp,iret)
360 if(me==0)print*,'MP_PHYSICS= ',imp_physics
361
362! Initializes constants for Ferrier microphysics
363 if(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95)then
364 CALL microinit(imp_physics)
365 end if
366
367 varname='sf_surface_physi'
368 if(me == 0)then
369 call nemsio_getheadvar(nfile,trim(varname),isf_surface_physics,iret)
370 if (iret /= 0) then
371 print*,varname," not found in file-Assigned 2 for NOAH LSM as default"
372 isf_surface_physics=2
373 end if
374 end if
375 call mpi_bcast(isf_surface_physics,1,mpi_integer,0,mpi_comm_comp,iret)
376 if(me==0) print*,'SF_SURFACE_PHYSICS= ',isf_surface_physics
377
378! IVEGSRC=1 for IGBP and 0 for USGS
379 varname='IVEGSRC'
380 if(me == 0)then
381 call nemsio_getheadvar(nfile,trim(varname),ivegsrc,iret)
382 if (iret /= 0) then
383 print*,varname," not found in file-Assigned 1 for IGBP as default"
384 ivegsrc=1
385 end if
386 end if
387 call mpi_bcast(ivegsrc,1,mpi_integer,0,mpi_comm_comp,iret)
388 if(me==0) print*,'IVEGSRC= ',ivegsrc
389
390! set novegtype based on vegetation classification
391 if(ivegsrc==1)then
392 novegtype=20
393 else if(ivegsrc==0)then
394 novegtype=24
395 end if
396 if(me==0) print*,'novegtype= ',novegtype
397
398 varname='CU_PHYSICS'
399 if(me == 0)then
400 call nemsio_getheadvar(nfile,trim(varname),icu_physics,iret)
401 if (iret /= 0) then
402 print*,varname," not found in file-Assigned 2 for BMJ as default"
403 icu_physics=2
404 end if
405 end if
406 call mpi_bcast(icu_physics,1,mpi_integer,0,mpi_comm_comp,iret)
407 if(me==0) print*,'CU_PHYSICS= ',icu_physics
408
409
410 allocate(bufy(jm))
411 varname='DX'
412! if(me == 0)then
413! call nemsio_getheadvar(nfile,trim(VarName),bufy,iret)
414! if (iret /= 0) then
415! print*,VarName," not found in file-Assigned missing values"
416! dx=spval
417! end if
418! end if
419! call mpi_bcast(bufy,jm,MPI_REAL,0,mpi_comm_comp,iret)
420! do j=jsta,jend
421! do i=1,im
422! dx(i,j)=bufy(j)
423! end do
424! end do
425 if(debugprint)print*,'sample ',varname,' = ',dx(im/2,(jsta+jend)/2)
426
427 varname='DY'
428! if(me == 0)then
429! call nemsio_getheadvar(nfile,trim(VarName),bufy,iret)
430! if (iret /= 0) then
431! print*,VarName," not found in file-Assigned missing values"
432! dx=spval
433! end if
434! end if
435! call mpi_bcast(bufy,jm,MPI_REAL,0,mpi_comm_comp,iret)
436! do j=jsta,jend
437! do i=1,im
438! dy(i,j)=bufy(j)
439! end do
440! end do
441 if(debugprint)print*,'sample ',varname,' = ',dy(im/2,(jsta+jend)/2)
442 deallocate(bufy)
443
444 varname='dt'
445 if(me == 0)then
446 call nemsio_getheadvar(nfile,trim(varname),garb,iret)
447 if (iret /= 0) then
448 print*,varname," not found in file-Assigned missing values"
449 dt=spval
450 else
451 dt=garb
452 end if
453 end if
454 call mpi_bcast(dt,1,mpi_real,0,mpi_comm_comp,iret)
455
456 varname='dphd'
457 if(me == 0)then
458 call nemsio_getheadvar(nfile,trim(varname),garb,iret)
459 if (iret /= 0) then
460 print*,varname," not found in file-Assigned missing values"
461 dyval=spval
462 else
463 dyval=garb*gdsdegr
464 end if
465 end if
466 call mpi_bcast(dyval,1,mpi_real,0,mpi_comm_comp,iret)
467! dyval=106 ! hard wire for AQ domain testing
468
469 varname='dlmd'
470 if(me == 0)then
471 call nemsio_getheadvar(nfile,trim(varname),garb,iret)
472 if (iret /= 0) then
473 print*,varname," not found in file-Assigned missing values"
474 dxval=spval
475 else
476 dxval=garb*gdsdegr
477 end if
478 end if
479 call mpi_bcast(dxval,1,mpi_real,0,mpi_comm_comp,iret)
480! dxval=124 ! hard wire for AQ domain testing
481
482 if(me==0) print*,'DX, DY, DT=',dxval,dyval,dt
483
484 varname='TPH0D'
485 if(me == 0)then
486 call nemsio_getheadvar(nfile,trim(varname),garb,iret)
487 if (iret /= 0) then
488 print*,varname," not found in file-Assigned missing values"
489 cenlat=spval
490 else
491 cenlat=nint(garb*gdsdegr)
492 end if
493 end if
494 call mpi_bcast(cenlat,1,mpi_integer,0,mpi_comm_comp,iret)
495
496 varname='TLM0D'
497 if(me == 0)then
498 call nemsio_getheadvar(nfile,trim(varname),garb,iret)
499 if (iret /= 0) then
500 print*,varname," not found in file-Assigned missing values"
501 cenlon=spval
502 else
503 if(grib=="grib2") then
504 cenlon=nint((garb+360.)*gdsdegr)
505 endif
506 end if
507 end if
508 call mpi_bcast(cenlon,1,mpi_integer,0,mpi_comm_comp,iret)
509
510! VarName='TRUELAT1'
511! call retrieve_index(index,VarName,varname_all,nrecs,iret)
512! if (iret /= 0) then
513! print*,VarName," not found in file"
514! else
515! call mpi_file_read_at(iunit,file_offset(index)+5*4 &
516! ,garb,1,mpi_real4, mpi_status_ignore, ierr)
517! if (ierr /= 0) then
518! print*,"Error reading ", VarName," using MPIIO"
519! else
520! print*,VarName, ' from MPIIO READ= ',garb
521! TRUELAT1=nint(garb*1000.)
522! write(6,*) 'truelat1= ', TRUELAT1
523! end if
524! end if
525
526! VarName='TRUELAT2'
527! call retrieve_index(index,VarName,varname_all,nrecs,iret)
528! if (iret /= 0) then
529! print*,VarName," not found in file"
530! else
531! call mpi_file_read_at(iunit,file_offset(index)+5*4 &
532! ,garb,1,mpi_real4, mpi_status_ignore, ierr)
533! if (ierr /= 0) then
534! print*,"Error reading ", VarName," using MPIIO"
535! else
536! print*,VarName, ' from MPIIO READ= ',garb
537! TRUELAT2=nint(garb*1000.)
538! write(6,*) 'truelat2= ', TRUELAT2
539! end if
540! end if
541
542! VarName='MAP_PROJ'
543! if(me == 0)then
544! call nemsio_getheadvar(nfile,trim(VarName),maptype,iret)
545! if (iret /= 0) then
546! print*,VarName," not found in file-Assigned 1000"
547! maptype=1000
548! end if
549! end if
550! call mpi_bcast(maptype,1,MPI_INTEGER,0,mpi_comm_comp,iret)
551 IF(.not. global)THEN
552 maptype=205 ! for Arakawa-B grid
553 gridtype='B'
554 ELSE
555 maptype=0 ! for global NMMB on latlon grid
556 gridtype='A' ! will put wind on mass point for now to make regular latlon
557 END IF
558 if(me==0) print*,'maptype and gridtype= ',maptype,gridtype
559
560 hbm2=1.0
561
562 varname='glat'
563 vcoordname='sfc'
564 l=1
565 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
566 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
567 ,l,impf,jmpf,nframe,gdlat)
568
569 call collect_loc(gdlat,dummy)
570! decides whether or not to convert to degree
571 if(me==0)then
572 if(maxval(abs(dummy))<pi)then ! convert from radian to degree
573 if(debugprint)print*,'convert from radian to degree'
574 dummy=dummy*180./pi
575 convert_rad_to_deg=.true.
576 end if
577 end if
578! print *,' gdlat1=',gdlat(1,:)
579! print *,' gdlatim=',gdlat(im,:)
580
581 call mpi_bcast(convert_rad_to_deg,1,mpi_logical,0,mpi_comm_comp,iret)
582 if(convert_rad_to_deg)call mpi_scatterv(dummy(1,1),icnt,idsp,mpi_real &
583 ,gdlat(1,jsta),icnt(me),mpi_real,0,mpi_comm_comp,iret)
584 if(debugprint)print*,'sample ',varname,' = ',gdlat(im/2,(jsta+jend)/2)
585 if(debugprint)print*,'max min lat=',maxval(gdlat),minval(gdlat),'im=',im, &
586 'jsta_2l=',jsta_2l,'jend_2u=',jend_2u
587 !call collect_loc(gdlat,dummy)
588 if(me==0.and.debugprint)print*,'after collect lat=',dummy(1,1),dummy(im,jm)
589 if(me==0)then
590 ii=(1+im)/2
591 jj=(1+jm)/2
592 latstart=nint(dummy(1,1)*gdsdegr)
593 latlast=nint(dummy(im,jm)*gdsdegr)
594 latnw=nint(dummy(1,jm)*gdsdegr)
595 latse=nint(dummy(im,1)*gdsdegr)
596! dyval=nint((dummy(1,2)-dummy(1,1))*1000.)
597! dyval=106 ! hard wire for AQ domain testing
598 if(mod(im,2)==0)then
599! cenlat=nint((dummy(ii,jj)+dummy(ii+1,jj)+dummy(ii+1,jj+1)+dummy(ii,jj+1))/4.0*1000.)
600 else
601! cenlat=nint(dummy(ii,jj)*1000.)
602 end if
603 print*,'latstart,latlast B bcast= ',latstart,latlast
604 end if
605 call mpi_bcast(latstart,1,mpi_integer,0,mpi_comm_comp,iret)
606 call mpi_bcast(latlast,1,mpi_integer,0,mpi_comm_comp,iret)
607! call mpi_bcast(dyval,1,MPI_INTEGER,0,mpi_comm_comp,iret)
608! call mpi_bcast(cenlat,1,MPI_INTEGER,0,mpi_comm_comp,iret)
609 if(debugprint) write(6,*) 'latstart,latlast,me A calling bcast=',latstart,latlast,me
610 if(me==0) print*,'dyval, cenlat= ',dyval, cenlat
611
612!$omp parallel do private(i,j)
613 do j=jsta,jend
614 do i=1,im
615 f(i,j) = 1.454441e-4*sin(gdlat(i,j)*dtr) ! 2*omeg*sin(phi)
616 end do
617 end do
618
619 varname='glon'
620 vcoordname='sfc'
621 l=1
622 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
623 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
624 ,l,impf,jmpf,nframe,gdlon)
625 if(convert_rad_to_deg)gdlon=gdlon*180./pi
626 if(global)then
627! do j=jsta,jend
628! do i=1,im
629! if(gdlon(i,j)<0.)gdlon(i,j)=360.+gdlon(i,j)
630! end do
631! end do
632 if(gdlon(1,jsta)>0. .and. gdlon(2,jsta)<0.)then
633 do j=jsta,jend
634 gdlon(1,j) = gdlon(1,j)-360.0
635 end do
636 end if
637 end if
638 if(debugprint)print*,'sample ',varname,' = ',(gdlon(i,(jsta+jend)/2),i=1,im,8)
639 if(debugprint)print*,'max min lon=',maxval(gdlon),minval(gdlon)
640 call collect_loc(gdlon,dummy)
641 if(me==0)then
642 if(grib=='grib2') then
643 if(dummy(1,1)<0) dummy(1,1)=dummy(1,1)+360.
644 if(dummy(im,jm)<0) dummy(im,jm)=dummy(im,jm)+360.
645 endif
646 lonstart=nint(dummy(1,1)*gdsdegr)
647 lonlast=nint(dummy(im,jm)*gdsdegr)
648 lonnw=nint(dummy(1,jm)*gdsdegr)
649 lonse=nint(dummy(im,1)*gdsdegr)
650! dxval=nint((dummy(2,1)-dummy(1,1))*1000.)
651! dxval=124 ! hard wire for AQ domain testing
652 if(mod(im,2)==0)then
653! cenlon=nint((dummy(ii,jj)+dummy(ii+1,jj)+dummy(ii+1,jj+1)+dummy(ii,jj+1))/4.0*1000.)
654 else
655! cenlon=nint(dummy(ii,jj)*1000.)
656 end if
657 end if
658 call mpi_bcast(lonstart,1,mpi_integer,0,mpi_comm_comp,iret)
659 call mpi_bcast(lonlast,1,mpi_integer,0,mpi_comm_comp,iret)
660! call mpi_bcast(dxval,1,MPI_INTEGER,0,mpi_comm_comp,iret)
661! call mpi_bcast(cenlon,1,MPI_INTEGER,0,mpi_comm_comp,iret)
662 write(6,*)'lonstart,lonlast A calling bcast=',lonstart,lonlast
663 print*,'dxval, cenlon= ',dxval, cenlon
664
665 convert_rad_to_deg=.false.
666 varname='vlat'
667 vcoordname='sfc'
668 l=1
669 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
670 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
671 ,l,impf,jmpf,nframe,buf)
672 if(debugprint)print*,'sample ',varname,' = ',buf(im/2,(jsta+jend)/2)
673 if(debugprint)print*,'max min vlat=',maxval(buf),minval(buf)
674 call collect_loc(buf,dummy)
675 if(me==0)then
676 if(maxval(abs(dummy))<pi)then ! convert from radian to degree
677 dummy(1,1)=dummy(1,1)*180./pi
678 dummy(im,jm)=dummy(im,jm)*180./pi
679 convert_rad_to_deg=.true.
680 end if
681 latstartv=nint(dummy(1,1)*gdsdegr)
682 latlastv=nint(dummy(im,jm)*gdsdegr)
683! cenlatv=nint(dummy(ii,jj)*1000.)
684! print*,'latstartv,cenlatv B bcast= ',latstartv,cenlatv
685 end if
686 call mpi_bcast(latstartv,1,mpi_integer,0,mpi_comm_comp,iret)
687 call mpi_bcast(latlastv,1,mpi_integer,0,mpi_comm_comp,iret)
688! call mpi_bcast(cenlatv,1,MPI_INTEGER,0,mpi_comm_comp,iret)
689 cenlatv=cenlat
690 if(debugprint)write(6,*) 'latstartv,cenlatv,latlastv,me A calling bcast=', &
691 latstartv,cenlatv,latlastv,me
692
693 varname='vlon'
694 vcoordname='sfc'
695 l=1
696 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
697 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
698 ,l,impf,jmpf,nframe,buf)
699 if(debugprint)print*,'sample ',varname,' = ',buf(im/2,(jsta+jend)/2)
700 if(debugprint)print*,'max min vlon=',maxval(buf),minval(buf)
701 call collect_loc(buf,dummy)
702 if(me==0)then
703 if(convert_rad_to_deg)then
704 dummy(1,1)=dummy(1,1)*180./pi
705 dummy(im,jm)=dummy(im,jm)*180./pi
706 end if
707 if(grib=='grib2') then
708 if(dummy(1,1)<0) dummy(1,1)=dummy(1,1)+360.
709 endif
710 lonstartv=nint(dummy(1,1)*gdsdegr)
711 lonlastv=nint(dummy(im,jm)*gdsdegr)
712! cenlonv=nint(dummy(ii,jj)*1000.)
713! print*,'lonstartv,cenlonv B bcast= ',lonstartv,cenlonv
714 end if
715 call mpi_bcast(lonstartv,1,mpi_integer,0,mpi_comm_comp,iret)
716 call mpi_bcast(lonlastv,1,mpi_integer,0,mpi_comm_comp,iret)
717! call mpi_bcast(cenlonv,1,MPI_INTEGER,0,mpi_comm_comp,iret)
718 cenlonv=cenlon
719 write(6,*) 'lonstartv,cenlonv,lonlastv,me A calling bcast=', &
720 lonstartv,cenlonv,lonlastv,me
721
722 varname='sm'
723 vcoordname='sfc'
724 l=1
725 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
726 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
727 ,l,impf,jmpf,nframe,sm)
728 if(debugprint)print*,'sample ',varname,' = ',sm(im/2,(jsta+jend)/2)
729
730 varname='sice'
731 vcoordname='sfc'
732 l=1
733 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
734 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
735 ,l,impf,jmpf,nframe,sice)
736 if(debugprint)print*,'sample ',varname,' = ',sice(im/2,(jsta+jend)/2)
737
738
739 varname='dpres'
740 vcoordname='hybrid sig lev'
741 l=1
742 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
743 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
744 ,l,impf,jmpf,nframe,pd)
745 if(debugprint)print*,'sample ',varname,' = ',pd(im/2,(jsta+jend)/2)
746
747! do j = jsta_2l, jend_2u
748! do i = 1, im
749! PD(I,J)=DUMMY2(I,J)
750! enddo
751! enddo
752
753 varname='hgt'
754 vcoordname='sfc'
755 l=1
756 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
757 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
758 ,l,impf,jmpf,nframe,fis)
759 if(debugprint)print*,'sample ',varname,' = ',fis(im/2,(jsta+jend)/2)
760 where(fis /= spval)fis=fis*g ! convert to geopotential
761
762 varname='tmp'
763 vcoordname='mid layer'
764 do l=1,lm
765! ll=lm-l+1
766 ll=l
767 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
768 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
769 ,l,impf,jmpf,nframe,t(1,jsta_2l,ll))
770 if(debugprint)print*,'sample l ',varname,' = ',ll,t(im/2,(jsta+jend)/2,ll)
771! if(debugprint)print*,'i=1 and im ',VarName,' = ',ll,t(1,(jsta+jend)/2,ll), &
772! t(im,(jsta+jend)/2,ll)
773 do i=1,im
774 do j=jsta,jend
775 if(t(i,j,ll)<150.)print*,'abnormal incoming T ',i,j,ll,t(i,j,ll)
776 end do
777 end do
778 end do ! do loop for l
779
780! Chuang: start mods to read variable for other mp schemes
781! model level q
782 if(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95 .or. imp_physics==99)then
783 varname='spfh'
784 vcoordname='mid layer'
785 do l=1,lm
786! ll=lm-l+1
787 ll=l
788 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
789 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
790 ,l,impf,jmpf,nframe,q(1,jsta_2l,ll))
791 if(debugprint)print*,'sample l ',varname,' = ',ll,q(im/2,(jsta+jend)/2,ll)
792 end do ! do loop for l
793 end if
794
795! model level u
796 varname='ugrd'
797 vcoordname='mid layer'
798 do l=1,lm
799! ll=lm-l+1
800 ll=l
801 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
802 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
803 ,l,impf,jmpf,nframe,uh(1,jsta_2l,ll))
804 if(debugprint)print*,'sample l ',varname,' = ',ll,uh(im/2,(jsta+jend)/2,ll)
805! put u on h point for global nmm
806 if(global)then
807 buf(:,:)=uh(:,:,ll)
808 call exch(buf(1,jsta_2l))
809 if(debugprint)print*,'sample l u = ',ll,buf(im/2,(jsta+jend)/2)
810 do j=jsta,jend
811 do i=1,im
812 im1=i-1
813 if(im1<1)im1=im1+im
814 jm1=j-1
815 if(j==1)then
816 ii=i+im/2
817 iim1=ii-1
818 if(iim1<1)iim1=iim1+im
819 if (ii > im) ii = ii - im
820 uh(i,j,ll)=(buf(i,j)+buf(im1,j)+buf(ii,j)+buf(iim1,j))/4.0
821 else
822 uh(i,j,ll)=(buf(i,j)+buf(im1,j)+buf(im1,jm1)+buf(i,jm1))/4.0
823 end if
824 end do
825 end do
826 end if ! end of wind interpolation for global NMM
827 end do ! do loop for l
828
829! model level v
830 varname='vgrd'
831 vcoordname='mid layer'
832 do l=1,lm
833! ll=lm-l+1
834 ll=l
835 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
836 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
837 ,l,impf,jmpf,nframe,vh(1,jsta_2l,ll))
838 if(debugprint)print*,'sample l ',varname,' = ',ll,vh(im/2,(jsta+jend)/2,ll)
839! put u on h point for global nmm
840 if(global)then
841 buf(:,:)=vh(:,:,ll)
842 call exch(buf(1,jsta_2l))
843 if(debugprint)print*,'sample l v = ',ll,buf(im/2,(jsta+jend)/2)
844 do j=jsta,jend
845 do i=1,im
846 im1=i-1
847 if(im1<1)im1=im1+im
848 jm1=j-1
849 if(j==1)then
850 ii=i+im/2
851 iim1=ii-1
852 if(iim1<1)iim1=iim1+im
853 if (ii > im) ii = ii - im
854 vh(i,j,ll)=(buf(i,j)+buf(im1,j)+buf(ii,j)+buf(iim1,j))/4.0
855 else
856 vh(i,j,ll)=(buf(i,j)+buf(im1,j)+buf(im1,jm1)+buf(i,jm1))/4.0
857 end if
858 end do
859 end do
860 end if ! end of wind interpolation for global NMM
861
862 end do ! do loop for l
863
864 varname='sg1'
865 if(me==0)then
866 call nemsio_getheadvar(nfile,trim(varname),eta1,iret)
867 if (iret /= 0) then
868 print*,varname," not found in file-Assigned missing values"
869 eta1=spval
870 end if
871 end if
872 call mpi_bcast(eta1,lm+1,mpi_real,0,mpi_comm_comp,iret)
873
874 varname='sg2'
875 if(me==0)then
876 call nemsio_getheadvar(nfile,trim(varname),eta2,iret)
877 if (iret /= 0) then
878 print*,varname," not found in file-Assigned missing values"
879 eta2=spval
880 end if
881 end if
882 call mpi_bcast(eta2,lm+1,mpi_real,0,mpi_comm_comp,iret)
883
884 open(75,file='ETAPROFILE.txt',form='formatted', &
885 status='unknown')
886 DO l=1,lm+1
887 write(75,1020)l, eta1(lm+2-l), eta2(lm+2-l)
888 END DO
889 1020 format(i3,2e17.10)
890 close (75)
891
892 varname='pdtop'
893 if(me==0)then
894 call nemsio_getheadvar(nfile,trim(varname),pdtop,iret)
895 if (iret /= 0) then
896 print*,varname," not found in file-Assigned missing values"
897 pdtop=spval
898 end if
899 end if
900 call mpi_bcast(pdtop,1,mpi_real,0,mpi_comm_comp,iret)
901
902 varname='pt'
903 if(me==0)then
904 call nemsio_getheadvar(nfile,trim(varname),pt,iret)
905 if (iret /= 0) then
906 print*,varname," not found in file-Assigned missing values"
907 pt=spval
908 end if
909 end if
910 call mpi_bcast(pt,1,mpi_real,0,mpi_comm_comp,iret)
911 if(me==0) print*,'PT, PDTOP= ',pt,pdtop
912
913 varname='pblh'
914 vcoordname='sfc'
915 l=1
916 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
917 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
918 ,l,impf,jmpf,nframe,pblh)
919 if(debugprint)print*,'sample ',varname,' = ',pblh(im/2,(jsta+jend)/2)
920
921 varname='mixht'
922 vcoordname='sfc'
923 l=1
924 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
925 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
926 ,l,impf,jmpf,nframe,mixht)
927 if(debugprint)print*,'sample ',varname,' = ',mixht(im/2,(jsta+jend)/2)
928
929 varname='uustar'
930 vcoordname='sfc'
931 l=1
932 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
933 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
934 ,l,impf,jmpf,nframe,ustar)
935 if(debugprint)print*,'sample ',varname,' = ',ustar(im/2,(jsta+jend)/2)
936
937 varname='zorl'
938 vcoordname='sfc'
939 l=1
940 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
941 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
942 ,l,impf,jmpf,nframe,z0)
943 if(debugprint)print*,'sample ',varname,' = ',z0(im/2,(jsta+jend)/2)
944
945 varname='ths'
946 vcoordname='sfc'
947 l=1
948 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
949 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
950 ,l,impf,jmpf,nframe,ths)
951 if(debugprint)print*,'sample ',varname,' = ',ths(im/2,(jsta+jend)/2)
952
953 varname='qsh'
954 vcoordname='sfc'
955 l=1
956 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
957 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
958 ,l,impf,jmpf,nframe,qs)
959 if(debugprint)print*,'sample ',varname,' = ',qs(im/2,(jsta+jend)/2)
960
961 varname='twbs'
962 vcoordname='sfc'
963 l=1
964 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
965 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
966 ,l,impf,jmpf,nframe,twbs)
967 if(debugprint)print*,'sample ',varname,' = ',twbs(im/2,(jsta+jend)/2)
968
969 varname='qwbs'
970 vcoordname='sfc'
971 l=1
972 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
973 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
974 ,l,impf,jmpf,nframe,qwbs)
975 if(debugprint)print*,'sample ',varname,' = ',qwbs(im/2,(jsta+jend)/2)
976
977 varname='prec' ! instantaneous precip rate?
978 vcoordname='sfc'
979 l=1
980 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
981 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
982 ,l,impf,jmpf,nframe,prec)
983 if(debugprint)print*,'sample ',varname,' = ',prec(im/2,(jsta+jend)/2)
984
985 varname='acprec' ! accum total precip
986 vcoordname='sfc'
987 l=1
988 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
989 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
990 ,l,impf,jmpf,nframe,acprec)
991 if(debugprint)print*,'sample ',varname,' = ',acprec(im/2,(jsta+jend)/2)
992
993 varname='cuprec' ! accum cumulus precip
994 vcoordname='sfc'
995 l=1
996 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
997 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
998 ,l,impf,jmpf,nframe,cuprec)
999 if(debugprint)print*,'sample ',varname,' = ',cuprec(im/2,(jsta+jend)/2)
1000
1001! compute grid scale precip
1002 do j=jsta,jend
1003 do i=1,im
1004 ancprc(i,j)=acprec(i,j)-cuprec(i,j)
1005 end do
1006 end do
1007
1008 varname='lspa'
1009 vcoordname='sfc'
1010 l=1
1011 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1012 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1013 ,l,impf,jmpf,nframe,lspa)
1014 if(debugprint)print*,'sample ',varname,' = ',lspa(im/2,(jsta+jend)/2)
1015
1016 varname='sno'
1017 vcoordname='sfc'
1018 l=1
1019 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1020 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1021 ,l,impf,jmpf,nframe,sno)
1022 if(debugprint)print*,'sample ',varname,' = ',sno(im/2,(jsta+jend)/2)
1023
1024 varname='snoavg'
1025 vcoordname='sfc'
1026 l=1
1027 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1028 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1029 ,l,impf,jmpf,nframe,snoavg)
1030 if(debugprint)print*,'sample ',varname,' = ',snoavg(im/2,(jsta+jend)/2)
1031
1032 varname='psfcavg'
1033 vcoordname='sfc'
1034 l=1
1035 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1036 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1037 ,l,impf,jmpf,nframe,psfcavg)
1038 if(debugprint)print*,'sample ',varname,' = ',psfcavg(im/2,(jsta+jend)/2)
1039
1040 varname='t10avg'
1041 vcoordname='10 m above gnd'
1042 l=1
1043 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1044 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1045 ,l,impf,jmpf,nframe,t10avg)
1046 if(debugprint)print*,'sample ',varname,' = ',t10avg(im/2,(jsta+jend)/2)
1047
1048 varname='t10'
1049 vcoordname='10 m above gnd'
1050 l=1
1051 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1052 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1053 ,l,impf,jmpf,nframe,t10m)
1054 if(debugprint)print*,'sample ',varname,' = ',t10m(im/2,(jsta+jend)/2)
1055
1056 varname='akhsavg'
1057 vcoordname='sfc'
1058 l=1
1059 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1060 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1061 ,l,impf,jmpf,nframe,akhsavg)
1062 if(debugprint)print*,'sample ',varname,' = ',akhsavg(im/2,(jsta+jend)/2)
1063
1064 varname='akmsavg'
1065 vcoordname='sfc'
1066 l=1
1067 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1068 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1069 ,l,impf,jmpf,nframe,akmsavg)
1070 if(debugprint)print*,'sample ',varname,' = ',akmsavg(im/2,(jsta+jend)/2)
1071
1072 varname='refdmax'
1073 vcoordname='sfc'
1074 l=1
1075 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1076 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1077 ,l,impf,jmpf,nframe,refd_max)
1078 if(debugprint)print*,'sample ',varname,' = ',refd_max(im/2,(jsta+jend)/2)
1079
1080 varname='upvvelmax'
1081 vcoordname='sfc' ! wrong
1082 l=1
1083 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1084 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1085 ,l,impf,jmpf,nframe,w_up_max)
1086 if(debugprint)print*,'sample ',varname,' = ',w_up_max(im/2,(jsta+jend)/2)
1087
1088 varname='dnvvelmax'
1089 vcoordname='sfc' ! wrong
1090 l=1
1091 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1092 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1093 ,l,impf,jmpf,nframe,w_dn_max)
1094 if(debugprint)print*,'sample ',varname,' = ',w_dn_max(im/2,(jsta+jend)/2)
1095
1096 varname='uphlmax'
1097 vcoordname='sfc' ! wrong
1098 l=1
1099 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1100 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1101 ,l,impf,jmpf,nframe,up_heli_max)
1102 if(debugprint)print*,'sample ',varname,' = ',up_heli_max(im/2,(jsta+jend)/2)
1103
1104 varname='pratemax' ! max hourly precip rate
1105 vcoordname='sfc' ! wrong
1106 l=1
1107 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1108 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1109 ,l,impf,jmpf,nframe,prate_max)
1110 if(debugprint)print*,'sample ',varname,' = ',prate_max(im/2,(jsta+jend)/2)
1111
1112 varname='fpratemax' ! max hourly frozen precip rate
1113 vcoordname='sfc' ! wrong
1114 l=1
1115 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1116 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1117 ,l,impf,jmpf,nframe,fprate_max)
1118 if(debugprint)print*,'sample ',varname,' = ',fprate_max(im/2,(jsta+jend)/2)
1119
1120 varname='si'
1121 vcoordname='sfc'
1122 l=1
1123 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1124 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1125 ,l,impf,jmpf,nframe,si)
1126 if(debugprint)print*,'sample ',varname,' = ',si(im/2,(jsta+jend)/2)
1127
1128 varname='cldefi'
1129 vcoordname='sfc'
1130 l=1
1131 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1132 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1133 ,l,impf,jmpf,nframe,cldefi)
1134 if(debugprint)print*,'sample ',varname,' = ',cldefi(im/2,(jsta+jend)/2)
1135
1136 varname='th10'
1137 vcoordname='10 m above gnd'
1138 l=1
1139 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1140 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1141 ,l,impf,jmpf,nframe,th10)
1142 if(debugprint)print*,'sample ',varname,' = ',th10(im/2,(jsta+jend)/2)
1143
1144 varname='q10'
1145 vcoordname='10 m above gnd'
1146 l=1
1147 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1148 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1149 ,l,impf,jmpf,nframe,q10)
1150 if(debugprint)print*,'sample ',varname,' = ',q10(im/2,(jsta+jend)/2)
1151
1152 varname='pshltr'
1153 vcoordname='sfc'
1154 l=1
1155 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1156 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1157 ,l,impf,jmpf,nframe,pshltr)
1158 if(debugprint)print*,'sample ',varname,' = ',pshltr(im/2,(jsta+jend)/2), &
1159 'max=',maxval(pshltr(1:im,jsta:jend)),minval(pshltr(1:im,jsta:jend))
1160
1161 varname='tshltr'
1162 vcoordname='sfc'
1163 l=1
1164 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1165 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1166 ,l,impf,jmpf,nframe,tshltr)
1167 if(debugprint)print*,'sample ',varname,' = ',tshltr(im/2,(jsta+jend)/2)
1168
1169 varname='qshltr'
1170 vcoordname='sfc'
1171 l=1
1172 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1173 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1174 ,l,impf,jmpf,nframe,qshltr)
1175 if(debugprint)print*,'sample ',varname,' = ',qshltr(im/2,(jsta+jend)/2)
1176
1177 tmaxmin=1.
1178! call mpi_bcast(TMAXMIN,1,MPI_REAL,0,mpi_comm_comp,iret)
1179
1180 varname='t02max'
1181 vcoordname='sfc'
1182 l=1
1183 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1184 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1185 ,l,impf,jmpf,nframe,maxtshltr)
1186 if(debugprint)print*,'sample ',varname,' = ',maxtshltr(im/2,(jsta+jend)/2)
1187
1188 varname='t02min'
1189 vcoordname='sfc'
1190 l=1
1191 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1192 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1193 ,l,impf,jmpf,nframe,mintshltr)
1194 if(debugprint)print*,'sample ',varname,' = ',mintshltr(im/2,(jsta+jend)/2)
1195
1196 varname='rh02max'
1197 vcoordname='sfc'
1198 l=1
1199 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1200 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1201 ,l,impf,jmpf,nframe,maxrhshltr)
1202 if(debugprint)print*,'sample ',varname,' = ',maxrhshltr(im/2,(jsta+jend)/2)
1203
1204 varname='rh02min'
1205 vcoordname='sfc'
1206 l=1
1207 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1208 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1209 ,l,impf,jmpf,nframe,minrhshltr)
1210 if(debugprint)print*,'sample ',varname,' = ',minrhshltr(im/2,(jsta+jend)/2)
1211
1212! model level q2
1213 varname='q2'
1214 vcoordname='mid layer'
1215 do l=1,lm
1216! ll=lm-l+1
1217 ll=l
1218 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1219 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1220 ,l,impf,jmpf,nframe,q2(1,jsta_2l,ll))
1221 if(debugprint)print*,'sample l ',varname,' = ',ll,q2(im/2,(jsta+jend)/2,ll)
1222 end do ! do loop for l
1223
1224 varname='akhs_out'
1225 vcoordname='sfc'
1226 l=1
1227 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1228 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1229 ,l,impf,jmpf,nframe,akhs)
1230 if(debugprint)print*,'sample ',varname,' = ',akhs(im/2,(jsta+jend)/2)
1231
1232 varname='akms_out'
1233 vcoordname='sfc'
1234 l=1
1235 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1236 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1237 ,l,impf,jmpf,nframe,akms)
1238 if(debugprint)print*,'sample ',varname,' = ',akms(im/2,(jsta+jend)/2)
1239
1240 varname='albase'
1241 vcoordname='sfc'
1242 l=1
1243 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1244 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1245 ,l,impf,jmpf,nframe,albase)
1246 if(debugprint)print*,'sample ',varname,' = ',albase(im/2,(jsta+jend)/2)
1247
1248 varname='albedo'
1249 vcoordname='sfc'
1250 l=1
1251 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1252 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1253 ,l,impf,jmpf,nframe,albedo)
1254 if(debugprint)print*,'sample ',varname,' = ',albedo(im/2,(jsta+jend)/2)
1255
1256 varname='czen'
1257 vcoordname='sfc'
1258 l=1
1259 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1260 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1261 ,l,impf,jmpf,nframe,czen)
1262 if(debugprint)print*,'sample ',varname,' = ',czen(im/2,(jsta+jend)/2)
1263
1264 varname='czmean'
1265 vcoordname='sfc'
1266 l=1
1267 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1268 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1269 ,l,impf,jmpf,nframe,czmean)
1270 if(debugprint)print*,'sample ',varname,' = ',czmean(im/2,(jsta+jend)/2)
1271
1272 varname='mxsnal'
1273 vcoordname='sfc'
1274 l=1
1275 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1276 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1277 ,l,impf,jmpf,nframe,mxsnal)
1278 if(debugprint)print*,'sample ',varname,' = ',mxsnal(im/2,(jsta+jend)/2)
1279
1280 varname='radot'
1281 vcoordname='sfc'
1282 l=1
1283 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1284 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1285 ,l,impf,jmpf,nframe,radot)
1286 if(debugprint)print*,'sample ',varname,' = ',radot(im/2,(jsta+jend)/2)
1287
1288 varname='sigt4'
1289 vcoordname='sfc'
1290 l=1
1291 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1292 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1293 ,l,impf,jmpf,nframe,sigt4)
1294 if(debugprint)print*,'sample ',varname,' = ',sigt4(im/2,(jsta+jend)/2)
1295
1296 varname='tg'
1297 vcoordname='sfc'
1298 l=1
1299 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1300 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1301 ,l,impf,jmpf,nframe,tg)
1302 if(debugprint)print*,'sample ',varname,' = ',tg(im/2,(jsta+jend)/2)
1303
1304! Chuang: mods to read variables from other mp schemes
1305 if(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95 .or. imp_physics==99)then
1306! model level cwm
1307 varname='clwmr'
1308 vcoordname='mid layer'
1309 do l=1,lm
1310! ll=lm-l+1
1311 ll=l
1312 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1313 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1314 ,l,impf,jmpf,nframe,cwm(1,jsta_2l,ll))
1315 if(debugprint)print*,'sample l ',varname,' = ',ll,cwm(im/2,(jsta+jend)/2,ll)
1316 end do ! do loop for l
1317
1318 varname='f_ice'
1319 vcoordname='mid layer'
1320 do l=1,lm
1321! ll=lm-l+1
1322 ll=l
1323 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1324 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1325 ,l,impf,jmpf,nframe,f_ice(1,jsta_2l,ll))
1326 if(debugprint)print*,'sample l ',varname,' = ',ll,f_ice(im/2,(jsta+jend)/2,ll)
1327 if(debugprint)print*,'max min ',varname,' = ',ll,maxval(f_ice(:,:,ll)),minval(f_ice(:,:,ll))
1328 end do ! do loop for l
1329
1330 varname='f_rain'
1331 vcoordname='mid layer'
1332 do l=1,lm
1333! ll=lm-l+1
1334 ll=l
1335 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1336 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1337 ,l,impf,jmpf,nframe,f_rain(1,jsta_2l,ll))
1338 if(debugprint)print*,'sample l ',varname,' = ',ll,f_rain(im/2,(jsta+jend)/2,ll)
1339 if(debugprint)print*,'max min ',varname,' = ',ll,maxval(f_rain(:,:,ll)),minval(f_rain(:,:,ll))
1340 end do ! do loop for l
1341
1342 if(imp_physics /= 99)then
1343 varname='f_rimef'
1344 vcoordname='mid layer'
1345 do l=1,lm
1346! ll=lm-l+1
1347 ll=l
1348 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1349 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1350 ,l,impf,jmpf,nframe,f_rimef(1,jsta_2l,ll))
1351 if(debugprint)print*,'sample l ',varname,' = ',ll,f_rimef(im/2,(jsta+jend)/2,ll)
1352 if(debugprint)print*,'max min ',varname,' = ',ll,maxval(f_rimef(:,:,ll)),minval(f_rimef(:,:,ll))
1353 end do ! do loop for l
1354 endif
1355
1356 if(imp_physics==5)then
1357 varname='refl_10cm'
1358 vcoordname='mid layer'
1359 do l=1,lm
1360 ll=l
1361 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1362 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1363 ,l,impf,jmpf,nframe,ref_10cm(1,jsta_2l,ll))
1364 if(debugprint)print*,'sample l ',varname,' = ' &
1365 ,ll,ref_10cm(im/2,(jsta+jend)/2,ll)
1366 if(debugprint)print*,'max min ',varname,' = ' &
1367 ,ll,maxval(ref_10cm(:,:,ll)),minval(ref_10cm(:,:,ll))
1368 end do ! do loop for l
1369 endif
1370
1371 else ! retrieve hydrometeo fields directly for non-Ferrier
1372 if(imp_physics==8 .or. imp_physics==9)then
1373 varname='ni'
1374 vcoordname='mid layer'
1375 do l=1,lm
1376 ll=l
1377 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1378 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1379 ,l,impf,jmpf,nframe,qqni(1,jsta_2l,ll))
1380 if(debugprint)print*,'sample l ',varname,' = ' &
1381 ,ll,qqni(im/2,(jsta+jend)/2,ll)
1382 if(debugprint)print*,'max min ',varname,' = ' &
1383 ,ll,maxval(qqni(:,:,ll)),minval(qqni(:,:,ll))
1384 end do ! do loop for l
1385
1386 varname='nr'
1387 vcoordname='mid layer'
1388 do l=1,lm
1389 ll=l
1390 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1391 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1392 ,l,impf,jmpf,nframe,qqnr(1,jsta_2l,ll))
1393 if(debugprint)print*,'sample l ',varname,' = ' &
1394 ,ll,qqnr(im/2,(jsta+jend)/2,ll)
1395 if(debugprint)print*,'max min ',varname,' = ' &
1396 ,ll,maxval(qqnr(:,:,ll)),minval(qqnr(:,:,ll))
1397 end do ! do loop for l
1398 end if !imp_physics==8 or 9
1399
1400! water vapor, will be converted to specific humidity
1401 varname='qv'
1402 vcoordname='mid layer'
1403 do l=1,lm
1404 ll=l
1405 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1406 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1407 ,l,impf,jmpf,nframe,q(1,jsta_2l,ll))
1408!$omp parallel do private(i,j)
1409 do j=jsta_2l,jend_2u
1410 do i=1,im
1411 q(i,j,ll) = max(10e-12,q(i,j,ll))
1412 q(i,j,ll) = q(i,j,ll) / (1.0+q(i,j,ll))
1413 end do
1414 end do
1415 if(debugprint)print*,'sample l ',varname,' = ' &
1416 ,ll,q(im/2,(jsta+jend)/2,ll)
1417 if(debugprint)print*,'max min ',varname,' = ' &
1418 ,ll,maxval(q(:,:,ll)),minval(q(:,:,ll))
1419 end do ! ,j,l) do loop for l
1420
1421!$omp parallel do private(i,j,l)
1422 do l=1,lm
1423 do j=jsta_2l,jend_2u
1424 do i=1,im
1425 qqw(i,j,l) = spval
1426 qqr(i,j,l) = spval
1427 qqs(i,j,l) = spval
1428 qqi(i,j,l) = spval
1429 qqg(i,j,l) = spval
1430 cwm(i,j,l) = spval
1431 enddo
1432 enddo
1433 enddo
1434! cloud water
1435 if(imp_physics/=0)then
1436 varname='qc'
1437 vcoordname='mid layer'
1438 do l=1,lm
1439 ll=l
1440 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1441 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1442 ,l,impf,jmpf,nframe,qqw(1,jsta_2l,ll))
1443 if(imp_physics==3)then !WSM3
1444!$omp parallel do private(i,j)
1445 do j=jsta_2l,jend_2u
1446 do i=1,im
1447 if(t(i,j,ll)<tfrz) qqi(i,j,ll) = qqw(i,j,ll)
1448 end do
1449 end do
1450 end if
1451!$omp parallel do private(i,j)
1452 do j=jsta_2l,jend_2u
1453 do i=1,im
1454 cwm(i,j,ll) = qqw(i,j,ll)
1455 end do
1456 end do
1457 if(debugprint)print*,'sample l ',varname,' = ' &
1458 ,ll,qqw(im/2,(jsta+jend)/2,ll)
1459 if(debugprint)print*,'max min ',varname,' = ' &
1460 ,ll,maxval(qqw(:,:,ll)),minval(qqw(:,:,ll))
1461 end do ! do loop for l
1462 end if ! for non-zero mp_physics
1463
1464! cloud ice
1465 if(imp_physics/=0 .and. imp_physics/=1 .and. imp_physics/=3)then
1466 varname='qi'
1467 vcoordname='mid layer'
1468 do l=1,lm
1469 ll=l
1470 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1471 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1472 ,l,impf,jmpf,nframe,qqi(1,jsta_2l,ll))
1473!$omp parallel do private(i,j)
1474 do j=jsta_2l,jend_2u
1475 do i=1,im
1476 cwm(i,j,ll) = cwm(i,j,ll) + qqi(i,j,ll)
1477 end do
1478 end do
1479 if(debugprint)print*,'sample l ',varname,' = ' &
1480 ,ll,qqi(im/2,(jsta+jend)/2,ll)
1481 if(debugprint)print*,'max min ',varname,' = ' &
1482 ,ll,maxval(qqi(:,:,ll)),minval(qqi(:,:,ll))
1483 end do ! do loop for l
1484 end if
1485! rain
1486 if(imp_physics/=0)then
1487 varname='qr'
1488 vcoordname='mid layer'
1489 do l=1,lm
1490 ll=l
1491 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1492 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1493 ,l,impf,jmpf,nframe,qqr(1,jsta_2l,ll))
1494 if(imp_physics==3)then !WSM3
1495!$omp parallel do private(i,j)
1496 do j=jsta_2l,jend_2u
1497 do i=1,im
1498 if(t(i,j,ll)<tfrz) qqs(i,j,ll) = qqr(i,j,ll)
1499 end do
1500 end do
1501 end if
1502!$omp parallel do private(i,j)
1503 do j=jsta_2l,jend_2u
1504 do i=1,im
1505 cwm(i,j,ll) = cwm(i,j,ll) + qqr(i,j,ll)
1506 end do
1507 end do
1508 if(debugprint)print*,'sample l ',varname,' = ' &
1509 ,ll,qqr(im/2,(jsta+jend)/2,ll)
1510 if(debugprint)print*,'max min ',varname,' = ' &
1511 ,ll,maxval(qqr(:,:,ll)),minval(qqr(:,:,ll))
1512 end do ! do loop for l
1513 end if ! for non-zero mp_physics
1514!snow
1515 if(imp_physics/=0 .and. imp_physics/=1 .and. imp_physics/=3)then
1516 varname='qs'
1517 vcoordname='mid layer'
1518 do l=1,lm
1519 ll=l
1520 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1521 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1522 ,l,impf,jmpf,nframe,qqs(1,jsta_2l,ll))
1523!$omp parallel do private(i,j)
1524 do j=jsta_2l,jend_2u
1525 do i=1,im
1526 cwm(i,j,ll) = cwm(i,j,ll) + qqs(i,j,ll)
1527 end do
1528 end do
1529 if(debugprint)print*,'sample l ',varname,' = ' &
1530 ,ll,qqs(im/2,(jsta+jend)/2,ll)
1531 if(debugprint)print*,'max min ',varname,' = ' &
1532 ,ll,maxval(qqs(:,:,ll)),minval(qqs(:,:,ll))
1533 end do ! do loop for l
1534 end if
1535!graupel
1536 if(imp_physics==2 .or. imp_physics==6 .or. imp_physics==8)then
1537 varname='qg'
1538 vcoordname='mid layer'
1539 do l=1,lm
1540 ll=l
1541 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1542 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1543 ,l,impf,jmpf,nframe,qqg(1,jsta_2l,ll))
1544!$omp parallel do private(i,j)
1545 do j=jsta_2l,jend_2u
1546 do i=1,im
1547 cwm(i,j,ll) = cwm(i,j,ll) + qqg(i,j,ll)
1548 end do
1549 end do
1550 if(debugprint)print*,'sample l ',varname,' = ' &
1551 ,ll,qqg(im/2,(jsta+jend)/2,ll)
1552 if(debugprint)print*,'max min ',varname,' = ' &
1553 ,ll,maxval(qqg(:,:,ll)),minval(qqg(:,:,ll))
1554 end do ! do loop for l
1555 end if
1556! radar ref and effective radius for Thompson
1557 if(imp_physics==8)then
1558 varname='refl_10cm'
1559 vcoordname='mid layer'
1560 do l=1,lm
1561 ll=l
1562 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1563 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1564 ,l,impf,jmpf,nframe,ref_10cm(1,jsta_2l,ll))
1565 if(debugprint)print*,'sample l ',varname,' = ' &
1566 ,ll,ref_10cm(im/2,(jsta+jend)/2,ll)
1567 if(debugprint)print*,'max min ',varname,' = ' &
1568 ,ll,maxval(ref_10cm(:,:,ll)),minval(ref_10cm(:,:,ll))
1569 end do ! do loop for l
1570! effective radius
1571 varname='re_cloud'
1572 vcoordname='mid layer'
1573 do l=1,lm
1574 ll=l
1575 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1576 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1577 ,l,impf,jmpf,nframe,radius_cloud(1,jsta_2l,ll))
1578 if(debugprint)print*,'sample l ',varname,' = ' &
1579 ,ll,radius_cloud(im/2,(jsta+jend)/2,ll)
1580 if(debugprint)print*,'max min ',varname,' = ' &
1581 ,ll,maxval(radius_cloud(:,:,ll)),minval(radius_cloud(:,:,ll))
1582 end do ! do loop for l
1583
1584 varname='re_ice'
1585 vcoordname='mid layer'
1586 do l=1,lm
1587 ll=l
1588 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1589 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1590 ,l,impf,jmpf,nframe,radius_ice(1,jsta_2l,ll))
1591 if(debugprint)print*,'sample l ',varname,' = ' &
1592 ,ll,radius_ice(im/2,(jsta+jend)/2,ll)
1593 if(debugprint)print*,'max min ',varname,' = ' &
1594 ,ll,maxval(radius_ice(:,:,ll)),minval(radius_ice(:,:,ll))
1595 end do ! do loop for l
1596
1597 varname='re_snow'
1598 vcoordname='mid layer'
1599 do l=1,lm
1600 ll=l
1601 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1602 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1603 ,l,impf,jmpf,nframe,radius_snow(1,jsta_2l,ll))
1604 if(debugprint)print*,'sample l ',varname,' = ' &
1605 ,ll,radius_snow(im/2,(jsta+jend)/2,ll)
1606 if(debugprint)print*,'max min ',varname,' = ' &
1607 ,ll,maxval(radius_snow(:,:,ll)),minval(radius_snow(:,:,ll))
1608 end do ! do loop for l
1609
1610 end if
1611 ! end of retrieving hydrometeo for non Ferrier's MP options
1612 end if ! end of retrieving hydrometeo for all MP options
1613
1614
1615 varname='cldfra'
1616 vcoordname='mid layer'
1617 do l=1,lm
1618! ll=lm-l+1
1619 ll=l
1620 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1621 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1622 ,l,impf,jmpf,nframe,cfr(1,jsta_2l,ll))
1623 if(debugprint)print*,'sample l ',varname,' = ',ll,cfr(im/2,(jsta+jend)/2,ll)
1624 end do ! do loop for l
1625
1626 varname='sr'
1627 vcoordname='sfc'
1628 l=1
1629 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1630 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1631 ,l,impf,jmpf,nframe,sr)
1632 if(debugprint)print*,'sample ',varname,' = ',sr(im/2,(jsta+jend)/2)
1633
1634 varname='cfrach'
1635 vcoordname='sfc'
1636 l=1
1637 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1638 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1639 ,l,impf,jmpf,nframe,cfrach)
1640 if(debugprint)print*,'sample ',varname,' = ',cfrach(im/2,(jsta+jend)/2)
1641
1642 varname='cfracl'
1643 vcoordname='sfc'
1644 l=1
1645 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1646 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1647 ,l,impf,jmpf,nframe,cfracl)
1648 if(debugprint)print*,'sample ',varname,' = ',cfracl(im/2,(jsta+jend)/2)
1649
1650 varname='cfracm'
1651 vcoordname='sfc'
1652 l=1
1653 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1654 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1655 ,l,impf,jmpf,nframe,cfracm)
1656 if(debugprint)print*,'sample ',varname,' = ',cfracm(im/2,(jsta+jend)/2)
1657
1658 varname='islope' !???
1659 vcoordname='sfc'
1660 l=1
1661 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1662 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1663 ,l,impf,jmpf,nframe,islope)
1664 if(debugprint)print*,'sample ',varname,' = ',islope(im/2,(jsta+jend)/2)
1665
1666! either assign SLDPTH to be the same as eta (which is original
1667! setup in WRF LSM) or extract thickness of soil layers from wrf
1668! output
1669
1670! or get SLDPTH from wrf output
1671 varname='sldpth'
1672 if(me==0)then
1673 call nemsio_getheadvar(nfile,trim(varname),sldpth,iret)
1674! if (iret /= 0) then
1675! print*,VarName," not found in file-Assigned NOAH default values"
1676! SLDPTH(1)=0.10
1677! SLDPTH(2)=0.3
1678! SLDPTH(3)=0.6
1679! SLDPTH(4)=1.0
1680! end if
1681 end if
1682 call mpi_bcast(sldpth,4,mpi_real,0,mpi_comm_comp,iret)
1683 if(me==0)print*,'SLDPTH= ',(sldpth(n),n=1,nsoil)
1684
1685 varname='cmc'
1686 vcoordname='sfc'
1687 l=1
1688 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1689 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1690 ,l,impf,jmpf,nframe,cmc)
1691 if(debugprint)print*,'sample ',varname,' = ',cmc(im/2,(jsta+jend)/2)
1692
1693 varname='grnflx'
1694 vcoordname='sfc'
1695 l=1
1696 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1697 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1698 ,l,impf,jmpf,nframe,grnflx)
1699 if(debugprint)print*,'sample ',varname,' = ',grnflx(im/2,(jsta+jend)/2)
1700
1701 varname='pctsno'
1702 vcoordname='sfc'
1703 l=1
1704 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1705 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1706 ,l,impf,jmpf,nframe,pctsno)
1707 if(debugprint)print*,'sample ',varname,' = ',pctsno(im/2,(jsta+jend)/2)
1708
1709 varname='soiltb'
1710 vcoordname='sfc'
1711 l=1
1712 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1713 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1714 ,l,impf,jmpf,nframe,soiltb)
1715 if(debugprint)print*,'sample ',varname,' = ',soiltb(im/2,(jsta+jend)/2)
1716 if(debugprint)then
1717 do j=jsta,jend
1718 do i=1,im
1719 if(soiltb(i,j)>350.)print*,'large soiltb='
1720 end do
1721 end do
1722 end if
1723
1724 varname='vegfrc'
1725 vcoordname='sfc'
1726 l=1
1727 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1728 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1729 ,l,impf,jmpf,nframe,vegfrc)
1730 if(debugprint)print*,'sample ',varname,' = ',vegfrc(im/2,(jsta+jend)/2)
1731
1732 varname='sh2o'
1733 vcoordname='soil layer'
1734 l=1
1735 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1736 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1737 ,l,impf,jmpf,nframe,sh2o(1,jsta_2l,l))
1738 if(debugprint)print*,'sample l ',varname,' = ',l,sh2o(im/2,(jsta+jend)/2,l)
1739
1740 varname='sh2o'
1741 vcoordname='soil layer'
1742 l=2
1743 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1744 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1745 ,l,impf,jmpf,nframe,sh2o(1,jsta_2l,l))
1746 if(debugprint)print*,'sample l ',varname,' = ',l,sh2o(im/2,(jsta+jend)/2,l)
1747
1748 varname='sh2o'
1749 vcoordname='soil layer'
1750 l=3
1751 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1752 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1753 ,l,impf,jmpf,nframe,sh2o(1,jsta_2l,l))
1754 if(debugprint)print*,'sample l ',varname,' = ',l,sh2o(im/2,(jsta+jend)/2,l)
1755
1756 varname='sh2o'
1757 vcoordname='soil layer'
1758 l=4
1759 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1760 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1761 ,l,impf,jmpf,nframe,sh2o(1,jsta_2l,l))
1762 if(debugprint)print*,'sample l ',varname,' = ',l,sh2o(im/2,(jsta+jend)/2,l)
1763
1764 varname='smc'
1765 vcoordname='soil layer'
1766 l=1
1767 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1768 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1769 ,l,impf,jmpf,nframe,smc(1,jsta_2l,l))
1770 if(debugprint)print*,'sample l ',varname,' = ',l,smc(im/2,(jsta+jend)/2,l)
1771
1772 varname='smc'
1773 vcoordname='soil layer'
1774 l=2
1775 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1776 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1777 ,l,impf,jmpf,nframe,smc(1,jsta_2l,l))
1778 if(debugprint)print*,'sample l ',varname,' = ',l,smc(im/2,(jsta+jend)/2,l)
1779
1780 varname='smc'
1781 vcoordname='soil layer'
1782 l=3
1783 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1784 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1785 ,l,impf,jmpf,nframe,smc(1,jsta_2l,l))
1786 if(debugprint)print*,'sample l ',varname,' = ',l,smc(im/2,(jsta+jend)/2,l)
1787
1788 varname='smc'
1789 vcoordname='soil layer'
1790 l=4
1791 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1792 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1793 ,l,impf,jmpf,nframe,smc(1,jsta_2l,l))
1794 if(debugprint)print*,'sample l ',varname,' = ',l,smc(im/2,(jsta+jend)/2,l)
1795
1796 varname='stc'
1797 vcoordname='soil layer'
1798 l=1
1799 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1800 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1801 ,l,impf,jmpf,nframe,stc(1,jsta_2l,l))
1802 if(debugprint)print*,'sample l ',varname,' = ',l,stc(im/2,(jsta+jend)/2,l)
1803
1804 varname='stc'
1805 vcoordname='soil layer'
1806 l=2
1807 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1808 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1809 ,l,impf,jmpf,nframe,stc(1,jsta_2l,l))
1810 if(debugprint)print*,'sample l ',varname,' = ',l,stc(im/2,(jsta+jend)/2,l)
1811
1812 varname='stc'
1813 vcoordname='soil layer'
1814 l=3
1815 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1816 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1817 ,l,impf,jmpf,nframe,stc(1,jsta_2l,l))
1818 if(debugprint)print*,'sample l ',varname,' = ',l,stc(im/2,(jsta+jend)/2,l)
1819
1820 varname='stc'
1821 vcoordname='soil layer'
1822 l=4
1823 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1824 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1825 ,l,impf,jmpf,nframe,stc(1,jsta_2l,l))
1826 if(debugprint)print*,'sample l ',varname,' = ',l,stc(im/2,(jsta+jend)/2,l)
1827
1828 varname='pres'
1829 vcoordname='layer'
1830 do l=1,lp1
1831! ll=lp1-l+1
1832 ll=l
1833 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1834 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1835 ,l,impf,jmpf,nframe,pint(1,jsta_2l,ll))
1836 if(debugprint)print*,'sample l ',varname,' = ',ll,pint(im/2,(jsta+jend)/2,ll)
1837 if(l /= 1)then ! assuming post counts from top down
1838 do j=jsta,jend
1839 do i=1,im
1840 alpint(i,j,ll)=alog(pint(i,j,ll))
1841 end do
1842 end do
1843 end if
1844 end do ! do loop for l
1845
1846! do l = 1, lp1
1847 l=1
1848 do j = jsta, jend
1849 do i = 1, im
1850 if(pint(i,j,l) /= 0.0)then
1851 alpint(i,j,l)=alog(pint(i,j,l))
1852 else
1853 alpint(i,j,l)=spval
1854 end if
1855 end do
1856 end do
1857! end do
1858
1859 do l = 2, lp1
1860 do j = jsta_2l, jend_2u
1861 do i = 1, im
1862 pmid(i,j,l-1 ) = (pint(i,j,l-1)+ &
1863 pint(i,j,l))*0.5 ! representative of what model does
1864 end do
1865 end do
1866 if(debugprint)print*,'sample l, PMID = ',l-1,pmid(im/2,(jsta+jend)/2,l-1)
1867 end do
1868
1869 if(gridtype=='E')then
1870 do l = 1, lm
1871 call exch(pmid(1:im,jsta_2l:jend_2u,l))
1872 do j = jsta, jend
1873 do i = 1, im-mod(j,2)
1874 IF(j == 1 .AND. i < im)THEN !SOUTHERN BC
1875 pmidv(i,j,l)=0.5*(pmid(i,j,l)+pmid(i+1,j,l))
1876 ELSE IF(j==jm .AND. i<im)THEN !NORTHERN BC
1877 pmidv(i,j,l)=0.5*(pmid(i,j,l)+pmid(i+1,j,l))
1878 ELSE IF(i == 1 .AND. mod(j,2) == 0) THEN !WESTERN EVEN BC
1879 pmidv(i,j,l)=0.5*(pmid(i,j-1,l)+pmid(i,j+1,l))
1880 ELSE IF(i == im .AND. mod(j,2) == 0 &
1881 .AND. j < jm) THEN !EASTERN EVEN BC
1882 pmidv(i,j,l)=0.5*(pmid(i,j-1,l)+pmid(i,j+1,l))
1883 ELSE IF (mod(j,2) < 1) THEN
1884 pmidv(i,j,l)=0.25*(pmid(i,j,l)+pmid(i-1,j,l) &
1885 +pmid(i,j+1,l)+pmid(i,j-1,l))
1886 ELSE
1887 pmidv(i,j,l)=0.25*(pmid(i,j,l)+pmid(i+1,j,l) &
1888 +pmid(i,j+1,l)+pmid(i,j-1,l))
1889 END IF
1890 end do
1891 end do
1892 end do
1893 else if(gridtype=='B')then
1894 do l = 1, lm
1895 call exch(pmid(1:im,jsta_2l:jend_2u,l))
1896 do j = jsta, jend_m
1897 do i = 1, im-1
1898 pmidv(i,j,l)=0.25*(pmid(i,j,l) + pmid(i+1,j,l) &
1899 + pmid(i,j+1,l) + pmid(i+1,j+1,l))
1900 end do
1901 end do
1902 end do
1903 end if
1904
1905!!!!! COMPUTE Z
1906 allocate(fi(im,jsta_2l:jend_2u,2))
1907 do j = jsta, jend
1908 do i = 1, im
1909 zint(i,j,lm+1) = fis(i,j)/g
1910 if (debugprint .and. i == im/2 .and. j ==(jsta+jend)/2 ) then
1911 write(6,*) 'G,ZINT: ', g,zint(i,j,lm+1)
1912 endif
1913 fi(i,j,1) = fis(i,j)
1914 end do
1915 end do
1916
1917! SECOND, INTEGRATE HEIGHT HYDROSTATICLY
1918 ii = im/2
1919 jj = (jsta+jend)/2
1920 DO l=lm,1,-1
1921!$omp parallel do private(i,j)
1922 do j = jsta, jend
1923 do i = 1, im
1924 fi(i,j,2) = htm(i,j,l)*t(i,j,l)*(q(i,j,l)*d608+1.0)*rd* &
1925 (alpint(i,j,l+1)-alpint(i,j,l))+fi(i,j,1)
1926 zint(i,j,l) = fi(i,j,2)/g
1927 fi(i,j,1) = fi(i,j,2)
1928 ENDDO
1929 ENDDO
1930 if(debugprint)print*,'L,sample HTM,T,Q,ALPINT(L+1),ALPINT(l)', &
1931 ',ZINT= ',l,htm(ii,jj,l),t(ii,jj,l), &
1932 q(ii,jj,l),alpint(ii,jj,l+1), &
1933 alpint(ii,jj,l),zint(ii,jj,l)
1934 END DO
1935 if (me==0) then
1936 print*,'finish deriving geopotential in nmm'
1937 write(*,*)' after ZINT lm=',lm,' js=',js,' je=',je,' im=',im
1938 write(*,*)' zmid lbounds=',lbound(zmid),' ubounds=',ubound(zmid)
1939 write(*,*)' zint lbounds=',lbound(zint),' ubounds=',ubound(zint)
1940 write(*,*)' pmid lbounds=',lbound(pmid),' ubounds=',ubound(pmid)
1941 write(*,*)' pint lbounds=',lbound(pint),' ubounds=',ubound(pint)
1942 endif
1943 deallocate(fi)
1944!
1945 DO l=1,lm
1946! write(*,*)' zmid l=',l
1947!$omp parallel do private(i,j,fact)
1948 DO j=jsta,jend
1949! write(*,*)' zmid j=',j
1950 DO i=1,im
1951! write(*,*)' zmid i=',i
1952! ZMID(I,J,L)=(ZINT(I,J,L+1)+ZINT(I,J,L))*0.5 ! ave of z
1953! write(*,*)' pmid=',pmid(i,j,l)
1954! write(*,*)' pint=',pint(i,j,l),pint(i,j,l+1)
1955! write(*,*)' zint=',zint(i,j,l),zint(i,j,l+1)
1956 fact = (log(pmid(i,j,l))-log(pint(i,j,l))) &
1957 / (log(pint(i,j,l+1))-log(pint(i,j,l)))
1958 zmid(i,j,l) = zint(i,j,l) + (zint(i,j,l+1)-zint(i,j,l))*fact
1959 ENDDO
1960 ENDDO
1961 ENDDO
1962
1963 varname='vvel'
1964 vcoordname='mid layer'
1965 do l=1,lm
1966! ll=lm-l+1
1967 ll=l
1968 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1969 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1970 ,l,impf,jmpf,nframe,wh(1,jsta_2l,ll))
1971 if(debugprint)print*,'sample l ',varname,' = ',ll,wh(im/2,(jsta+jend)/2,ll)
1972 end do ! do loop for l
1973
1974 varname='acfrcv'
1975 vcoordname='sfc'
1976 l=1
1977 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1978 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1979 ,l,impf,jmpf,nframe,acfrcv)
1980 if(debugprint)print*,'sample ',varname,' = ',acfrcv(im/2,(jsta+jend)/2)
1981
1982 varname='acfrst'
1983 vcoordname='sfc'
1984 l=1
1985 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1986 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1987 ,l,impf,jmpf,nframe,acfrst)
1988 if(debugprint)print*,'sample ',varname,' = ',acfrst(im/2,(jsta+jend)/2)
1989
1990!insert-mp
1991 varname='ssroff'
1992 vcoordname='sfc'
1993 l=1
1994 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
1995 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
1996 ,l,impf,jmpf,nframe,ssroff)
1997 if(debugprint)print*,'sample ',varname,' = ',ssroff(im/2,(jsta+jend)/2)
1998
1999! reading UNDERGROUND RUNOFF
2000 varname='bgroff'
2001 vcoordname='sfc'
2002 l=1
2003 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2004 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2005 ,l,impf,jmpf,nframe,bgroff)
2006 if(debugprint)print*,'sample ',varname,' = ',bgroff(im/2,(jsta+jend)/2)
2007
2008 varname='rlwin'
2009 vcoordname='sfc'
2010 l=1
2011 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2012 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2013 ,l,impf,jmpf,nframe,rlwin)
2014 if(debugprint)print*,'sample ',varname,' = ',rlwin(im/2,(jsta+jend)/2)
2015
2016 varname='rlwtoa'
2017 vcoordname='sfc'
2018 l=1
2019 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2020 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2021 ,l,impf,jmpf,nframe,rlwtoa)
2022 if(debugprint)print*,'sample ',varname,' = ',rlwtoa(im/2,(jsta+jend)/2)
2023
2024 varname='alwin'
2025 vcoordname='sfc'
2026 l=1
2027 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2028 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2029 ,l,impf,jmpf,nframe,alwin)
2030 if(debugprint)print*,'sample ',varname,' = ',alwin(im/2,(jsta+jend)/2)
2031
2032 varname='alwout'
2033 vcoordname='sfc'
2034 l=1
2035 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2036 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2037 ,l,impf,jmpf,nframe,alwout)
2038 if(debugprint)print*,'sample ',varname,' = ',alwout(im/2,(jsta+jend)/2)
2039
2040 varname='alwtoa'
2041 vcoordname='sfc'
2042 l=1
2043 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2044 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2045 ,l,impf,jmpf,nframe,alwtoa)
2046 if(debugprint)print*,'sample ',varname,' = ',alwtoa(im/2,(jsta+jend)/2)
2047
2048 varname='rswin'
2049 vcoordname='sfc'
2050 l=1
2051 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2052 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2053 ,l,impf,jmpf,nframe,rswin)
2054 if(debugprint)print*,'sample ',varname,' = ',rswin(im/2,(jsta+jend)/2)
2055
2056 varname='rswinc'
2057 vcoordname='sfc'
2058 l=1
2059 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2060 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2061 ,l,impf,jmpf,nframe,rswinc)
2062 if(debugprint)print*,'sample ',varname,' = ',rswinc(im/2,(jsta+jend)/2)
2063
2064 varname='rswout'
2065 vcoordname='sfc'
2066 l=1
2067 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2068 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2069 ,l,impf,jmpf,nframe,rswout)
2070 if(debugprint)print*,'sample ',varname,' = ',rswout(im/2,(jsta+jend)/2)
2071
2072 varname='aswin'
2073 vcoordname='sfc'
2074 l=1
2075 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2076 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2077 ,l,impf,jmpf,nframe,aswin)
2078 if(debugprint)print*,'sample ',varname,' = ',aswin(im/2,(jsta+jend)/2)
2079
2080 varname='aswout'
2081 vcoordname='sfc'
2082 l=1
2083 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2084 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2085 ,l,impf,jmpf,nframe,aswout)
2086 if(debugprint)print*,'sample ',varname,' = ',aswout(im/2,(jsta+jend)/2)
2087
2088 varname='aswtoa'
2089 vcoordname='sfc'
2090 l=1
2091 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2092 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2093 ,l,impf,jmpf,nframe,aswtoa)
2094 if(debugprint)print*,'sample ',varname,' = ',aswtoa(im/2,(jsta+jend)/2)
2095
2096 varname='sfcshx'
2097 vcoordname='sfc'
2098 l=1
2099 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2100 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2101 ,l,impf,jmpf,nframe,sfcshx)
2102 if(debugprint)print*,'sample ',varname,' = ',sfcshx(im/2,(jsta+jend)/2)
2103
2104 varname='sfclhx'
2105 vcoordname='sfc'
2106 l=1
2107 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2108 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2109 ,l,impf,jmpf,nframe,sfclhx)
2110 if(debugprint)print*,'sample ',varname,' = ',sfclhx(im/2,(jsta+jend)/2)
2111
2112 varname='subshx'
2113 vcoordname='sfc'
2114 l=1
2115 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2116 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2117 ,l,impf,jmpf,nframe,subshx)
2118 if(debugprint)print*,'sample ',varname,' = ',subshx(im/2,(jsta+jend)/2)
2119
2120 varname='snopcx'
2121 vcoordname='sfc'
2122 l=1
2123 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2124 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2125 ,l,impf,jmpf,nframe,snopcx)
2126 if(debugprint)print*,'sample ',varname,' = ',snopcx(im/2,(jsta+jend)/2)
2127
2128 varname='sfcuvx'
2129 vcoordname='sfc'
2130 l=1
2131 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2132 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2133 ,l,impf,jmpf,nframe,sfcuvx)
2134 if(debugprint)print*,'sample ',varname,' = ',sfcuvx(im/2,(jsta+jend)/2)
2135
2136 varname='potevp'
2137 vcoordname='sfc'
2138 l=1
2139 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2140 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2141 ,l,impf,jmpf,nframe,potevp)
2142 if(debugprint)print*,'sample ',varname,' = ',potevp(im/2,(jsta+jend)/2)
2143
2144 varname='rlwtt'
2145 vcoordname='mid layer'
2146 do l=1,lm
2147! ll=lm-l+1
2148 ll=l
2149 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2150 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2151 ,l,impf,jmpf,nframe,rlwtt(1,jsta_2l,ll))
2152 if(debugprint)print*,'sample l ',varname,' = ',ll,rlwtt(im/2,(jsta+jend)/2,ll)
2153 end do ! do loop for l
2154
2155 varname='rswtt'
2156 vcoordname='mid layer'
2157 do l=1,lm
2158! ll=lm-l+1
2159 ll=l
2160 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2161 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2162 ,l,impf,jmpf,nframe,rswtt(1,jsta_2l,ll))
2163 if(debugprint)print*,'sample l ',varname,' = ',ll,rswtt(im/2,(jsta+jend)/2,ll)
2164 end do ! do loop for l
2165 where(rlwtt/=spval .and. rswtt/=spval)ttnd=rswtt+rlwtt
2166
2167 varname='tcucn'
2168 vcoordname='mid layer'
2169 do l=1,lm
2170! ll=lm-l+1
2171 ll=l
2172 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2173 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2174 ,l,impf,jmpf,nframe,tcucn(1,jsta_2l,ll))
2175 if(debugprint)print*,'sample l ',varname,' = ',ll,tcucn(im/2,(jsta+jend)/2,ll)
2176 end do ! do loop for l
2177
2178 varname='train'
2179 vcoordname='mid layer'
2180 do l=1,lm
2181! ll=lm-l+1
2182 ll=l
2183 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2184 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2185 ,l,impf,jmpf,nframe,train(1,jsta_2l,ll))
2186 if(debugprint)print*,'sample l ',varname,' = ',ll,train(im/2,(jsta+jend)/2,ll)
2187 end do ! do loop for l
2188
2189 varname='CFRCV'
2190 vcoordname='sfc'
2191 l=1
2192 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2193 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2194 ,l,impf,jmpf,nframe,ncfrcv)
2195 if(debugprint)print*,'sample ',varname,' = ',ncfrcv(im/2,(jsta+jend)/2)
2196
2197 varname='CFRST'
2198 vcoordname='sfc'
2199 l=1
2200 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2201 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2202 ,l,impf,jmpf,nframe,ncfrst)
2203 if(debugprint)print*,'sample ',varname,' = ',ncfrst(im/2,(jsta+jend)/2)
2204
2205! set default to not empty buket
2206 nsrfc=0
2207 nrdlw=0
2208 nrdsw=0
2209 nheat=0
2210 nclod=0
2211 nprec=0
2212 nphs=0
2213
2214 varname='nprec'
2215 if(me==0)then
2216 call nemsio_getheadvar(nfile,trim(varname),nprec,iret)
2217 if (iret /= 0) then
2218 print*,varname," not found in file-Assigned zero"
2219 end if
2220 end if
2221 call mpi_bcast(nprec,1,mpi_integer,0,mpi_comm_comp,iret)
2222 if(debugprint)print*,'sample ',varname,' = ',nprec
2223
2224 varname='nphs'
2225 if(me==0)then
2226 call nemsio_getheadvar(nfile,trim(varname),nphs,iret)
2227 if (iret /= 0) then
2228 print*,varname," not found in file-Assigned zero"
2229 end if
2230 end if
2231 call mpi_bcast(nphs,1,mpi_integer,0,mpi_comm_comp,iret)
2232 if(debugprint)print*,'sample ',varname,' = ',nphs
2233
2234 varname='nclod'
2235 if(me==0)then
2236 call nemsio_getheadvar(nfile,trim(varname),nclod,iret)
2237 if (iret /= 0) then
2238 print*,varname," not found in file-Assigned zero"
2239 end if
2240 end if
2241 call mpi_bcast(nclod,1,mpi_integer,0,mpi_comm_comp,iret)
2242 if(debugprint)print*,'sample ',varname,' = ',nclod
2243
2244 varname='nheat'
2245 if(me==0)then
2246 call nemsio_getheadvar(nfile,trim(varname),nheat,iret)
2247 if (iret /= 0) then
2248 print*,varname," not found in file-Assigned zero"
2249 end if
2250 end if
2251 call mpi_bcast(nheat,1,mpi_integer,0,mpi_comm_comp,iret)
2252 if(debugprint)print*,'sample ',varname,' = ',nheat
2253
2254 varname='nrdlw'
2255 if(me==0)then
2256 call nemsio_getheadvar(nfile,trim(varname),nrdlw,iret)
2257 if (iret /= 0) then
2258 print*,varname," not found in file-Assigned zero"
2259 end if
2260 end if
2261 call mpi_bcast(nrdlw,1,mpi_integer,0,mpi_comm_comp,iret)
2262 if(debugprint)print*,'sample ',varname,' = ',nrdlw
2263
2264 varname='nrdsw'
2265 if(me==0)then
2266 call nemsio_getheadvar(nfile,trim(varname),nrdsw,iret)
2267 if (iret /= 0) then
2268 print*,varname," not found in file-Assigned zero"
2269 end if
2270 end if
2271 call mpi_bcast(nrdsw,1,mpi_integer,0,mpi_comm_comp,iret)
2272 if(debugprint)print*,'sample ',varname,' = ',nrdsw
2273
2274 varname='nsrfc'
2275 if(me==0)then
2276 call nemsio_getheadvar(nfile,trim(varname),nsrfc,iret)
2277 if (iret /= 0) then
2278 print*,varname," not found in file-Assigned zero"
2279 end if
2280 end if
2281 call mpi_bcast(nsrfc,1,mpi_integer,0,mpi_comm_comp,iret)
2282 if(debugprint)print*,'sample ',varname,' = ',nsrfc
2283
2284! VarName='avrain'
2285! if(me==0)then
2286! call nemsio_getheadvar(nfile,trim(varname),avrain,iret)
2287! if (iret /= 0) then
2288! print*,VarName," not found in file-Assigned zero"
2289! end if
2290! end if
2291! call mpi_bcast(avrain,1,MPI_REAL,0,mpi_comm_comp,iret)
2292! if(debugprint)print*,'sample ',VarName,' = ',avrain
2293
2294! VarName='avcnvc'
2295! if(me==0)then
2296! call nemsio_getheadvar(nfile,trim(varname),avcnvc,iret)
2297! if (iret /= 0) then
2298! print*,VarName," not found in file-Assigned zero"
2299! end if
2300! end if
2301! call mpi_bcast(avcnvc,1,MPI_REAL,0,mpi_comm_comp,iret)
2302! if(debugprint)print*,'sample ',VarName,' = ',avcnvc
2303
2304! VarName='ardlw'
2305! if(me==0)then
2306! call nemsio_getheadvar(nfile,trim(varname),ardlw,iret)
2307! if (iret /= 0) then
2308! print*,VarName," not found in file-Assigned zero"
2309! end if
2310! end if
2311! call mpi_bcast(ardlw,1,MPI_REAL,0,mpi_comm_comp,iret)
2312! if(debugprint)print*,'sample ',VarName,' = ',ardlw
2313
2314! VarName='ardsw'
2315! if(me==0)then
2316! call nemsio_getheadvar(nfile,trim(varname),ardsw,iret)
2317! if (iret /= 0) then
2318! print*,VarName," not found in file-Assigned zero"
2319! end if
2320! end if
2321! call mpi_bcast(ardsw,1,MPI_REAL,0,mpi_comm_comp,iret)
2322! if(debugprint)print*,'sample ',VarName,' = ',ardsw
2323
2324! VarName='asrfc'
2325! if(me==0)then
2326! call nemsio_getheadvar(nfile,trim(varname),asrfc,iret)
2327! if (iret /= 0) then
2328! print*,VarName," not found in file-Assigned zero"
2329! end if
2330! end if
2331! call mpi_bcast(asrfc,1,MPI_REAL,0,mpi_comm_comp,iret)
2332! if(debugprint)print*,'sample ',VarName,' = ',asrfc
2333
2334! ! setting all counters to 1 because Ratko has divided all fluxes by counters within NEMS model
2335! avrain=1.0
2336! avcnvc=1.0
2337! ardlw=1.0
2338! ardsw=1.0
2339! asrfc=1.0
2340! if(debugprint)print*,'sample avrain, avcnvc, ardlw, ardsw, asrfc = ', &
2341! avrain, avcnvc, ardlw, ardsw, asrfc
2342
2343!-- Changes to NMMB to allow for counters to vary during the forecast by making them
2344! 2D arrays in order to get around an ESMF limitation (Ferrier 13 Aug 2009)
2345
2346 varname='AVRAIN'
2347 vcoordname='sfc'
2348 l=1
2349 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2350 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2351 ,l,impf,jmpf,nframe,buf)
2352 avrain=buf(im/2,(jsta+jend)/2)
2353 if(debugprint)print*,'sample ',varname,' = ',avrain
2354
2355 varname='AVCNVC'
2356 vcoordname='sfc'
2357 l=1
2358 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2359 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2360 ,l,impf,jmpf,nframe,buf)
2361 avcnvc=buf(im/2,(jsta+jend)/2)
2362 if(debugprint)print*,'sample ',varname,' = ',avcnvc
2363
2364 varname='ARDLW'
2365 vcoordname='sfc'
2366 l=1
2367 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2368 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2369 ,l,impf,jmpf,nframe,buf)
2370 ardlw=buf(im/2,(jsta+jend)/2)
2371 if(debugprint)print*,'sample ',varname,' = ',ardlw
2372
2373 varname='ARDSW'
2374 vcoordname='sfc'
2375 l=1
2376 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2377 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2378 ,l,impf,jmpf,nframe,buf)
2379 ardsw=buf(im/2,(jsta+jend)/2)
2380 if(debugprint)print*,'sample ',varname,' = ',ardsw
2381
2382 varname='ASRFC'
2383 vcoordname='sfc'
2384 l=1
2385 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2386 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2387 ,l,impf,jmpf,nframe,buf)
2388 asrfc=buf(im/2,(jsta+jend)/2)
2389 if(debugprint)print*,'sample ',varname,' = ',asrfc
2390
2391! reading 10 m wind
2392 varname='u10'
2393 vcoordname='10 m above gnd'
2394 l=1
2395 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2396 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2397 ,l,impf,jmpf,nframe,u10h)
2398! Chuang Aug 2012: 10 m winds are computed on mass points in the model
2399! post interpolates them onto V points because copygb interpolates
2400! wind points differently and 10 m winds are identified as 33/34
2401 call h2u(u10h,u10)
2402 if(debugprint)print*,'sample ',varname,' = ',u10(im/2,(jsta+jend)/2)
2403
2404 varname='v10'
2405 vcoordname='10 m above gnd'
2406 l=1
2407 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2408 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2409 ,l,impf,jmpf,nframe,v10h)
2410 call h2u(v10h,v10)
2411 if(debugprint)print*,'sample ',varname,' = ',v10(im/2,(jsta+jend)/2)
2412
2413 varname='u10max'
2414 vcoordname='10 m above gnd'
2415 l=1
2416 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2417 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2418 ,l,impf,jmpf,nframe,u10max)
2419 if(debugprint)print*,'sample ',varname,' = ',u10max(im/2,(jsta+jend)/2)
2420
2421 varname='v10max'
2422 vcoordname='10 m above gnd'
2423 l=1
2424 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2425 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2426 ,l,impf,jmpf,nframe,v10max)
2427 if(debugprint)print*,'sample ',varname,' = ',v10max(im/2,(jsta+jend)/2)
2428
2429! reading SMSTAV
2430 varname='smstav'
2431 vcoordname='sfc'
2432 l=1
2433 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2434 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2435 ,l,impf,jmpf,nframe,smstav)
2436 if(debugprint)print*,'sample ',varname,' = ',smstav(im/2,(jsta+jend)/2)
2437 if(debugprint)print*,'MAX/MIN ',varname,' = ' &
2438 ,maxval(smstav),minval(smstav)
2439
2440 varname='smstot'
2441 vcoordname='sfc'
2442 l=1
2443 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2444 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2445 ,l,impf,jmpf,nframe,smstot)
2446 if(debugprint)print*,'sample ',varname,' = ',smstot(im/2,(jsta+jend)/2)
2447
2448! reading VEGETATION TYPE
2449 varname='VGTYP'
2450 vcoordname='sfc'
2451 l=1
2452 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2453 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2454 ,l,impf,jmpf,nframe,sfcevp) ! temporary use sfcevp because it's real in nemsio
2455
2456! do j=jsta,jend
2457! do i=1,im
2458! if(sfcevp(i,j)> 27.0 .or. sfcevp(i,j)<1.0)print*, &
2459! 'bad vegtype=',i,j,sfcevp(i,j)
2460! end do
2461! end do
2462
2463 where(sfcevp /= spval)ivgtyp=nint(sfcevp)
2464 if(debugprint)print*,'sample ',varname,' = ',ivgtyp(im/2,(jsta+jend)/2)
2465
2466 varname='SLTYP'
2467 vcoordname='sfc'
2468 l=1
2469 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2470 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2471 ,l,impf,jmpf,nframe,sfcevp)
2472 where(sfcevp /= spval)isltyp=nint(sfcevp)
2473 if(debugprint)print*,'sample ',varname,' = ',isltyp(im/2,(jsta+jend)/2)
2474
2475! call retrieve_index(index,VarName,varname_all,nrecs,iret)
2476! if (iret /= 0) then
2477! print*,VarName," not found in file-Assigned missing values"
2478! ISLTYP=NINT(SPVAL)
2479! else
2480! this_offset=file_offset(index+1)+(jsta_2l-1)*4*im
2481! this_length=im*(jend_2u-jsta_2l+1)
2482! call mpi_file_read_at(iunit,this_offset &
2483! ,isltyp,this_length,mpi_integer4, mpi_status_ignore, ierr)
2484! if (ierr /= 0) then
2485! print*,"Error reading ", VarName,"Assigned missing values"
2486! ISLTYP=NINT(SPVAL)
2487! end if
2488! end if
2489
2490 varname='sfcevp'
2491 vcoordname='sfc'
2492 l=1
2493 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2494 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2495 ,l,impf,jmpf,nframe,sfcevp)
2496 if(debugprint)print*,'sample ',varname,' = ',sfcevp(im/2,(jsta+jend)/2)
2497
2498 varname='sfcexc'
2499 vcoordname='sfc'
2500 l=1
2501 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2502 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2503 ,l,impf,jmpf,nframe,sfcexc)
2504 if(debugprint)print*,'sample ',varname,' = ',sfcexc(im/2,(jsta+jend)/2)
2505 if(debugprint)print*,'MAX/MIN ',varname,' = ' &
2506 ,maxval(sfcexc),minval(sfcexc)
2507
2508 varname='acsnow'
2509 vcoordname='sfc'
2510 l=1
2511 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2512 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2513 ,l,impf,jmpf,nframe,acsnow)
2514 if(debugprint)print*,'sample ',varname,' = ',acsnow(im/2,(jsta+jend)/2)
2515
2516 varname='acsnom'
2517 vcoordname='sfc'
2518 l=1
2519 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2520 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2521 ,l,impf,jmpf,nframe,acsnom)
2522 if(debugprint)print*,'sample ',varname,' = ',acsnom(im/2,(jsta+jend)/2)
2523
2524 varname='tsea'
2525 vcoordname='sfc'
2526 l=1
2527 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2528 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2529 ,l,impf,jmpf,nframe,sst)
2530 if(debugprint)print*,'sample ',varname,' = ',sst(im/2,(jsta+jend)/2)
2531
2532! VarName='EL_PBL' ! not in nems io yet
2533 varname='xlen_mix'
2534 vcoordname='mid layer'
2535 do l=1,lm
2536! ll=lm-l+1
2537 ll=l
2538 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2539 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2540 ,l,impf,jmpf,nframe,el_pbl(1,jsta_2l,ll))
2541 if(debugprint)print*,'sample l ',varname,' = ',ll,el_pbl(im/2,(jsta+jend)/2,ll)
2542 end do ! do loop for l
2543
2544 varname='exch_h'
2545 vcoordname='mid layer'
2546 do l=1,lm
2547! ll=lm-l+1
2548 ll=l
2549 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2550 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2551 ,l,impf,jmpf,nframe,exch_h(1,jsta_2l,ll))
2552 if(debugprint)print*,'sample l ',varname,' = ',ll,exch_h(im/2,(jsta+jend)/2,ll)
2553 end do ! do loop for l
2554
2555 varname='thz0'
2556 vcoordname='sfc'
2557 l=1
2558 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2559 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2560 ,l,impf,jmpf,nframe,thz0)
2561 if(debugprint)print*,'sample ',varname,' = ',thz0(im/2,(jsta+jend)/2)
2562
2563 varname='qz0'
2564 vcoordname='sfc'
2565 l=1
2566 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2567 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2568 ,l,impf,jmpf,nframe,qz0)
2569 if(debugprint)print*,'sample ',varname,' = ',qz0(im/2,(jsta+jend)/2)
2570
2571 varname='uz0'
2572 vcoordname='sfc'
2573 l=1
2574 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2575 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2576 ,l,impf,jmpf,nframe,uz0)
2577 if(debugprint)print*,'sample ',varname,' = ',uz0(im/2,(jsta+jend)/2)
2578
2579 varname='vz0'
2580 vcoordname='sfc'
2581 l=1
2582 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2583 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2584 ,l,impf,jmpf,nframe,vz0)
2585 if(debugprint)print*,'sample ',varname,' = ',vz0(im/2,(jsta+jend)/2)
2586
2587!
2588! Very confusing story ...
2589!
2590! Retrieve htop and hbot => They are named CNVTOP, CNVBOT in the model and
2591! with HBOTS,HTOPS (shallow conv) and HBOTD,HTOPD (deep conv) represent
2592! the 3 sets of convective cloud base/top arrays tied to the frequency
2593! that history files are written.
2594!
2595! IN THE *MODEL*, arrays HBOT,HTOP are similar to CNVTOP,CNVBOT but are
2596! used in radiation and are tied to the frequency of radiation updates.
2597!
2598! For historical reasons model arrays CNVTOP,CNVBOT are renamed HBOT,HTOP
2599! and manipulated throughout the post.
2600
2601! retrieve htop and hbot
2602! VarName='HTOP'
2603 varname='cnvtop'
2604 vcoordname='sfc'
2605 l=1
2606 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2607 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2608 ,l,impf,jmpf,nframe,htop)
2609 where(htop /= spval)htop=float(lm)-htop+1.0
2610! where(htop /= spval .and. htop > lm)htop=lm*1.0
2611 if(debugprint)print*,'sample ',varname,' = ',htop(im/2,(jsta+jend)/2)
2612
2613! VarName='HBOT'
2614 varname='cnvbot'
2615 vcoordname='sfc'
2616 l=1
2617 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2618 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2619 ,l,impf,jmpf,nframe,hbot)
2620 where(hbot /= spval)hbot=float(lm)-hbot+1.0
2621! where(hbot /= spval .and. hbot > lm)hbot=lm*1.0
2622 if(debugprint)print*,'sample ',varname,' = ',hbot(im/2,(jsta+jend)/2)
2623
2624 varname='htopd'
2625 vcoordname='sfc'
2626 l=1
2627 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2628 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2629 ,l,impf,jmpf,nframe,htopd)
2630 where(htopd /= spval)htopd=float(lm)-htopd+1.0
2631! where(htopd /= spval .and. htopd > lm)htopd=lm*1.0
2632 if(debugprint)print*,'sample ',varname,' = ',htopd(im/2,(jsta+jend)/2)
2633
2634 varname='hbotd'
2635 vcoordname='sfc'
2636 l=1
2637 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2638 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2639 ,l,impf,jmpf,nframe,hbotd)
2640 where(hbotd /= spval)hbotd=float(lm)-hbotd+1.0
2641! where(hbotd /= spval .and. hbotd > lm)hbotd=lm*1.0
2642 if(debugprint)print*,'sample ',varname,' = ',hbotd(im/2,(jsta+jend)/2)
2643
2644 varname='htops'
2645 vcoordname='sfc'
2646 l=1
2647 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2648 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2649 ,l,impf,jmpf,nframe,htops)
2650 where(htops /= spval)htops=float(lm)-htops+1.0
2651! where(htops /= spval .and. htops > lm)htops=lm*1.0
2652 if(debugprint)print*,'sample ',varname,' = ',htops(im/2,(jsta+jend)/2)
2653
2654 varname='hbots'
2655 vcoordname='sfc'
2656 l=1
2657 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2658 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2659 ,l,impf,jmpf,nframe,hbots)
2660 where(hbots /= spval)hbots=float(lm)-hbots+1.0
2661! where(hbots /= spval .and. hbots > lm)hbots=lm*1.0
2662 if(debugprint)print*,'sample ',varname,' = ',hbots(im/2,(jsta+jend)/2)
2663
2664 varname='cuppt'
2665 vcoordname='sfc'
2666 l=1
2667 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2668 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2669 ,l,impf,jmpf,nframe,cuppt)
2670 if(debugprint)print*,'sample ',varname,' = ',cuppt(im/2,(jsta+jend)/2)
2671
2672 varname='cprate'
2673 vcoordname='sfc'
2674 l=1
2675 call getnemsandscatter(me,nfile,im,jm,jsta,jsta_2l &
2676 ,jend_2u,mpi_comm_comp,icnt,idsp,spval,varname,vcoordname &
2677 ,l,impf,jmpf,nframe,cprate)
2678 if(debugprint)print*,'sample ',varname,' = ',cprate(im/2,(jsta+jend)/2)
2679
2680!!!! DONE GETTING
2681
2682!$omp parallel do private(i,j,l)
2683 do l = 1, lm
2684 do j = jsta, jend
2685 do i = 1, im
2686 IF(abs(t(i,j,l))>1.0e-3 .and. (wh(i,j,1) < spval)) &
2687 omga(i,j,l) = -wh(i,j,l)*pmid(i,j,l)*g &
2688 / (rd*t(i,j,l)*(1.+d608*q(i,j,l)))
2689 end do
2690 end do
2691 end do
2692
2693 thl=210.
2694 plq=70000.
2695
2696 CALL table(ptbl,ttbl,pt, &
2697 rdq,rdth,rdp,rdthe,pl,thl,qs0,sqs,sthe,the0)
2698
2699 CALL tableq(ttblq,rdpq,rdtheq,plq,thl,stheq,the0q)
2700
2701!
2702!
2703 IF(me==0)THEN
2704 WRITE(6,*)' SPL (POSTED PRESSURE LEVELS) BELOW: '
2705 WRITE(6,51) (spl(l),l=1,lsm)
2706 50 FORMAT(14(f4.1,1x))
2707 51 FORMAT(8(f8.1,1x))
2708 ENDIF
2709!
2710! COMPUTE DERIVED TIME STEPPING CONSTANTS.
2711!
2712!MEB need to get DT
2713! DT = 120. !MEB need to get DT
2714! NPHS = 4 !MEB need to get physics DT
2715 dtq2 = dt * nphs !MEB need to get physics DT
2716 tsph = 3600./dt !MEB need to get DT
2717
2718 IF (pthresh > 0.) THEN
2719 pthresh = 0.01*dtq2/3.6e6 !-- Precip rate >= 0.01 mm/h
2720! PTHRESH = 0.01*DTQ2/(3600.*39.37) !-- Precip rate >= 0.01 inches/h
2721 ENDIF
2722
2723 tsrfc=float(nsrfc)/tsph
2724 if(me==0)write(6,*)'tsfrc ',tsrfc,nsrfc,tsph
2725 IF(nsrfc==0)tsrfc=float(ifhr) !in case buket does not get emptied
2726 trdlw=float(nrdlw)/tsph
2727 IF(nrdlw==0)trdlw=float(ifhr) !in case buket does not get emptied
2728 trdsw=float(nrdsw)/tsph
2729 IF(nrdsw==0)trdsw=float(ifhr) !in case buket does not get emptied
2730 theat=float(nheat)/tsph
2731 IF(nheat==0)theat=float(ifhr) !in case buket does not get emptied
2732 tclod=float(nclod)/tsph
2733 IF(nclod==0)tclod=float(ifhr) !in case buket does not get emptied
2734 tprec=float(nprec)/tsph
2735 IF(nprec==0)tprec=float(ifhr) !in case buket does not get emptied
2736! TPREC=float(ifhr)
2737 if(me==0)print*,'TSRFC TRDLW TRDSW THEAT TCLOD TPREC= ' &
2738 ,tsrfc, trdlw, trdsw, theat, tclod, tprec
2739!MEB need to get DT
2740
2741!how am i going to get this information?
2742! NPREC = INT(TPREC *TSPH+D50)
2743! NHEAT = INT(THEAT *TSPH+D50)
2744! NCLOD = INT(TCLOD *TSPH+D50)
2745! NRDSW = INT(TRDSW *TSPH+D50)
2746! NRDLW = INT(TRDLW *TSPH+D50)
2747! NSRFC = INT(TSRFC *TSPH+D50)
2748!how am i going to get this information?
2749!
2750! IF(ME==0)THEN
2751! WRITE(6,*)' '
2752! WRITE(6,*)'DERIVED TIME STEPPING CONSTANTS'
2753! WRITE(6,*)' NPREC,NHEAT,NSRFC : ',NPREC,NHEAT,NSRFC
2754! WRITE(6,*)' NCLOD,NRDSW,NRDLW : ',NCLOD,NRDSW,NRDLW
2755! ENDIF
2756!
2757! COMPUTE DERIVED MAP OUTPUT CONSTANTS.
2758 DO l = 1,lsm
2759 alsl(l) = alog(spl(l))
2760 END DO
2761!
2762!HC WRITE IGDS OUT FOR WEIGHTMAKER TO READ IN AS KGDSIN
2763 if(me==0)then
2764 print*,'writing out igds'
2765 igdout=110
2766! open(igdout,file='griddef.out',form='unformatted'
2767! + ,status='unknown')
2768 IF(maptype==203)THEN !A STAGGERED E-GRID
2769 WRITE(igdout)203
2770 WRITE(igdout)im
2771 WRITE(igdout)jm
2772 WRITE(igdout)latstart
2773 WRITE(igdout)lonstart
2774 WRITE(igdout)136
2775 WRITE(igdout)cenlat
2776 WRITE(igdout)cenlon
2777 WRITE(igdout)dxval
2778 WRITE(igdout)dyval
2779 WRITE(igdout)64
2780 WRITE(igdout)0
2781 WRITE(igdout)0
2782 WRITE(igdout)0
2783 WRITE(igdout)latlast
2784 WRITE(igdout)lonlast
2785 ELSE IF(maptype==205)THEN !A STAGGERED B-GRID
2786 WRITE(igdout)205
2787 WRITE(igdout)im
2788 WRITE(igdout)jm
2789 WRITE(igdout)latstart
2790 WRITE(igdout)lonstart
2791 WRITE(igdout)136
2792 WRITE(igdout)cenlat
2793 WRITE(igdout)cenlon
2794 WRITE(igdout)dxval
2795 WRITE(igdout)dyval
2796 WRITE(igdout)64
2797 WRITE(igdout)latlast
2798 WRITE(igdout)lonlast
2799 WRITE(igdout)0
2800 WRITE(igdout)0
2801 WRITE(igdout)0
2802 END IF
2803 open(111,file='copygb_gridnav.txt',form='formatted' &
2804 ,status='unknown')
2805 IF(maptype==203)THEN !A STAGGERED E-GRID
2806 write(111,1000) 2*im-1,jm,latstart,lonstart,cenlon, &
2807 nint(dxval*107.),nint(dyval*110.),cenlat,cenlat
2808 ELSE IF(maptype==205)THEN !A STAGGERED B-GRID
2809 if(grib=="grib2") then
2810 write(111,1000) im,jm,latstart/1000,lonstart/1000,cenlon/1000, &
2811 nint(dxval*107.)/1000,nint(dyval*110.)/1000, &
2812 cenlat/1000,cenlat/1000, &
2813 latlast/1000,lonlast/1000
2814 endif
2815 END IF
28161000 format('255 3 ',2(i4,x),i6,x,i7,x,'8 ',i7,x,2(i6,x),'0 64', &
2817 3(x,i6),x,i7)
2818 close(111)
2819!
2820 IF (maptype==205)THEN !A STAGGERED B-GRID
2821 open(112,file='latlons_corners.txt',form='formatted' &
2822 ,status='unknown')
2823 if(grib=="grib2") then
2824 write(112,1001)latstart/1000,(lonstart/1000)-360000, &
2825 latse/1000, &
2826 lonse/1000,latnw/1000,lonnw/1000,latlast/1000, &
2827 (lonlast/1000)-360000
2828 endif
28291001 format(4(i6,x,i7,x))
2830 close(112)
2831 ENDIF
2832
2833 end if
2834
2835! close all files
2836 call nemsio_close(nfile,iret=status)
2837!
2838 if(me==0)write(*,*)'end of INIT_NEMS'
2839
2840 RETURN
2841 END
collect_loc
subroutine collect_loc(a, b)
COLLECT_LOC.
Definition COLLECT_LOC.f:23
exch
subroutine exch(a)
exch() Subroutine that exchanges one halo row.
Definition EXCH.f:27
initpost_nems
subroutine initpost_nems(nrec, nfile)
INITPOST_NEMS This routine initializes constants and variables at the start of an NEMS model or post ...
Definition INITPOST_NEMS.f:28