@@ -272,6 +272,9 @@ program trakmain
272272c wrapping was being handled. I fixed the
273273c GM-wrapping in a few different areas in both
274274c subroutines.
275+ c
276+ c 23-07-31 Marchok Fixed an error in output_atcfunix where I had
277+ c an extra comma in the output record.
275278c
276279c
277280c Input files:
@@ -9046,7 +9049,7 @@ subroutine get_gen_diags (imax,jmax,inp,dx,dy
90469049 real fixlon(maxstorm,maxtime),fixlat(maxstorm,maxtime)
90479050 real clon(maxstorm,maxtime,maxtp)
90489051 real clat(maxstorm,maxtime,maxtp)
9049- real dx,dy,xcenlon,xcenlat,q850conv, q850_smooth
9052+ real dx,dy,xcenlon,xcenlat,q850_smooth
90509053 real rh_1000_925_smooth,rh_800_600_smooth,omega500_smooth
90519054 real divg,moist_divg,re,ri,xsmoothval
90529055 character :: already_computed_domain_wide_rh*1
@@ -9087,7 +9090,8 @@ subroutine get_gen_diags (imax,jmax,inp,dx,dy
90879090 !----------------------------------------------------------------
90889091 ! Now get a smoothed, barnes-averaged value of q850 at the center
90899092 ! point. Then multiply the 850 mb divg we just calculated by the
9090- ! smoothed q850 to get the 850 mb moisture convergence (q850conv).
9093+ ! smoothed q850 to get the 850 mb moisture convergence
9094+ ! (moist_divg).
90919095 !----------------------------------------------------------------
90929096
90939097 if (readgenflag(1)) then
@@ -9119,14 +9123,17 @@ subroutine get_gen_diags (imax,jmax,inp,dx,dy
91199123 write (6,126) date_time(5),date_time(6),date_time(7)
91209124 126 format (1x,'TIMING: gen_diag after smooth q850 ... ',i2.2,':'
91219125 & ,i2.2,':',i2.2)
9126+ write (6,228) igsvret,q850_smooth,igdret
9127+ 228 format (1x,' After get_smooth_value_at_pt for q850, igsvret= '
9128+ & ,i4,' q850_smooth= ',f12.8,' igdret= ',i4)
91229129 endif
91239130
91249131 endif
91259132
91269133 if (igdret == 0 .and. igsvret == 0) then
9127- q850conv = divg * q850_smooth
9134+ moist_divg = divg * q850_smooth
91289135 else
9129- q850conv = -9999.0
9136+ moist_divg = -9999.0
91309137 endif
91319138
91329139 !----------------------------------------------------------------
@@ -9141,6 +9148,9 @@ subroutine get_gen_diags (imax,jmax,inp,dx,dy
91419148 write (6,128) date_time(5),date_time(6),date_time(7)
91429149 128 format (1x,'TIMING: gen_diag before get_rh ... ',i2.2,':'
91439150 & ,i2.2,':',i2.2)
9151+ write (6,230) divg,q850_smooth,moist_divg
9152+ 230 format (1x,' divg= ',f14.10,' q850_smooth= ',f14.10
9153+ & ,' moist_divg= ',f14.10)
91449154 endif
91459155
91469156 igrhret = 0
@@ -11101,14 +11111,15 @@ subroutine output_aext (outlon,outlat,inp,ist
1110111111 & ,2(', ',i4),', ',i3,a27,2(', ',i3),a44
1110211112 & ,', THERMO PARAMS'
1110311113 & ,3(', ',i7),', ',a1,', ',i2,', DT, -999, SHR82, ',i4,', '
11104- & ,i3,', SST, ',i4,', ARMW',2(', ',i3),', ',i8,8(', ',i5))
11114+ & ,i3,', SST, ',i4,', ARMW',2(', ',i3),2(', ',i9)
11115+ & ,7(', ',i5))
1110511116 91 format (a2,', ',a4,', ',i10.10,', 03, ',a4,', ',i3.3,', ',i3,a1
1110611117 & ,', ',i4,a1,', ',i3,', ',i4,', ',a12,4(', ',i4.4)
1110711118 & ,2(', ',i4),', ',i3,a27,2(', ',i3),a44
1110811119 & ,', THERMO PARAMS'
1110911120 & ,3(', ',i7),', ',a1,', ',i2,', DT, -999, SHR82, ',i4,', '
11110- & ,i3,', SST, ',i4,', ARMW',2(', ',i3),', ',i8,8 (', ',i5 )
11111- & ,', ',a3)
11121+ & ,i3,', SST, ',i4,', ARMW',2(', ',i3),2 (', ',i9 )
11122+ & ,7(', ',i5), ', ',a3)
1111211123
1111311124c bug fix for IBM: flush the output stream so it actually writes
1111411125 flush(68)
@@ -12953,16 +12964,24 @@ subroutine output_atcf_gen (outlon,outlat,inp,ist
1295312964
1295412965 if (genflag == 'y' .or. genflag == 'Y') then
1295512966
12956- if (divg > -998.0) then
12957- idivg = int ((divg * 1e6) + 0.5)
12958- else
12967+ write (6,125) divg
12968+ 125 format (1x,' in output_atcf_gen, before scaling, divg= ',f16.8)
12969+ write (6,127) moist_divg
12970+ 127 format (1x,' in output_atcf_gen, before scaling, moist_divg= '
12971+ & ,f16.8)
12972+
12973+ if (divg > -9999.1 .and. divg < -9998.9) then
12974+ ! We have the original, initialized, undefined value of -9999
1295912975 idivg = -99
12976+ else
12977+ idivg = int ((divg * 1e6) + 0.5)
1296012978 endif
1296112979
12962- if (moist_divg > -998.0) then
12963- imoistdivg = int ((moist_divg * 1e6) + 0.5)
12964- else
12980+ if (moist_divg > -9999.1 .and. moist_divg < -9998.9) then
12981+ ! We have the original, initialized, undefined value of -9999
1296512982 imoistdivg = -99
12983+ else
12984+ imoistdivg = int ((moist_divg * 1e6) + 0.5)
1296612985 endif
1296712986
1296812987 if (rh_800_600_smooth > -998.0) then
@@ -12983,10 +13002,12 @@ subroutine output_atcf_gen (outlon,outlat,inp,ist
1298313002 irh_1000_925 = -99
1298413003 endif
1298513004
12986- if (omega500_smooth > -998.0) then
12987- iomega500 = int ((omega500_smooth * 100) + 0.5)
12988- else
13005+ if (omega500_smooth > -9999.1 .and. omega500_smooth < -9998.9)
13006+ & then
13007+ ! We have the original, initialized, undefined value of -9999
1298913008 iomega500 = -99
13009+ else
13010+ iomega500 = int ((omega500_smooth * 100) + 0.5)
1299013011 endif
1299113012
1299213013 if (sst_smooth > -998.0) then
@@ -13134,8 +13155,8 @@ subroutine output_atcf_gen (outlon,outlat,inp,ist
1313413155 & ,'_',a3,', ',i10.10,', 03, ',a4,', ',i3.3,', ',i3,a1
1313513156 & ,', ',i4,a1,', ',i3,', ',i4,', ',a12,4(', ',i4.4)
1313613157 & ,', ',3(i4,', '),3(i6,', '),a1,2(', ',i4),4(', ',i6)
13137- & ,', SHR82, ',i4,', ',i3,3(', ',i4),', ',i9
13138- & ,4 (', ',i4))
13158+ & ,', SHR82, ',i4,', ',i3,3(', ',i4),2( ', ',i9)
13159+ & ,3 (', ',i4))
1313913160
1314013161c bug fix for IBM: flush the output stream so it actually writes
1314113162 flush(66)
@@ -24854,6 +24875,12 @@ subroutine getdata_netcdf (ncfile_id,nc_lsmask_file_id,readflag
2485424875 & ,imax,jmax,nc_zero_ix,f8,igvret)
2485524876 f = f8
2485624877 endif
24878+ if (verb .ge. 3) then
24879+ print *,'After read of separate land_mask file, parm= '
24880+ & ,chparm(ip),' ifh= ',ifh
24881+ & ,' lead time index= ',ltix(ifh),' parm# (ip) = ',ip
24882+ & ,' nc_zero_ix= ',nc_zero_ix,' igvret= ',igvret
24883+ endif
2485724884 else
2485824885 call get_netcdf_real_type (nc_lsmask_file_id,chparm(ip)
2485924886 & ,xtype,ignrret)
@@ -24866,6 +24893,12 @@ subroutine getdata_netcdf (ncfile_id,nc_lsmask_file_id,readflag
2486624893 & ,imax,jmax,ncix,f8,igvret)
2486724894 f = f8
2486824895 endif
24896+ if (verb .ge. 3) then
24897+ print *,'After read of land_mask record from main file,'
24898+ & ,' parm= ',chparm(ip),' ifh= ',ifh
24899+ & ,' lead time index= ',ltix(ifh),' parm# (ip) = ',ip
24900+ & ,' ncix= ',ncix,' igvret= ',igvret
24901+ endif
2486924902 endif
2487024903 else
2487124904 print *,' '
@@ -24884,16 +24917,7 @@ subroutine getdata_netcdf (ncfile_id,nc_lsmask_file_id,readflag
2488424917 & ,imax,jmax,ncix,f8,igvret)
2488524918 f = f8
2488624919 endif
24887- endif
24888-
24889- if (verb .ge. 3) then
24890- print *,' '
24891- if (trkrinfo%read_separate_land_mask_file == 'y') then
24892- print *,'After separate land-sea mask file read, parm= '
24893- & ,chparm(ip),' ifh= ',ifh
24894- & ,' lead time index= ',ltix(ifh),' parm# (ip) = ',ip
24895- & ,' nc_zero_ix= ',nc_zero_ix,' igvret= ',igvret
24896- else
24920+ if (verb .ge. 3) then
2489724921 print *,'After read, parm= ',chparm(ip),' ifh= ',ifh
2489824922 & ,' lead time index= ',ltix(ifh),' parm# (ip) = ',ip
2489924923 & ,' ncix= ',ncix,' igvret= ',igvret
@@ -32949,10 +32973,20 @@ subroutine check_for_closed_wind_circulation (imax,jmax,ip,jp
3294932973 hemisphere = -1.0
3295032974 endif
3295132975
32976+ c print *,' '
32977+ c print *,' ***----------------------------------------------*** '
32978+ c print *,' LLC debug follows'
32979+ c print *,' ***----------------------------------------------*** '
32980+
3295232981 radiusloop: do idist = 1,numdist
3295332982
3295432983 azimuth_ct = 0
3295532984
32985+ c print *,' '
32986+ c print *,'llc1 idist= ',idist,' rdist(idist)= ',rdist(idist)
32987+ c print *,' xcandlon= ',xcandlon,' ycandlat= ',ycandlat
32988+ c print *,' '
32989+
3295632990 azimloop: do iazim = 1,numazim
3295732991
3295832992 bear = ((iazim-1) * 22.5) + 11.25
@@ -32981,12 +33015,22 @@ subroutine check_for_closed_wind_circulation (imax,jmax,ip,jp
3298133015 & ,dx,dy,imax,jmax,trkrinfo,1020,'v',xintrp_v
3298233016 & ,valid_pt,bimct,-99,ibiret2)
3298333017
33018+ c write (6,81) iazim,bear,targlat,targlon,xintrp_u,xintrp_v
33019+ c & ,ibiret1,ibiret2
33020+ c 81 format (1x,'iazim= ',i2,' bear= ',f8.2,' targlat= ',f7.2
33021+ c & ,' targlon= ',f7.2,' xintrp_u= ',f7.2
33022+ c & ,' xintrp_v= ',f7.2,' ibiret1= ',i3
33023+ c & ,' ibiret2= ',i3)
33024+
3298433025 if (ibiret1 == 0 .and. ibiret2 == 0) then
3298533026
3298633027 call getvrvt (xcandlon,ycandlat,targlon,targlat
3298733028 & ,xintrp_u,xintrp_v,vr
3298833029 & ,vt,ifh,igvtret)
3298933030
33031+ c write (6,83) vr,vt
33032+ c 83 format (1x,' vr= ',f7.2,' vt= ',f7.2)
33033+
3299033034 if (bear >= 0. .and. bear < 90.) then
3299133035 iq = 1
3299233036 elseif (bear >= 90. .and. bear < 180.) then
@@ -33001,6 +33045,7 @@ subroutine check_for_closed_wind_circulation (imax,jmax,ip,jp
3300133045 vtct(iq,idist) = vtct(iq,idist) + 1
3300233046
3300333047 if ((hemisphere*vt) >= 8.75) then
33048+ ! For the "free pass" check, use 8.75 m/s (17 kts).
3300433049 ! If cyclonic Vt exceeds 8.75 m/s (17 kts) at this
3300533050 ! azimuth, then increment the counter for this quad by 1.
3300633051 vt_exceed_17kts_ct(iq,idist) =
@@ -33027,12 +33072,18 @@ subroutine check_for_closed_wind_circulation (imax,jmax,ip,jp
3302733072 ! quadrant, but that's okay. What it is *not* able to do here
3302833073 ! is take that 'y' setting away that may have just been set in
3302933074 ! the IF statement above with two azimuths passing 17 kts.
33075+ ! For the check here, we use a slightly lower threshold of
33076+ ! 7 m/s (13.6 kts) than we did above with the free-pass
33077+ ! threshold of 8.75 m/s (17 kts).
3303033078
3303133079 do nq = 1,numquad
3303233080 ! We need at least 2 valid azimuths in order to get a proper
3303333081 ! mean Vt.
3303433082 if (vtct(nq,idist) >= 2) then
3303533083 vtavg = vtsum(nq,idist) / vtct(nq,idist)
33084+ c print *,' +++ nq= ',nq,' vtct(nq,idist)= ',vtct(nq,idist)
33085+ c & ,' vtsum(nq,idist)= ',vtsum(nq,idist)
33086+ c & ,' vtavg= ',vtavg
3303633087 if ((hemisphere*vtavg) >= full_vt_thresh) then
3303733088 ! The mean Vt averaged over the number of azimuths in this
3303833089 ! quadrant (ideally, the max number of azimuths per
@@ -33044,6 +33095,7 @@ subroutine check_for_closed_wind_circulation (imax,jmax,ip,jp
3304433095 quad_pass_half_vt_flag(nq) = 'y'
3304533096 endif
3304633097 else
33098+ c print *,' !! BAD nq= ',nq,' vtct(nq,idist)= ',vtct(nq,idist)
3304733099 vtavg = -9999.0
3304833100 endif
3304933101
@@ -33057,8 +33109,12 @@ subroutine check_for_closed_wind_circulation (imax,jmax,ip,jp
3305733109 endif
3305833110 endif
3305933111
33112+ c print *,' In nq loop, nq= ',nq
33113+ c & ,' vtquadmax(nq)= ',vtquadmax(nq)
33114+
3306033115 enddo
3306133116
33117+
3306233118 enddo radiusloop
3306333119
3306433120 if (verb >= 3) then
0 commit comments