9a10
>     !YuxiangHe+
18c19,20
<    
---
> 
>     LOGICAL :: file_exist   
20,21c22
<     INTEGER           :: YEAR,MONTH,DAY,DOY
<     CHARACTER(LEN=20) :: Tbase
---
>     INTEGER           :: YEAR,MONTH,DAY
29c30
<     iNTEGER,PARAMETER :: i_era5=1440,j_era5=721
---
>     iNTEGER,PARAMETER :: i_ech=1440,j_ech=721
31,32d31
<     LOGICAL :: file_exist1,file_exist2,file_exist3,file_exist4
<     LOGICAL :: file_exist5,file_exist6,file_exist7,file_exist8  
45a45
> 
51a52
>     CHARACTER(LEN=2)   :: stype
64,66c65,67
<     ! all levels
<     integer(kind = 4)  :: dim_time,dim_lev,dim_lat,dim_lon
<     REAL(dp),allocatable,dimension(:)       :: pp,zz,tt,qq,nn 
---
>     ! all levels and surface
>     integer(kind = 4)  :: dim_lev
>     REAL(dp),allocatable,dimension(:)       :: pp,zz,tt,qq,nn !pp is constant
68d68
<     REAL(dp),allocatable,dimension(:,:,:,:) :: Z2,T2,Q2,P2
72,73c72,73
<     CHARACTER(LEN=512) :: era5_dir
<     CHARACTER(LEN=512) :: era5_name(8)
---
>     CHARACTER(LEN=512) :: ECH_dir
>     CHARACTER(LEN=512) :: ECH_file(9)
75d74
<     CHARACTER(LEN=512) :: RO_dir
77,78c76,77
<     LOGICAL            :: cross_day,curr_day_read_flag &
<                                    ,next_day_read_flag
---
>     CHARACTER(LEN=512) :: RO_dir
>     LOGICAL            :: read_flag 
82,83c81,82
<     REAL(DP)  :: sJDF,   eJDF,   cJDF,   tJDF
<     INTEGER   :: sCDT(8),eCDT(8),cCDT(8),RO_CDT(8)
---
>     REAL(DP)  :: sJDF,   eJDF,   cJDF
>     INTEGER   :: sCDT(8),eCDT(8),cCDT(8),RO_CDT(8),cCDT_int(8)
97a97
>      
99,100c99
<     !! ERA5 3D Input Path: input data 
<     ! era5_dir    = "./data/ERA5/ERA5_hires"
---
>     if(iargc() .ne. 7)then
102d100
<     if(iargc() .ne. 6)then
107c105
<       write(*,*) "    comand Year Month Day EAR5_Dir COSMIC2_List Output_dir"
---
>       write(*,*) "    comand Year Month Day ECH_Dir COSMIC2_List Output_dir stype"
109,116c107,110
<       write(*,*) "    EAR5_Dir includes next day data (/):"
<       write(*,*) "      2019/e5.oper.an.pl.128_129_z.ll025sc.2019100100_2019100123.nc"  
<       write(*,*) "      2019/e5.oper.an.pl.128_130_t.ll025sc.2019100100_2019100123.nc" 
<       write(*,*) "      2019/e5.oper.an.pl.128_133_q.ll025sc.2019100100_2019100123.nc" 
<       write(*,*) "      2019/e5.oper.an.pl.128_129_z.ll025sc.2019100200_2019100223.nc" 
<       write(*,*) "      2019/e5.oper.an.pl.128_130_t.ll025sc.2019100200_2019100223.nc" 
<       write(*,*) "      2019/e5.oper.an.pl.128_133_q.ll025sc.2019100200_2019100223.nc" 
<       write(*,*) "    COSMIC2_List:"
---
>       write(*,*) "    ECH_Dir includes next day data (/):"
>       write(*,*) "      2026/01/UAD01070000010709001_ml.nc "
>       write(*,*) "      2026/01/UAD01071200010721001_ml.nc "
>       write(*,*) "    COSMIC2_List(/):"
118a113,115
>       write(*,*) "      ./output/" 
>       write(*,*) "    stype(/):" 
>       write(*,*) "      03; 09; "
120,125c117,119
<       write(*,*) "    1)e.g.: ./star_fm_ERA5_main_v2.exe 2019 10 1 -1 "  &
<                        //"/data/smcd5/xzhou/data_ERA/ERA5/ " &
<                        //"./output/list_2019.274.txt ./output/"
<       write(*,*) "    2)e.g.: ./star_fm_ERA5_main_v2.exe 2021 1 1 1 "  &
<                        //"/data/smcd5/xzhou/data_ERA/ERA5_cds/ " &
<                        //"./output/list_2021.001.txt ./output/"
---
>       write(*,*) "    e.g.: ./star_fm_ECH_main_v2.exe 2026 1 7 "  &
>                        //"/data2/xinjiaz/data_ECMWF/ " &
>                        //"./output/list_2026.006.txt ./output/ 09"
138a133
> 
140c135
<      era5_dir    = arg
---
>      ech_dir    = arg
144a140,141
>     call getarg(7,arg)
>      stype  = arg
152a150
> 
154d151
<     !print*,YEAR," ",MONTH," ",DAY
158a156
>     !eCDT = (/YEAR,MONTH,DAY+1,0,0,0,0,0/)
164c162
< ! Step2: daily loop 
---
> ! Step1: daily loop 
168,176d165
< 
<         WRITE ( Tbase, "(A2,I4.4)" ) "JD", cCDT(1)
<         CALL TimeSince ( cCDT, tJDF, 1, Base=Tbase )
< 
<         ! Include end date in calculation (1 day is added)
<         ! allow for JDFs from midday
<         DOY = INT ( tJDF + HalfDay ) + 1 
< 
<         !print*,DOY
179,244c168
<         !! DoY to CDT
<         curr_day_read_flag=.false.
<         next_day_read_flag=.false.
< 
<         !! get era5 file name 
<         !! module_read_write.f90
<         CALL get_era5hres_name(era5_dir,era5_name,cCDT)
<         print*,'era5-hire, file1:',trim(era5_name(1))
<         print*,'era5-hire, file2:',trim(era5_name(2))
<         print*,'era5-hire, file3:',trim(era5_name(3))
<         print*,'era5-hire, file4:',trim(era5_name(4))
<         print*,'era5-hire, file5:',trim(era5_name(5))
<         print*,'era5-hire, file6:',trim(era5_name(6))
<         print*,'era5-hire, file7:',trim(era5_name(7))
<         print*,'era5-hire, file8:',trim(era5_name(8))
< 
<         ! z,t,q,p
<         inquire( file=trim(era5_name(1)), exist=file_exist1 )  ! Works with gfortran, but not ifort
<         inquire( file=trim(era5_name(2)), exist=file_exist2 )  ! Works with gfortran, but not ifort
<         inquire( file=trim(era5_name(3)), exist=file_exist3 )  ! Works with gfortran, but not ifort
<         inquire( file=trim(era5_name(4)), exist=file_exist4 )  ! Works with gfortran, but not ifort
< 
<         inquire( file=trim(era5_name(5)), exist=file_exist5 )  ! Works with gfortran, but not ifort
<         inquire( file=trim(era5_name(6)), exist=file_exist6 )  ! Works with gfortran, but not ifort
<         inquire( file=trim(era5_name(7)), exist=file_exist7 )  ! Works with gfortran, but not ifort
<         inquire( file=trim(era5_name(8)), exist=file_exist8 )  ! Works with gfortran, but not ifort
< 
<         if(.not. file_exist1) then
<           write(*,*) trim(era5_name(1))," doesn't exist!" 
<           return
<         endif
< 
<         if(.not. file_exist2) then
<           write(*,*) trim(era5_name(2))," doesn't exist!" 
<           return
<         endif
< 
<         if(.not. file_exist3) then
<           write(*,*) trim(era5_name(3))," doesn't exist!" 
<           return
<         endif
< 
<         if(.not. file_exist4) then
<           write(*,*) trim(era5_name(4))," doesn't exist!" 
<           return
<         endif
< 
<         if(.not. file_exist5) then
<           write(*,*) trim(era5_name(5))," doesn't exist!" 
<           return
<         endif
< 
<         if(.not. file_exist6) then
<           write(*,*) trim(era5_name(6))," doesn't exist!" 
<           return
<         endif
< 
<         if(.not. file_exist7) then
<           write(*,*) trim(era5_name(7))," doesn't exist!" 
<           return
<         endif
< 
<         if(.not. file_exist8) then
<           write(*,*) trim(era5_name(8))," doesn't exist!" 
<           return
<         endif
---
>         cCDT_int = (/cCDT(1),cCDT(2),cCDT(3),0,0,0,0,0/)
245a170
>         read_flag=.false.
246a172
>         !! get ECH file name 
248,252c174,197
<         CALL get_era5hres_pztq_dim_memory(era5_name(2)& ! 2 or 3 is good
<           ,dim_time,dim_lev,dim_lat,dim_lon &
<           ,P1,Z1,T1,Q1 &
<           ,P2,Z2,T2,Q2 &
<           ,pp,zz,tt,qq,nn)
---
>         CALL get_ech_name(ech_dir,ECH_file,cCDT,stype)
> 
>         print*,'ECH(',stype,'), file1 (yesterday):',trim(ECH_file(1))
>         print*,'ECH(',stype,'), file2 (yesterday):',trim(ECH_file(2))
>         print*,'ECH(',stype,'), file3 (today 00h):',trim(ECH_file(3))
>         print*,'ECH(',stype,'), file4 (today 00h):',trim(ECH_file(4))
>         print*,'ECH(',stype,'), file5 (today 12h):',trim(ECH_file(5))
>         print*,'ECH(',stype,'), file6 (today 12h):',trim(ECH_file(6))
>         print*,'ECH(',stype,'), file7 (tomorrow) :',trim(ECH_file(7))
>         print*,'ECH(',stype,'), file8 (tomorrow) :',trim(ECH_file(8))
>         print*,'ECH(',stype,'), file9 (z from ERA):',trim(ECH_file(9))
> 
>         do i = 1, 9
>           inquire( file=trim(ECH_file(i)), exist=file_exist )
>           if (.not. file_exist) then
>             write(*,*) trim(ECH_file(i)), " doesn't exist!"
>             return
>           endif
>         enddo
> 
>         CALL get_ech_pztq_dim_memory(ECH_file(1) &
>           ,dim_lev &
>           ,P1,Z1,T1,Q1 &   ! raw data, in 4-dim
>           ,pp,zz,tt,qq,nn) ! interploated to RO, in 1-dim
284c229
<         
---
> 
288c233
<             if(cCDT(1) .eq. RO_CDT(1) .and. cCDT(2) .eq. RO_CDT(2))then
---
>             if(cCDT(1) .eq. RO_CDT(1) .and. cCDT(2) .eq. RO_CDT(2)) then
295c240,241
<             !print*,'before get_era5_index'
---
>             print*,"Time Found: ", trim(fileStamp1)
> 
297,298c243
<             cross_day = .false.
<             CALL get_era5_index(RO_CDT,cross_day &
---
>             CALL get_ech_index(RO_CDT,cCDT_int &
300c245
<               ,lonlat_coeff,time_coeff,lon_in,lat_in)
---
>               ,lonlat_coeff,time_coeff,lon_in,lat_in,stype)
302,306c247,250
<             !!! read NetCDF file to get zz(:),tt(:),qq(:),pp(:)
<             ! --- height from top to bottom in ERA5 file
<             ! module_read_write.f90
<             CALL read_era5hres_upLevels_ztq(era5_name,cross_day &
<                 ,zz,tt,qq,pp &    ! intepolated to the RO location
---
>             !!! read NetCDF file to get zz(:),tt(:),qq(:)
>             ! --- height from top to bottom in GFS file
>             CALL read_ech_3d(ECH_file&
>                 ,pp,zz,tt,qq &
309,311c253,254
<                 ,Z1,T1,Q1,P1,Z2,T2,Q2,P2 &   ! raw data
<                 ,curr_day_read_flag,next_day_read_flag &
<                 ,dim_lev)
---
>                 ,P1,Z1,T1,Q1 &
>                 ,read_flag, dim_lev)
314c257
<             call calculate_refractivity(tt,pp/100.0D0,nn,q=qq) ! K,hPa,kg/kg
---
>             call calculate_refractivity(tt,pp/100.0D0,nn,rh=qq) ! K,hPa,%(kg/kg)
316,321c259,263
<             !print*,curr_day_read_flag,next_day_read_flag
<             !print*,'zz',zz
<             !print*,'tt',tt
<             !print*,'qq',qq
<             !print*,'pp',pp
<             !print*,'nn',nn
---
>             !print*,'pp:',pp
>             !print*,'zz:',zz
>             !print*,'tt:',tt
>             !print*,'qq:',qq
>             !print*,'nn:',nn
323d264
<             ! ---------------------------
333c274
<   
---
> 
336,337c277,278
<             ! zz, tt, pp, nn, model_z are input.
<             ! model_t, model_p, model_n are output
---
>             !!! star interpolation procedures
>             !pp: Pa, model_p: Pa, model_h:meter
339c280
<                        zz,model_z &  ! geopotential
---
>                        zz,model_z &
342c283
<                      , nn,model_n)
---
>                      , nn,model_n )
344,348c285,289
<             !print*,'model_h',model_h
<             !print*,'model_z',model_z
<             !print*,'model_t',model_t
<             !print*,'model_p',model_p
<             !print*,'model_n',model_n
---
>             !print*,'model_h:',model_h
>             !print*,'model_z:',model_z
>             !print*,'model_t:',model_t
>             !print*,'model_p:',model_p
>             !print*,'model_n:',model_n
362,363d302
<             !!! doesn't work for positive undulation, bending angle output are NaN
<             !model_impact(:) = (model_n(:)*(1d-6)+1.d0)*(model_h(:)+Rearth+undulation/1000.0)
371c310
<           
---
>             
373c312
< 
---
>             
382a322,323
>             !print*,'model_b',model_b
> 
386c327,329
<             fname_alp  = trim(output_dir )//"era5hres_bend_"//trim(fileStamp1)// ".nc"
---
>                                            ! ech09_bend_GN05.2026.019.16.42.E25.nc
>             fname_alp  = trim(output_dir )//"ech"//stype//"_bend_"//trim(fileStamp1)// ".nc"
> 
394c337,338
<         enddo 
---
> 
>         enddo
