! module calculate_ecf_wx ! ! Coded by Tong Dang, August, 2017 ! This program is to calculate the rest euv flux during eclipse period ! use params_module,only: nlon,nlonp4,nlat ! use shr_kind_mod, only: r8=>shr_kind_r8 implicit none ! contains !----------------------------------------------------------------------- subroutine calculate_ecf(ncol,nlev,lat,lon,year,month,day,uthour, | chi,eclipfactor_EUV, eclipfactor, newmag_factor) use time_manager, only: get_curr_calday, get_curr_date use spmd_utils, only: masterproc use cam_logfile, only : iulog ! ! coded by Dang, 2017 ! use params_module,only: glon,glat ! use hist_module, only: modeltime ! use chapman_module, only: chi ! Args: ! integer,intent(in) :: lat,lon0,lon1,lev0,lev1 ! real,dimension(lev0:lev1,lon0:lon1),intent(out) :: eclipfactor ! | ,newmag_factor integer,intent(in) :: ncol, nlev, year, month, day real(r8),intent(in) :: uthour real(r8),dimension(ncol),intent(in) :: lat, lon, chi real(r8),dimension(ncol,nlev),intent(out) :: eclipfactor_EUV real(r8),dimension(ncol,nlev),intent(out) :: eclipfactor | ,newmag_factor ! Local: real(r8) :: longitude,latitude,height,ut real(r8) :: newmag,EclipseFactor,mag real(r8) :: t1_SLT,t4_SLT,tm_SLT,SLT real(r8) :: calday ! current calendar day integer :: k,i integer :: year_eclip, month_eclip, day_eclip, todsec, ieclipse ! year_eclip=2021 ! month_eclip=12 ! day_eclip=4 year_eclip=2006 month_eclip=03 day_eclip=29 eclipfactor(1:ncol,1:nlev) = 1.0_r8 eclipfactor_EUV(1:ncol,1:nlev) = 1.0_r8 newmag_factor(1:ncol,1:nlev) = 0.0_r8 ! ! Check if this is the eclipse day ! if (year /= year_eclip .or. month /= month_eclip .or. | day /= day_eclip) return ! Currently height is fixed at 300km for all altitudes height=300. ! ut=modeltime(2)+modeltime(3)/60.+modeltime(4)/3600. ! ! if(modeltime(1) .gt. 233 .or. ut .gt. 21.2 .or. ut .lt. 15.6 .or. ! . modeltime(1) .lt. 233) then ! eclipfactor(lev0:lev1,lon0:lon1)=1.0 ! newmag_factor(lev0:lev1,lon0:lon1)=0.0 ! else calday = get_curr_calday() ! write(iulog,*) 'In calculate_ecf 180.-asind(6371.0/(6371.0+height))', ! | 180.-asind(6371.0/(6371.0+height)) ieclipse = 0 do i=1,ncol ! call get_curr_date(year, month, day, todsec) ! uthour = todsec / 60.0_r8 / 60.0_r8 ! if(i .ne. 1 .and. i .ne. 2 .and. i .ne. nlon+3 .and. ! . i .ne. nlon+4) then ! longitude=glon(i-2) ! latitude=glat(lat) ! SLT=modeltime(2)+modeltime(3)/60.+modeltime(4)/3600. ! . + longitude/15. longitude=lon(i) ! write(iulog,*) 'icol,longitude before conversion ', i,longitude if (longitude > 180.0_r8) longitude = longitude - 360.0_r8 ! write(iulog,*) 'icol,longitude after conversion ', i,longitude latitude=lat(i) SLT = uthour + longitude/15. if(SLT .LT. 0) then SLT = SLT + 24. endif if(SLT .GT. 24) then SLT = SLT - 24. endif ! write(iulog,*) 'In calculate_ecf i, uthour,slt,chi(i),longitude,latitude', ! | i, uthour,slt,chi(i),longitude,latitude ! if (latitude < -36._r8 .and. latitude > -37._r8 .and. ! | longitude < -72._r8 .and. longitude > -73._r8) ! | write(iulog,*) 'In calculate_ecf i, uthour,slt,chi(i),longitude,latitude', ! | i, uthour,slt,chi(i)/3.14159_r8*180._r8,longitude,latitude if(chi(i)/3.14159*180. .le. (180.-asind(6371.0/( . 6371.0+height)))) then ! if (masterproc) write(iulog,*) 'Calling eclip ' ! if (masterproc) write(iulog,*) 'Before eclip uthour', uthour call eclip(uthour,SLT,longitude,latitude,height,newmag, . year,month,day,t1_SLT,t4_SLT,tm_SLT,mag) ! if (masterproc) write(iulog,*) 'Done with eclip ' ! if (masterproc) write(iulog,*) 'newmag ', newmag call Eclipse(newmag,EclipseFactor) eclipfactor_EUV(i,1:nlev)=EclipseFactor newmag_factor(i,1:nlev)=newmag if (EclipseFactor < 1.0_r8) | eclipfactor(i,1:nlev)=EclipseFactor / 1.125_r8 ! write(iulog,*) 'In calculate_ecf1 i,longitude,latitude, ! | eclipfactor_EUV(i,1:nlev), ! | eclipfactor(i,1:nlev)', ! | i,longitude,latitude, ! | eclipfactor_EUV(i,1:nlev), ! | eclipfactor(i,1:nlev) ! if (latitude < -36._r8 .and. latitude > -37._r8 .and. ! | longitude < -72._r8 .and. longitude > -73._r8) ! | write(iulog,*) 'In calculate_ecf1 i, uthour,slt,chi(i),longitude,latitude,newmag,EclipseFactor', ! | i, uthour,slt,chi(i)/3.14159_r8*180._r8,longitude,latitude,newmag,EclipseFactor,eclipfactor(i,1) ! eclipfactor(1:ncol,1:nlev)=1.0 ! newmag_factor(1:ncol,1:nlev)=0.0 else eclipfactor(i,1:nlev)=1.0_r8 eclipfactor_EUV(i,1:nlev)=1.0_r8 newmag_factor(i,1:nlev)=0.0_r8 endif ! if (latitude < -36._r8 .and. latitude > -37._r8 .and. ! | longitude < -72._r8 .and. longitude > -73._r8) ! | write(iulog,*) 'In calculate_ecf2 i, uthour,slt,chi(i),longitude,latitude,newmag,EclipseFactor', ! | i, uthour,slt,chi(i)/3.14159_r8*180._r8,longitude,latitude,newmag,EclipseFactor,eclipfactor(i,1) ! if (latitude < -36._r8 .and. latitude > -37._r8 .and. ! | longitude < -72._r8 .and. longitude > -73._r8) ! | ieclipse = i ! if (masterproc .and. eclipfactor(i,1) < 1.0_r8) write(iulog,*) ! if (eclipfactor(i,1) < 1.0_r8) write(iulog,*) ! |'In calculate_ecf after Eclipse,i,uthour,SLT,longitude,latitude,eclipfactor(i,1), chi(i) ', ! | i,uthour,SLT,longitude,latitude,eclipfactor(i,1),chi(i) ! if (masterproc) write(iulog,*) 'Done with Eclipse ' ! if (masterproc) write(iulog,*) 'MINVAL(eclipfactor) ', ! . MINVAL(eclipfactor) ! if (masterproc) write(iulog,*) 'MAXVAL(eclipfactor) ', ! . MAXVAL(eclipfactor) ! if (masterproc .and. ANY(eclipfactor < 0.9) .and. ANY(eclipfactor > 0.1)) ! . write(iulog,*) 'eclipfactorall ', eclipfactor ! endif ! if (ieclipse > 0) ! | write(iulog,*) 'In calculate_ecf end of column loop, ieclipse, i, eclipfactor(ieclipse,1) ', ! | ieclipse, i, eclipfactor(ieclipse,1) enddo ! if(lon0 .eq. 1) then ! eclipfactor(lev0:lev1,1)=eclipfactor(lev0:lev1,3) ! eclipfactor(lev0:lev1,2)=eclipfactor(lev0:lev1,3) ! endif ! if(lon1 .eq. nlon+4)then ! eclipfactor(lev0:lev1,nlon+3)=eclipfactor(lev0:lev1,nlon+2) ! eclipfactor(lev0:lev1,nlon+4)=eclipfactor(lev0:lev1,nlon+2) ! endif ! ! endif ! write(iulog,*) 'End of calculate_ecf, ieclipse, eclipfactor(ieclipse,1) ', ! | ieclipse, eclipfactor(ieclipse,1) end subroutine calculate_ecf !----------------------------------------------------------------------- Subroutine Eclipse(ecf,EclipseFactor) CC calculation of the ratio of φTi / φT0 = (12.3Spi+Sci)/(12.3sp0 + Sc0) CC calculation of the unmasked area for photosphere for P and C regions CC during the solar eclipse. CC P region is an inner central disk over the "photosphere"; CC C region is an outer ring between the "coronal" and "phosphere". CC reference to J. J Curto et al.(2006), Modeling the geomagnetic effects CC caused by the solar eclipse of 11 August 1999. CC initialize implicit real*8 (A-Z) Rp=1.0; Rc=2.1; pi=3.1415926; Sp0=pi*Rp*Rp; Sc0=pi*(Rc*Rc-Rp*Rp); CC dl should be smaller than Rp+Rc; that is dl <= 3.1 C v=3.1*2/Duration !! for duration=240 mins v= 3.1/120 mins =0.025833 C !! for duration=180 mins v= 3.1/90 mins =0.034444 C dl=dabs(3.1-v*eclipseTime); dl=2.0-ecf*2; if(dl<=Rc-Rp) then A1=0; E=2*Rp*Rp*acos(dl/2/Rp) - dl*sqrt(Rp*Rp - dl*dl*0.25); end if if(dl>Rc-Rp .AND. dl<2*Rp) then E=2*Rp*Rp*acos(dl/2/Rp) - dl*sqrt(Rp*Rp - dl*dl*0.25); Xd=(Rc*Rc - Rp*Rp + dl*dl)/2/dl; A2=Rc*Rc*acos(Xd/Rc) - Xd*sqrt(Rc*Rc -Xd*Xd); A3=Rp*Rp*acos((dl-Xd)/Rp) - (dl-Xd)* & sqrt(Rp*Rp-(dl-Xd)*(dl-Xd)); A1=pi*Rp*Rp - A2 -A3; endif if(dl>=2*Rp .AND. dl<=3.1) then E=0; Xd=(Rc*Rc - Rp*Rp + dl*dl)/2/dl; A2=Rc*Rc*acos(Xd/Rc) - Xd*sqrt(Rc*Rc -Xd*Xd); A3=Rp*Rp*acos((dl-Xd)/Rp) - (dl-Xd)* & sqrt(Rp*Rp-(dl-Xd)*(dl-Xd)); A1=pi*Rp*Rp - A2 -A3; endif Sc=pi*(Rc*Rc - 2*Rp*Rp) + E + A1; Sp=pi*Rp*Rp - E; EclipseFactor=(12.3*Sp+Sc)/(12.3*Sp0+Sc0); if(ecf==0) EclipseFactor=1.0 End Subroutine Eclipse !----------------------------------------------------------------------- end module calculate_ecf_wx