/home/tgcm/bin/Diffsrc executing from /home/foster/tiegcm/amie/modsrc_oct03 at Thu Apr 27 11:17:46 MDT 2006 Using tgcmroot from TGCMROOT env var: /home/tgcm Source code directory = /home/tgcm/tiegcm1 ======================================================================== Diff of /home/tgcm/tiegcm1/advance.F and /home/foster/tiegcm/amie/modsrc_oct03/advance.F 2,5d1 < ! Include the option of setting the lower boundary by using GSWM output < ! this affects subroutine addiag < ! < ! 11a8,9 > use heelis_module,only: heelis,colath > use weimer_module,only: weimer01 13c11,13 < | f107,f107a,ctpoten,power,calendar_advance --- > | f107,f107a,ctpoten,power,calendar_advance,ntimes_ctpoten, > | ctpoten_time,ntimes_power,power_time,potential_model,iamie, > | iamie_bkg,ilow_proton 16c16 < use cons_module,only: dt,re,cs,racs,kut --- > use cons_module,only: dt,re,cs,racs,kut,ylatm,rtd 21c21 < use aurora_module,only: aurora_cons --- > use aurora_module,only: aurora_cons,cusp2d,drzl2d,alfa2d,nflx2d 24,25c24,27 < use dynamo_module,only: prep_dynamo,dynamo < use dynamo_old_module,only: prep_dynamo_old,dynamo_old --- > use dynamo_module,only: prep_dynamo,prep_dynamo_dyn0, dynamo, > | nodynamo,phihm,nmlat0,nmlat > use dynamo_old_module,only: prep_dynamo_old,dynamo_old, > | phihm_old=>phihm 30c32,34 < use params_module,only: nlon,nlat,nmlev --- > use params_module,only: nlonp4,nlat > use amie_module,only: getamie,tiepot > use low_proton_module,only: getproton 45a50 > | rtcmp_periodic_f2d , 174a180,187 > ! Interpolate power and/or ctpoten to current model time, if time-dependent > ! values were read from input: > if (ntimes_ctpoten > 0) > | call set_index(ctpoten_time,ntimes_ctpoten,nsecs,ctpoten, > | 'ctpoten') > if (ntimes_power > 0) > | call set_index(power_time,ntimes_power,nsecs,power,'power') > ! 214a228,243 > ! Get AMIE data > if (iamie==1) then > iprint = 0 > if (istep==1) iprint = 1 > if (iprint>0) write(6,"('advance calling getamie...')") > call getamie(iyear,iday,int(secs),iamie_bkg,iprint) > endif > ! > ! Get auroral proton data > if (ilow_proton==1) then > iprint = 0 > if (istep==1) iprint = 1 > if (iprint>0) write(6,"('advance calling getproton...')") > call getproton(iyear,iday,int(secs),iprint) > endif > ! 225,227d253 < if (igetgpi > 0) write(6,"('gpi: f107=',f7.2,' f107a=',f7.2, < | ' ctpoten=',f7.2,' power=',f7.2)") < | f107,f107a,ctpoten,power 232,234d257 < if (igetgpi > 0) write(6,"('gpi: f107=',f7.2,' f107a=',f7.2, < | ' ctpoten=',f7.2,' power=',f7.2)") < | f107,f107a,ctpoten,power 235a259,263 > if (igetgpi > 0) write(6,"('gpi: f107=',f7.2,' f107a=',f7.2, > | ' ctpoten=',f7.2,' power=',f7.2)") > | f107,f107a,ctpoten,power > if (ntimes_ctpoten > 0.or.ntimes_power > 0) > | write(6,"('ctpoten=',f8.2,' power=',f8.2)") ctpoten,power 388a417,449 > ! Dynamo calls Heelis (heelis.F), Weimer (wei01gcm.F), or neither > ! for high latitude electric potential, according to user-provided > ! "model_potential". > ! Get high latitude (Heelis or other) colatitudes, NH pfrac, and poten phihm. > ! If Weimer is used, then theta0,phid etc is changed before use in aurora > ! in dynamics. > ! If an amie run, ignore potential_model, and use electric potential from > ! amie module. > ! > if (iamie <= 0) then > if (potential_model == 'WEIMER') then > call weimer01 > elseif (potential_model == 'HEELIS') then > call heelis > elseif (potential_model == 'NONE') then > do j=1,nmlat0 > do i=1,nmlonp1 > phihm(i,j) = 0. > phihm_old(i,j) = 0. > enddo ! i=1,nmlonp1 > enddo ! j=1,nmlat0 > call colath > endif > ! > ! If an amie run, use tiepot from amie.F: > ! (phihm and tiepot are dimensioned (nmlonp1,nmlat)) > else > phihm = tiepot ! whole array op > ! call fminmax(phihm,nmlonp1*nmlat,fmin,fmax) > ! write(6,"('advance: amie phihm min,max=',2e12.4)") fmin,fmax > call colath > endif > ! 400a462,472 > ! Save Cusp, drizzle, alfa and nflux to secondary history to plot later > call addfsech_ij ('CUSP','2D CUSP','0-1',cusp2d,1,nlonp4, > | 1,nlat) > call addfsech_ij ('DRIZZLE','2D DRIZZLE','0-1',drzl2d,1, > | nlonp4,1,nlat) > call addfsech_ij ('ALFA','2D ALFA','keV',alfa2d,1, > | nlonp4,1,nlat) > call addfsech_ij ('NFLUX','2D NFLUX','#/cm3?',nflx2d,1, > | nlonp4,1,nlat) > > ! 403c475 < if (idynamo > 0) then --- > if (idynamo >= 0) then 418a491,495 > else if (idynamo==0) then > call prep_dynamo_dyn0( > | z (levd0,lond0,latd0,itp), > | 1,nlevp1,lon0,lon1,lat0,lat1) > if (debug) write(6,"('advance after prep_dynamo_dyn0')") 441a519,524 > ! > ! Call only part of the new dynamo to get the high-lat potential into E > ! fields if idynamo = 0 (high lat potential can be Heelis, Weimer, or None) > else if (idynamo==0) then ! no dynamo > call nodynamo > if (debug) write(6,"('advance after nodynamo')") 501a585,588 > ! Saving extraoutput here from master only - G. Lu, Aug. 27, 2003 > ! call prep_output(1,nlevp1,lon0,lon1,lat0,lat1) > ! if (mytid==0) call extraoutput > ! end of saving extraoutput 597a685,686 > write(6,"('Total rtc secs in mp_periodic_f2d = ',f8.2)") > | rtcmp_periodic_f2d 665a755,831 > !----------------------------------------------------------------------- > subroutine set_index(rindex,ntimes,msecs,outindex,name) > ! > ! User has provided time-dependent ctpoten (ctpoten_time) and/or > ! power (power_time) via namelist input. This routine interpolates > ! these inputs to the current model time msecs, returning outindex. > ! This routine is called separately (from advance) for ctpoten and power, > ! (i.e., rindex will be either ctpoten_time or power_time from input). > ! Note times in seconds are 8-byte integers. > ! > use params_module,only: > | mxind_time ! max number of time-dependent solar index points > implicit none > ! > ! Args: > real,intent(in) :: > | rindex(4,mxind_time) ! user input times and values (day,hour,min,value) > integer,intent(in) :: > | ntimes ! number of valid time/values in rindex(:,1:ntimes) > integer(kind=8),intent(in) :: > | msecs ! current model time in seconds > real,intent(out) :: outindex ! output interpolated value > character(len=*),intent(in) :: name > ! > ! Local: > integer :: i > integer(kind=8) :: nsec0,nsec1 > ! > ! External: > integer(kind=8),external :: mtime_to_nsec > real,external :: finterp_bigints > ! > ! If model time is beyond last rindex time, use last rindex time: > nsec1 = mtime_to_nsec(int(rindex(1:3,ntimes))) > if (msecs > nsec1) then > outindex = rindex(4,ntimes) > goto 100 > endif > ! > ! Bracket model time: > do i=1,ntimes-1 > nsec0 = mtime_to_nsec(int(rindex(1:3,i))) > nsec1 = mtime_to_nsec(int(rindex(1:3,i+1))) > ! > ! If model time is at a provided time, interpolation is not necessary: > if (nsec0 == msecs) then > outindex = rindex(4,i) > goto 100 > endif > if (nsec1 == msecs) then > outindex = rindex(4,i+1) > goto 100 > endif > ! > ! Interpolate to model time msecs: > if (msecs >= nsec0 .and. msecs <= nsec1) then > outindex = finterp_bigints(rindex(4,i),rindex(4,i+1),nsec0, > | nsec1,msecs) > goto 100 > endif > enddo ! i=1,ntimes-1 > ! > ! Error if model time could not be bracketed. This should not happen, > ! but you never know... > write(6,"(/,'>>> set_index: could not bracket model time ', > | i10)") msecs > write(6,"('ntimes=',i3)") ntimes > do i=1,ntimes > write(6,"(' i=',i3,' ntimes=',i3,' day,hr,min=',3f7.2, > | ' value=',e12.4)") i,ntimes,rindex(1:3,i),rindex(4,i) > enddo > stop 'set_index' > ! > ! Report to stdout: > 100 continue > ! write(6,"('set_index: ',a,' = ',e12.4)") name,outindex > end subroutine set_index ======================================================================== Diff of /home/tgcm/tiegcm1/advec.F and /home/foster/tiegcm/amie/modsrc_oct03/advec.F 8c8 < | dphi_1div12,re_inv --- > | dphi_1div12,dphi_2div3,re_inv 69a70,84 > ! 2nd order finite difference in latitude (j+1,j-1): > ! > ! (10/21/03 btf: 2nd order diffs (next 2 statements) were missing > ! in tgcm13,14,15, and tiegcm1 until this date.) > ! > ! VBARP = (V(J+1)+V(J-1))/2 > vbarp(k,i) = (vn(k,i,lat+1,itp)+ > | vn(k,i,lat-1,itp))*0.5 > ! > ! D2X = dphi_2div3*(X(J+1)-X(J-1)) > d2x(k,i) = ((f(k,i,lat+1)-f(k,i,lat-1))*dphi_2div3)* > | vbarp(k,i) > ! > ! Fourth order finite difference in latitude (j+2,j-2): > ! 71d85 < ! (comment in tgcm15 apparently wrong: VBARP = (V(J+1)+V(J-1))/2) ======================================================================== >>> WARNING: Cannot find source file /home/tgcm/tiegcm1/amie.F ======================================================================== Diff of /home/tgcm/tiegcm1/aurora.F and /home/foster/tiegcm/amie/modsrc_oct03/aurora.F 2,3d1 < ! Pre-processor macros for grid resolution: < ! 55c53,54 < | pi ! 4.*atan(1.) --- > | pi, ! 4.*atan(1.) > | crit ! may be changed from default if amie run 56a56,58 > use amie_module,only: crad,phida,ekvg,efxg,hpi_sh_amie, > | hpi_nh_amie,pcp_sh_amie,pcp_nh_amie > use input_module,only: iamie 77c79,80 < | ed = 0.5 --- > | ed = 0.5, > | emin = 0.50 169a173,175 > ! TEMP > real :: cusp2d(nlonp4,nlat),drzl2d(nlonp4,nlat), > | alfa2d(nlonp4,nlat),nflx2d(nlonp4,nlat) 197c203 < real :: arad(2) --- > real :: arad(2),hp_amie,c25,c35 333a340,361 > ! Recalculate some parameters if an amie run: > if (iamie > 0) then > theta0(:) = crad(:) ! radians from amie module > phid(:) = phida(:) ! radians from amie module > hp_amie = max(hpi_sh_amie,hpi_nh_amie) > plevel = 2.09*alog(hp_amie) > h1 = min(2.35,0.83+0.33*plevel) > h2 = 2.87 + 0.15 * plevel > h0 = 0.5 * (h1 + h2) * dtr > roth = (12.18- 0.89*plevel)*15. > rote = (2.62 - 0.55*plevel)*15. > rroth = roth * dtr > rrote = rote * dtr > ! > ! Set critical colatitude crit(2) 40 deg -- G. Lu 5/11/98 > ! Therefore, crit2 = 40, crit1 = 25-30, depending on theta0 > ! > c35 = min(30.,(theta0(1)+theta0(2))*0.5 + 5.0) > c25 = max(25.,c35) > crit(1) = c25 * dtr > endif > ! 336c364 < write(6,"(/,'aurora_cons:')") --- > write(6,"(/,'aurora_cons (iamie=',i2,'):')") iamie 353c381 < subroutine aurora(tn,o2,o1,barm,lev0,lev1,lon0,lon1,lat) --- > subroutine aurora(z,tn,o2,o1,barm,lev0,lev1,lon0,lon1,lat) 374a403 > | z, ! geopotential 449a479,487 > > ! TEMP > ! Save aurora oval on exit to full arrays so can save later to secondary history > do i=lon0,lon1 > cusp2d(i,lat)=cusp(i) > drzl2d(i,lat)=drizl(i) > alfa2d(i,lat)=alfa(i) > nflx2d(i,lat)=flux(i) > enddo 462c500 < | tn,o2,o1,barm,zsb,dz,lev0,lev1,lon0,lon1,lat) --- > | z,tn,o2,o1,barm,zsb,dz,lev0,lev1,lon0,lon1,lat) 475a514,515 > ! (if this is an amie run, then theta0 and phid have been taken > ! from amie, see aurora_cons) 485c525 < | pi_cusp = 3.1415927 ! hard-wired by tgcm15 --- > | pi_cusp = 3.1415927 ! hard-wired by tgcm15 517c557 < ! write(6,"('aurora_cusp: lat=',i2,' cusp=',/,(6e12.4))") --- > ! write(6,"('aurora_cusp: lat=',i3,' cusp=',/,(6e12.4))") 535a576,577 > real,parameter :: s10=0.174532925 > real :: clat 594a637,656 > ! > ! Recalculate alfa, flux, drizl if an amie run (as in old sub amiepa): > if (iamie > 0) then > ! > ! Insure ekvg mean energy >= 1.: > ekvg(i,lat) = max(ekvg(i,lat),1.) > alfa(i+2) = ekvg(i,lat)/2. > alfa2(i+2) = alfa20 > clat = acos(sin(abs(dlat_aur(i)))) > drizl(i+2) = exp(-((clat-crad(ihem)+ > | abs(clat-crad(ihem)))/(2.*s10))**2) > ! drizl(i+2) = exp(-((colat(i)-crad(ihem)+ > ! | abs(colat(i)-crad(ihem)))/s10)**2) > flux(i+2) = max(efxg(i,lat)/(2.*alfa(i+2)*1.602e-9),1.e-20) > flux2(i+2) = e20*(1.-re2*coslamda(i))* > | exp(-((clat-crad(ihem))/ > ! | exp(-((colat(i)-crad(ihem))/ > | halfwidth(i))**2) / (2.*alfa2(i+2)*1.602E-9) > > endif 614a677,685 > > ! write(6,"('aurora_heat: lat=',i3,' alfa=',/,(6e12.4))") lat,alfa > ! write(6,"('aurora_heat: lat=',i3,' flux=',/,(6e12.4))") lat,flux > ! write(6,"('aurora_heat: lat=',i3,' drizl=',/,(6e12.4))") lat,drizl > ! write(6,"('aurora_heat: lat=',i3,' alfa2=',/,(6e12.4))") lat,alfa2 > ! write(6,"('aurora_heat: lat=',i3,' alfa3=',/,(6e12.4))") lat,alfa3 > ! write(6,"('aurora_heat: lat=',i3,' flux2=',/,(6e12.4))") lat,flux2 > ! write(6,"('aurora_heat: lat=',i3,' flux3=',/,(6e12.4))") lat,flux3 > 619c690 < | flux3,tn,o2,o1,barm,zpbot,dzp,lev0,lev1,lon0,lon1,lat) --- > | flux3,z,tn,o2,o1,barm,zpbot,dzp,lev0,lev1,lon0,lon1,lat) 623a695 > use input_module,only: ilow_proton 653c725 < | tn,o2,o1,barm --- > | z,tn,o2,o1,barm 668c740,741 < real :: qia(5) ! low energy proton source (not in use, 1/02) --- > real,dimension(lev0:lev1,lon0:lon1,5) :: > | qia ! low energy proton source (not in use, 1/02) 671c744 < qia(:) = 0. --- > qia(:,:,:) = 0. 727a801,809 > ! > ! Contribution of low energy protons (From Marina Galand's original code): > ! qia(k,i,5) is output ion production for N2+, O2+, O+, NO+, N+ > ! (top level zkmxp not defined) > ! Modified by G. Lu - Oct. 13, 2003 > ! > if (ilow_proton > 0) > | call low_proton(z,tn,o1,o2,barm,qia,lev0,lev1,lon0,lon1,lat) > ! 771c853,854 < qo2p_aur(k,i) = qsum(k,i)*o2(k,i)/(rmass_o2*denom(k,i))+qia(2) --- > qo2p_aur(k,i) = qsum(k,i)*o2(k,i)/(rmass_o2*denom(k,i))+ > | qia(k,i,2) 776c859,860 < | 0.56*o1(k,i)*rmassinv_o1)/denom(k,i) + qia(3) --- > | 0.56*o1(k,i)*rmassinv_o1)/denom(k,i) + > | qia(k,i,3) 781c865,866 < qn2p_aur(k,i) = qsum(k,i)*0.7*xn2/(rmass_n2*denom(k,i))+qia(1) --- > qn2p_aur(k,i) = qsum(k,i)*0.7*xn2/(rmass_n2*denom(k,i))+ > | qia(k,i,1) 784a870,872 > ! write(6,"('aurora: k=14,i=lon0,lat=',i3, > ! | ' QIA(k,i,lat) = ',f12.4,'QN2P_AUR = ',f12.4)") > ! | lat,qia(14,lon0,1),qn2p_aur(14,lon0) 800,806c888,900 < ! call addfsech('QO2P_AUR',' ',' ',qo2p_aur,lon0,lon1,nlevs,nlevs-1, < ! | lat) < ! call addfsech('QOP_AUR' ,' ',' ',qop_aur ,lon0,lon1,nlevs,nlevs-1, < ! | lat) < ! call addfsech('QN2P_AUR',' ',' ',qn2p_aur,lon0,lon1,nlevs,nlevs-1, < ! | lat) < --- > call addfsech('QO2P_AUR',' ',' ',qo2p_aur,lon0,lon1,nlevs,nlevs-1, > | lat) > call addfsech('QOP_AUR' ,' ',' ',qop_aur ,lon0,lon1,nlevs,nlevs-1, > | lat) > call addfsech('QN2P_AUR',' ',' ',qn2p_aur,lon0,lon1,nlevs,nlevs-1, > | lat) > call addfsech('QIA_O2P',' ',' ',qia(:,:,2),lon0,lon1,nlevs, > | nlevs-1,lat) > call addfsech('QIA_N2P',' ',' ',qia(:,:,1),lon0,lon1,nlevs, > | nlevs-1,lat) > call addfsech('QIA_OP',' ',' ',qia(:,:,3),lon0,lon1,nlevs,nlevs-1, > | lat) > ! ======================================================================== Diff of /home/tgcm/tiegcm1/bndry_gswm.F and /home/foster/tiegcm/amie/modsrc_oct03/bndry_gswm.F 125,127c125,127 < expt=cexp(ci*freq_semidi*dt*iter) < expt = cexp(ci*.5*freq_semidi*dt*iter) < expt = 1. --- > expt = cexp(ci*freq_semidi*dt*iter) > expt2 = cexp(ci*.5*freq_semidi*dt*iter) > expta = 1. 134c134 < z(i,j) = z(i,j)+real(zba(j)*bnda(i)*expt) ! annual tide --- > z(i,j) = z(i,j)+real(zba(j)*bnda(i)*expta) ! annual tide 141,142c141,142 < z(i,j) = z(i,j)+real(zb2(j)*bnd2(i)*expt)! diurnal tide < z(i,j) = z(i,j)+real(zba(j)*bnda(i)*expt) ! annual tide --- > z(i,j) = z(i,j)+real(zb2(j)*bnd2(i)*expt2)! diurnal tide > z(i,j) = z(i,j)+real(zba(j)*bnda(i)*expta) ! annual tide 148,150c148,150 < z(i,j) = real(zb(j)*bnd(i)*expt) ! semidiurnal tide < z(i,j) = z(i,j)+ zdi_gswm(i,j) ! diurnal tide < z(i,j) = z(i,j)+real(zba(j)*bnda(i)*expt) ! annual tide --- > z(i,j) = real(zb(j)*bnd(i)*expt) ! semidiurnal tide > z(i,j) = z(i,j)+ zdi_gswm(i,j) ! diurnal tide > z(i,j) = z(i,j)+real(zba(j)*bnda(i)*expta) ! annual tide 156,158c156,158 < z(i,j) = real(zb(j)*bnd(i)*expt) ! semidiurnal tide < z(i,j) = z(i,j)+real(zb2(j)*bnd2(i)*expt) ! diurnal tide < z(i,j) = z(i,j)+real(zba(j)*bnda(i)*expt) ! annual tide --- > z(i,j) = real(zb(j)*bnd(i)*expt) ! semidiurnal tide > z(i,j) = z(i,j)+real(zb2(j)*bnd2(i)*expt2) ! diurnal tide > z(i,j) = z(i,j)+real(zba(j)*bnda(i)*expta) ! annual tide ======================================================================== Diff of /home/tgcm/tiegcm1/chemrates.F and /home/foster/tiegcm/amie/modsrc_oct03/chemrates.F 162a163 > use magfield_module,only: rlatm 217a219,230 > ! reduce rk1 and rk2 to 20% in order to increase F-region Ne > ! G. Lu, September 1, 2003 > ! if (rlatm(i,lat) >= 0.1) then > ! rk1(k,i,lat) = .20*rk1(k,i,lat) > ! rk2(k,i,lat) = .20*rk2(k,i,lat) > ! else > ! rk1(k,i,lat) = .58*rk1(k,i,lat) > ! rk2(k,i,lat) = .58*rk2(k,i,lat) > ! endif > ! end of change in rk1 and rk2 > ! > ! ======================================================================== Diff of /home/tgcm/tiegcm1/cons.F and /home/foster/tiegcm/amie/modsrc_oct03/cons.F 127,128c127 < real,parameter :: < | crit(2) = (/0.261799387, 0.523598775/) --- > real :: crit(2) 254a254,258 > ! Set default crit. If reading amie data, this may be changed > ! (see aurora.F). > crit(1) = 0.261799387 > crit(2) = 0.523598775 > ! ======================================================================== >>> WARNING: Cannot find source file /home/tgcm/tiegcm1/dispose.F ======================================================================== Diff of /home/tgcm/tiegcm1/dynamics.F and /home/foster/tiegcm/amie/modsrc_oct03/dynamics.F 29a30 > use dyndiag_module,only: dyndiag,dyndiag_sech 229a231 > | z (levd0,lond0,lat,itp), 305a308 > | z (levd0,lond0,lat,itp), 362a366 > ! 619a624,644 > ! > ! Calculate 2d (lon,lat) diagnostics for secondary histories: > ! (see call dyndiag_sech below after lat loop) > ! > call dyndiag( > | ped (levd0,lond0,lat), ! pedersen conductivity > | hall(levd0,lond0,lat), ! hall conductivity > | qji_tn(levd0,lond0,lat), ! ion Joule heating for tn > | z (levd0,lond0,lat,itc), ! updated Z from addiag > | un (levd0,lond0,lat,itp), ! neutral zonal velocity > | vn (levd0,lond0,lat,itp), ! neutral meridional velocity > | ui (levd0,lond0,lat), ! ion zonal velocity > | vi (levd0,lond0,lat), ! ion meridional velocity > | 1,nlevp1,lon0,lon1,lat) > ! > ! Save dyndiags to secondary histories: > ! (these were calculated by sub dyndiag, see call above) > ! do lat=lat0,lat1 > ! call dyndiag_sech(lon0,lon1,1,nlevp1,lat) > ! enddo > ! 640c665 < | qji_tn(levd0,lond0,latd0), ! joule heating (sub qjoule_tn) --- > | qji_tn(levd0,lond0,latd0), ! joule heating (sub qjoule_tn.F) ======================================================================== Diff of /home/tgcm/tiegcm1/dynamo.F and /home/foster/tiegcm/amie/modsrc_oct03/dynamo.F 16,18d15 < ! Dymamo calls heelis to add heelis potential at high latitudes. < ! See heelis routines in heelis.F. < ! 107a105,111 > ! The dynamo is symmetric about the magnetic equator, but the high latitude > ! is anti-symmetric in both hemispheres. However, since Mudpack uses the > ! NH potential pattern, then the SH potential pattern must be added > ! back into the 2-D phim before the call threed, and before it is > ! transformed to geographic coordinates. > ! jn index for NH part of potential (nmlat down to ~nmlat0) > ! jp index for NH pfrac (nmlat0 down to 1), the fraction of the dynamo 109a114 > integer :: jn,jp 185a191,197 > ! For dot products: > real,parameter :: unitvm(nmlon)=1., unitv(nlon)=1. > ! > ! Electric potential from heelis or weimer: > real,dimension(nmlonp1,nmlat0) :: pfrac ! NH fraction of potential > real,dimension(nmlonp1,nmlat) :: phihm ! potential in magnetic > ! 290a303,422 > subroutine nodynamo > ! > ! 12/02: This routine solves for the 3-D electric potential and fields > ! when the dynamo flag is off (idynamo=0). Hence, the low-latitude > ! potential is zero, while the high-latitude potential is set > ! with potential_model='HEELIS','WEIMER', or 'NONE'. > ! > use input_module,only: potential_model > use cons_module,only: dlatm,dlonm,pi_dyn,ylatm,rtd > use magfield_module,only: im,jm,dim,djm > use fields_module,only: dynpot,phim3d > ! > ! Local: > integer :: i,j,k > real,dimension(nlonp1,0:nlatp1) :: phih ! potential in geographic > ! > ! Set phim=phihm of the chosen high-latitude model > do i=1,nmlonp1 > do j=1,nmlat > phim(i,j) = phihm(i,j) > enddo ! j=1,nmlat > enddo ! i=1,nmlonp1 > ! > ! transform geopotential height to mag. grid (used in subroutine threed) > call transf_dyn0 > ! > ! TEMP > ! do j=1,21,2 > ! write (6,"(1x,'j phihm =',i2,5e12.4)")j,ylatm(j)*rtd, > ! | (phihm(i,j),i=1,nmlon,20) > ! enddo > ! do j=51,nmlat,2 > ! write (6,"(1x,'j phihm =',i2,5e12.4)")j,ylatm(j)*rtd, > ! | (phihm(i,j),i=1,nmlon,20) > ! enddo > ! > ! Save single-level high latitude potential: > if (potential_model == 'WEIMER') then > call addfsech_ij('PHIHM2D','2D WEIMER01 ELECTRIC POTENTIAL', > | 'VOLTS', phihm,1,nmlonp1,1,nmlat) > else > if (potential_model == 'HEELIS') then > call addfsech_ij('PHIHM2D','2D HEELIS ELECTRIC POTENTIAL', > | 'VOLTS', phihm,1,nmlonp1,1,nmlat) > else > call addfsech_ij('PHIHM2D','2D ZERO ELECTRIC POTENTIAL', > | 'VOLTS', phihm,1,nmlonp1,1,nmlat) > endif > endif > ! > ! Save correct single-level potential to secondary histories (mag grid): > call addfsech_ij('PHIM2D','2D ELECTRIC POTENTIAL','VOLTS', > | phim,1,nmlonp1,1,nmlat) > ! > ! am_02/02: calculate current for both hemisphere: > ! R**2 * J_mr / dt0dts /rcos0s > if(icalkqlam.eq.1) call nosocrrt > ! > ! Call threed to generate 3-d potential array in geomagnetic coordinates > ! from 2-d solver output phim, corrected for the SH potential. > ! phim3d(nmlonp1,nmlat,-2:nlevp1) is in fields.F. > ! > call threed > > ! Save 3d mag electric potential: > ! phim3d(nmlonp1,nmlat,-2:nlevp1) ! 3d electric potential magnetic > ! > do j=1,nmlat > call addfsech_ik('PHIM3D','ELECTRIC POTENTIAL (MAG)','VOLTS', > | phim3d(:,j,:),1,nmlonp1,nmlev,nmlev-1,j) > enddo > ! > ! am_02/02: calculate Je_1, K_(q,phi), K_(q,lam) > if(icalkqlam.eq.1) call nosocrdens > ! > ! Transform phim3d to geographic coordinates in dynpot (fields.F): > ! phim3d(nmlonp1,nmlat,-2:nlevp1) ! 3d electric potential magnetic > ! dynpot(nlonp1,0:nlatp1,nlevp1), ! 3d electric potential geographic > ! > do k=1,nlevp1 > call mag2geo(phim3d(1,1,k),dynpot(1,0,k),im(1,0),jm(1,0), > | dim(1,0),djm(1,0),nlonp1,nmlonp1,nlon,nlat+2,nmlon,nmlat) > enddo ! k=1,nlevp1 > ! > ! Periodic point: > do k=1,nlevp1 > do j=0,nlatp1 > dynpot(nlonp1,j,k) = dynpot(1,j,k) > enddo ! j=0,nlatp1 > enddo ! k=1,nlevp1 > ! > ! Save electric potential on geographic coords to secondary history: > do j=1,nlat > call addfsech_ik('DYNPOT',' ',' ',dynpot(:,j,:), > | 1,nlonp1,nlevp1,nlevp1-1,j) > enddo ! j=1,nlat > ! > ! Transform single-level heelis magnetic potential phihm to geographic > ! in phih: > call mag2geo(phihm(1,1),phih(1,0),im(1,0),jm(1,0),dim(1,0), > | djm(1,0),nlonp1,nmlonp1,nlon,nlat+2,nmlon,nmlat) > ! > ! Periodic point: > do j=0,nlatp1 > phih(nlonp1,j) = phih(1,j) > enddo ! j=0,nlatp1 > ! > ! Save 2d heelis potential on geographic grid: > if (potential_model == 'WEIMER') then > call addfsech_ij('PHIH2D','2D WEIMER01 POTENTIAL (GEOG)','VOLTS' > | ,phih,1,nlonp1,1,nlat) > elseif (potential_model == 'HEELIS') then > call addfsech_ij('PHIH2D','2D HEELIS POTENTIAL (GEOG)','VOLTS', > | phih,1,nlonp1,1,nlat) > else > call addfsech_ij('PHIH2D','2D ZERO POTENTIAL (GEOG)','VOLTS', > | phih,1,nlonp1,1,nlat) > endif > end subroutine nodynamo > !----------------------------------------------------------------------- 292c424 < use cons_module,only: dlatm,dlonm,pi_dyn --- > use cons_module,only: dlatm,dlonm,pi_dyn,ylatm,rtd 294a427 > use input_module,only: potential_model,iamie 308c441 < integer :: i,j,jj,jjj,j0,jntl,k,n,ncc,nmaglat,nmaglon --- > integer :: i,j,jj,jjj,j0,jntl,k,n,ncc,nmaglat,nmaglon,ier 310a444,446 > real,dimension(nlonp1,0:nlatp1) :: phih ! potential in geographic > real,dimension(nmlonp1,nmlat0) :: phisym ! symmetric amie potential > real,dimension(nmlonp1,nmlat) :: phiasym ! asymmetric amie potential 312,317d447 < ! Heelis potential, and related: < real,dimension(nmlonp1,nmlat0) :: colatc ! colatitudes < real,dimension(nmlon0 ,nmlat0) :: pfrac ! fraction of potential < real,dimension(nmlonp1,nmlat) :: phihm ! potential in magnetic < real,dimension(nlonp1,0:nlatp1) :: phih ! potential in geographic < ! 516c646 < ! Get heelis colatitudes, pfrac, and potential phihm: --- > ! At this point, phihm is heelis, weimer, or amie (see advance.F) 518c648,658 < call heelis(colatc,pfrac,phihm,nmlat0,nmlon0,nmlat,nmlon,nmlonp1) --- > ! If not amie: > ! Modify stencils and RHS so that the NH high lat potential is inserted at > ! high latitude. The SH high lat potential will be added back later. > ! pfrac = fraction of dynamo in solution in the NH. = 1 low lat, = 0 hi lat > ! cons_module: crit(1)=15, crit(2)=30 deg colats, or hi-lat > 75 deg, > ! dynamo < 60 deg, and combination between 60-75 mag lat. > ! The dynamo is symmetric about the magnetic equator, but the high latitude > ! is anti-symmetric in both hemispheres. However, since Mudpack uses the > ! NH potential pattern, then the SH potential pattern must be added > ! back into the 2-D phim before the call threed, and before it is > ! transformed to geographic coordinates. 520,522c660,675 < ! Save single-level heelis potential: < call addfsech_ij('PHIHM2D','2d HEELIS ELECTRIC POTENTIAL','VOLTS', < | phihm,1,nmlonp1,1,nmlat) --- > if (iamie <= 0) then ! not amie run (heelis or weimer) > ncc = 1 > nmaglon = nmlon0 > nmaglat = nmlat0 > do n=1,5 > if (isolve==2) then > call stenmd(nmaglon,nmaglat,cee(ncc),phihm(1,nmlat0),pfrac) > else > call stenmod(nmaglon,nmaglat,cee(ncc),phihm(1,nmlat0),pfrac) > endif > ncc = ncc+9*nmaglon*nmaglat > if (n==1) ncc = ncc+nmaglon*nmaglat ! rhs is in 10th slot > nmaglon = (nmaglon+1)/2 > nmaglat = (nmaglat+1)/2 > enddo ! n=1,5 > phiasym = 0. 524,525c677,678 < ! Modify stencils and RHS so that heelis potential is inserted at < ! high latitude. --- > ! If amie, then break phihm into symmetric (phisym) and asymmetric (phiasym) > ! parts. phihm(nmlonp1,nmlat), phisym(nmlonp1,nmlat0), phyasym(nmlonp1,nmlat) 527,540c680,691 < ncc = 1 < nmaglon = nmlon0 < nmaglat = nmlat0 < do n=1,5 < if (isolve==2) then < call stenmd(nmaglon,nmaglat,cee(ncc),phihm(1,nmlat0),pfrac) < else < call stenmod(nmaglon,nmaglat,cee(ncc),phihm(1,nmlat0),pfrac) < endif < ncc = ncc+9*nmaglon*nmaglat < if (n==1) ncc = ncc+nmaglon*nmaglat ! rhs is in 10th slot < nmaglon = (nmaglon+1)/2 < nmaglat = (nmaglat+1)/2 < enddo ! n=1,5 --- > else ! is amie > do i=1,nmlonp1 > phisym(i,1) = phihm(i,nmlat0) > phiasym(i,nmlat0) = 0. > do j=1,nmlat0-1 > jj = nmlat+1-j > j0 = nmlat0-j+1 > phisym(i,j0) = 0.5*(phihm(i,j)+phihm(i,jj)) > phiasym(i,j) = 0.5*(phihm(i,j)+phihm(i,jj))-phihm(i,j) > phiasym(i,jj) = 0.5*(phihm(i,j)+phihm(i,jj))-phihm(i,jj) > enddo ! j=1,nmlat0-1 > enddo ! i=1,nmlonp1 541a693,721 > ! Modify stencils of symmetric part of amie potential: > ncc = 1 > nmaglon = nmlon0 > nmaglat = nmlat0 > do n=1,5 > if (isolve==2) then > call stenmd(nmaglon,nmaglat,cee(ncc),phisym(1,1), > | pfrac) > else > call stenmod(nmaglon,nmaglat,cee(ncc),phisym(1,1), > | pfrac) > endif > ncc = ncc+9*nmaglon*nmaglat > if (n==1) ncc = ncc+nmaglon*nmaglat ! rhs is in 10th slot > nmaglon = (nmaglon+1)/2 > nmaglat = (nmaglat+1)/2 > enddo ! n=1,5 > ! > ! p=0 from pole to crit1, p=1 from crit2 to equator: > do i=1,nmlonp1 > do j=1,nmlat0-1 > jj = nmlat+1-j > j0 = nmlat0-j+1 > phiasym(i,j) = phiasym(i,j) *(1.-pfrac(i,j0)) > phiasym(i,jj) = phiasym(i,jj)*(1.-pfrac(i,j0)) > enddo ! j=1,nmlat0-1 > enddo ! i=1,nmlonp1 > endif ! iamie > ! 551a732 > ier = 0 553c734,738 < call mud(rim,jntl,isolve) ! solver in mud.F --- > call mud(rim,jntl,isolve,ier) ! solver in mud.F > if(ier < 0 ) then ! not converged > write(6,*) 'muh: use direct solver' > call muh(rim,jntl) ! solver in mud.F > endif 555c740 < call muh(rim,jntl) ! solver in muh2cr.F --- > call muh(rim,jntl) ! solver in muh2cr.F 557c742,746 < call mudmod(rim,jntl) ! solver in mudmod.F --- > call mudmod(rim,jntl,ier) ! solver in mudmod.F > if(ier < 0 ) then ! not converged > write(6,*) 'muh: use direct solver' > call muh(rim,jntl) ! solver in mud.F > endif 564,570c753,759 < do j=1,nmlat < do i=1,nmlonp1 < phim(i,j) = rim(i,j,1) < enddo ! i=1,nmlonp1 < ! write(6,"('dynamo: j=',i3,' phim(:,j)=',/,(6e12.4))") < ! | j,phim(:,j) < enddo ! j=1,nmlat --- > ! Correct the SH potential for the anti-symmetric imposed NH high lat poten > ! pfrac = fraction of dynamo in solution in the NH. = 1 low lat, = 0 hi lat > ! cons_module: crit(1)=15, crit(2)=30 deg colats, or hi-lat > 75 deg, > ! dynamo < 60 deg, and combination between 60-75 mag lat. > ! ie, need (1.-pfrac(i,nmlat0)) at phim(i,1), etc > ! jn index for NH part of potential (nmlat down to ~nmlat0) > ! jp index for NH pfrac (nmlat0 down to 1) 572,573c761,804 < ! Save single-level potential to secondary histories (mag grid): < call addfsech_ij('PHIM2D','2d ELECTRIC POTENTIAL','VOLTS', --- > if (iamie <= 0) then ! not amie (heelis or weimer) > do j=1,nmlat0 > jn = nmlat - j + 1 > jp = nmlat0 - j + 1 > do i=1,nmlonp1 > phim(i,j)=rim(i,j,1)+(1.-pfrac(i,jp))* > | (phihm(i,j)-phihm(i,jn)) > enddo ! i=1,nmlonp1 > enddo ! j=1,nmlat0 > do j=nmlat0+1,nmlat > jp = j - nmlat0 > do i=1,nmlonp1 > phim(i,j) = rim(i,j,1) > enddo ! i=1,nmlonp1 > enddo ! j=1,nmlat > ! > ! Amie (fractional partitioning was done above): > else ! amie > do j=1,nmlat > do i=1,nmlonp1 > phim(i,j) = rim(i,j,1)-phiasym(i,j) > enddo ! i=1,nmlonp1 > enddo ! j=1,nmlat > endif > ! > ! Save single-level high latitude potential: > if (potential_model == 'WEIMER') then > call addfsech_ij('PHIHM2D','2D WEIMER01 ELECTRIC POTENTIAL', > | 'VOLTS', phihm,1,nmlonp1,1,nmlat) > elseif (potential_model == 'HEELIS') then > call addfsech_ij('PHIHM2D','2D HEELIS ELECTRIC POTENTIAL', > | 'VOLTS', phihm,1,nmlonp1,1,nmlat) > else > call addfsech_ij('PHIHM2D','2D ZERO ELECTRIC POTENTIAL', > | 'VOLTS', phihm,1,nmlonp1,1,nmlat) > endif > ! > ! Save correct single-level potential to secondary histories (mag grid): > call addfsech_ij('RIM',' ',' ',rim(:,:,1),1,nmlonp1,1,nmlat) > call addfsech_ij('PHISYM','symmetric AMIE potential','Volts ', > | phisym,1,nmlonp1,1,nmlat) > call addfsech_ij('PHIASYM','asymmetric AMIE potential','Volts ', > | phiasym,1,nmlonp1,1,nmlat) > call addfsech_ij('PHIM2D','2D ELECTRIC POTENTIAL','VOLTS', 581,582c812,813 < ! from 2-d solver output phim. phim3d(nmlonp1,nmlat,-2:nlevp1) is in < ! fields.F. --- > ! from 2-d solver output phim, corrected for the SH potential. > ! phim3d(nmlonp1,nmlat,-2:nlevp1) is in fields.F. 629,631c860,870 < ! Save 2d heelis potential on geographic grid: < call addfsech_ij('PHIH2D','2D HEELIS POTENTIAL (GEOG)','VOLTS', < | phih,1,nlonp1,1,nlat) --- > ! Save 2d potential on geographic grid: > if (potential_model == 'WEIMER') then > call addfsech_ij('PHIH2D','2D WEIMER01 POTENTIAL (GEOG)','VOLTS' > | ,phih,1,nlonp1,1,nlat) > elseif (potential_model == 'HEELIS') then > call addfsech_ij('PHIH2D','2D HEELIS POTENTIAL (GEOG)','VOLTS', > | phih,1,nlonp1,1,nlat) > else > call addfsech_ij('PHIH2D','2D ZERO POTENTIAL (GEOG)','VOLTS', > | phih,1,nlonp1,1,nlat) > endif 649d887 < real,parameter :: unitv(nmlon)=1. 1281,1282c1519,1520 < | (4.*sddot(nmlon,unitv,zpotenm3d(1,2,k))- < | sddot(nmlon,unitv,zpotenm3d(1,3,k)))/(3.*float(nmlon)) --- > | (4.*sddot(nmlon,unitvm,zpotenm3d(1,2,k))- > | sddot(nmlon,unitvm,zpotenm3d(1,3,k)))/(3.*float(nmlon)) 1284,1285c1522,1523 < | (4.*sddot(nmlon,unitv,zpotenm3d(1,nmlat-1,k))- < | sddot(nmlon,unitv,zpotenm3d(1,nmlat-2,k)))/ --- > | (4.*sddot(nmlon,unitvm,zpotenm3d(1,nmlat-1,k))- > | sddot(nmlon,unitvm,zpotenm3d(1,nmlat-2,k)))/ 1302,1317c1540,1555 < zigm11(1, 1) = (4.*sddot(nmlon,unitv,zigm11(1, 2))- < 1 sddot(nmlon,unitv,zigm11(1, 3)))/(3.*float(nmlon)) < zigm11(1,nmlat) = (4.*sddot(nmlon,unitv,zigm11(1,nmlat-1))- < 1 sddot(nmlon,unitv,zigm11(1,nmlat-2)))/(3.*float(nmlon)) < zigmc(1, 1) = (4.*sddot(nmlon,unitv,zigmc(1, 2))- < 1 sddot(nmlon,unitv,zigmc(1, 3)))/(3.*float(nmlon)) < zigmc(1,nmlat) = (4.*sddot(nmlon,unitv,zigmc(1,nmlat-1))- < 1 sddot(nmlon,unitv,zigmc(1,nmlat-2)))/(3.*float(nmlon)) < zigm2(1, 1) = (4.*sddot(nmlon,unitv,zigm2(1, 2))- < 1 sddot(nmlon,unitv,zigm2(1, 3)))/(3.*float(nmlon)) < zigm2(1,nmlat) = (4.*sddot(nmlon,unitv,zigm2(1,nmlat-1))- < 1 sddot(nmlon,unitv,zigm2(1,nmlat-2)))/(3.*float(nmlon)) < zigm22(1, 1) = (4.*sddot(nmlon,unitv,zigm22(1, 2))- < 1 sddot(nmlon,unitv,zigm22(1, 3)))/(3.*float(nmlon)) < zigm22(1,nmlat) = (4.*sddot(nmlon,unitv,zigm22(1,nmlat-1))- < 1 sddot(nmlon,unitv,zigm22(1,nmlat-2)))/(3.*float(nmlon)) --- > zigm11(1, 1) = (4.*sddot(nmlon,unitvm,zigm11(1, 2))- > 1 sddot(nmlon,unitvm,zigm11(1, 3)))/(3.*float(nmlon)) > zigm11(1,nmlat) = (4.*sddot(nmlon,unitvm,zigm11(1,nmlat-1))- > 1 sddot(nmlon,unitvm,zigm11(1,nmlat-2)))/(3.*float(nmlon)) > zigmc(1, 1) = (4.*sddot(nmlon,unitvm,zigmc(1, 2))- > 1 sddot(nmlon,unitvm,zigmc(1, 3)))/(3.*float(nmlon)) > zigmc(1,nmlat) = (4.*sddot(nmlon,unitvm,zigmc(1,nmlat-1))- > 1 sddot(nmlon,unitvm,zigmc(1,nmlat-2)))/(3.*float(nmlon)) > zigm2(1, 1) = (4.*sddot(nmlon,unitvm,zigm2(1, 2))- > 1 sddot(nmlon,unitvm,zigm2(1, 3)))/(3.*float(nmlon)) > zigm2(1,nmlat) = (4.*sddot(nmlon,unitvm,zigm2(1,nmlat-1))- > 1 sddot(nmlon,unitvm,zigm2(1,nmlat-2)))/(3.*float(nmlon)) > zigm22(1, 1) = (4.*sddot(nmlon,unitvm,zigm22(1, 2))- > 1 sddot(nmlon,unitvm,zigm22(1, 3)))/(3.*float(nmlon)) > zigm22(1,nmlat) = (4.*sddot(nmlon,unitvm,zigm22(1,nmlat-1))- > 1 sddot(nmlon,unitvm,zigm22(1,nmlat-2)))/(3.*float(nmlon)) 1548d1785 < real,parameter :: unitv(nmlon)=1. 1597c1834 < | sddot(nmlon,unitv,rim(1,nmlat-1,2))/tint1(nmlat-1) --- > | sddot(nmlon,unitvm,rim(1,nmlat-1,2))/tint1(nmlat-1) 2352d2588 < real,parameter :: unitv(nmlon)=1. 2413a2650,2652 > > call addfsech_ij('ED1','2D ELECTRIC FIELD','VOLTS/M', > | ed1,1,nmlonp1,1,nmlat) 2477,2478c2716,2717 < phims = sddot(nmlon,unitv,phim(1,1))/float(nmlon) < phimn = sddot(nmlon,unitv,phim(1,nmlat))/float(nmlon) --- > phims = sddot(nmlon,unitvm,phim(1,1))/float(nmlon) > phimn = sddot(nmlon,unitvm,phim(1,nmlat))/float(nmlon) 2686a2926,3136 > subroutine transf_dyn0 > ! am_11/03 > ! used only if dynamo = 0 > ! need zpoten3d in magnetic coordinates for subroutine threed > ! > use cons_module,only: h0,r0 > use magfield_module,only: > | zb, ! (nlonp1,0:nlatp1) > | ig, ! (nmlonp1,nmlat) geog lon at each geomag grid point > | jg, ! (nmlonp1,nmlat) geog lat at each geomag grid point > | wt ! (4,nmlonp1,nmlat) interpolation weights > ! > ! Local: > integer :: i,ii,k,kk,j,jj,lat,n > real :: z0 > ! > ! Fields to be transformed to geomagnetic space (formerly in transmag.h): > ! (these will be input to the geographic to magnetic transformation). > real,dimension(nlonp1,0:nlatp1,-2:nlevp1) :: zz > ! > real,dimension(nlonp4,nlevp1) :: zz_plt > ! > ! External: > real,external :: sddot ! in util.F > ! > ! Set constants: > ! z0 is lowest level for start of field line integration set in h0 > ! (h0 in cons module) > ! > ! > ! Pack inputs from 3->nlon+3 to 1->nlon+1, as in tgcm15 (end of sub lamdas): > do lat=1,nlat > do i=1,nlon+1 > ii = i+2 > zpoten(i,lat,:) = zpoten(ii,lat,:) > enddo ! i=1,nlonp4-3 > enddo ! lat=1,nlat > ! > z0 = h0 > ! > ! Calculate quantities to be transformed to geomagnetic space: > do k=1,nlev > do j=1,nlat > do i=1,nlonp1 > zz(i,j,k) = zpoten(i,j,k) > enddo ! i=1,nlonp1 > enddo ! j=1,nlat > enddo ! k=1,nlev > do j=1,nlat > do i=1,nlonp1 > zz(i,j,nlevp1) = zpoten(i,j,nlevp1) > enddo ! i=1,nlonp1 > enddo ! j=1,nlat > ! > ! Extend fields down to 90 km inserting 3 extra levels. > ! Set three equally spaced levels for Z, take U, V, and W > ! to be constant, and extrapolate sigmas exponentially. > ! > do k=0,-2,-1 > do j=1,nlat > do i=1,nlonp1 > zz(i,j,k) = z0+float(k+2)*(zz(i,j,1)-z0)/3. > enddo ! i=1,nlonp1 > enddo ! j=1,nlat > enddo ! k=0,-2,-1 > ! > ! Values at poles: > do k=-2,nlev > zz(1,0,k) = > | (9.*sddot(nlon,unitv,zz(1,1,k))- > | sddot(nlon,unitv,zz(1,2,k)))/(8.*float(nlon)) > zz(1,nlatp1,k) = > | (9.*sddot(nlon,unitv,zz(1,nlat ,k))- > | sddot(nlon,unitv,zz(1,nlat-1,k)))/(8.*float(nlon)) > ! > ! Extend in longitude: > do i = 2,nlon > zz(i,0,k) = zz(1,0,k) > zz(i,nlatp1,k) = zz(1,nlatp1,k) > enddo ! i = 2,nlon > enddo ! k=-2,nlev > ! > ! Values at the poles: > zz(1,0,nlevp1)= (9.*sddot(nlon,unitv,zz(1,1,nlevp1))- > | sddot(nlon,unitv,zz(1,2,nlevp1)))/ > | (8.*float(nlon)) > zz(1,nlatp1,nlevp1) = (9.* > | sddot(nlon,unitv,zz(1,nlat,nlevp1))- > | sddot(nlon,unitv,zz(1,nlat-1,nlevp1)))/ > | (8.*float(nlon)) > ! > ! Extend in longitude: > do i = 2,nlon > zz(i,0,nlevp1) = zz(1,0,nlevp1) > zz(i,nlatp1,nlevp1)= zz(1,nlatp1,nlevp1) > enddo ! i = 2,nlon > ! > ! Periodic points: > do j = 0,nlatp1 > do k = -2,nlev > zz(nlonp1,j,k) = zz(1,j,k) > enddo ! k = -2,nlev > zz(nlonp1,j,nlevp1)= zz(1,j,nlevp1) > enddo ! j = 0,nlatp1 > ! > ! Transform needed fields to geomagnetic coordinate system, one latitude > ! at a time. > ! > ! subroutine geo2mag(fmag,fgeo,long,latg,wght,nlonp1_geo,nlonp1_mag, > !| nlon_mag,nlat_mag,lat) > ! > ! Magnetic latitude loop: > maglat_loop: do j=2,nmlat-1 > do k=-2,nlevp1 > call geo2mag(zpotenm(1,k),zz(1,0,k),ig,jg,wt,nlonp1, > | nmlonp1,nmlon,nmlat,j) > enddo ! k=-2,nlevp1 > ! Without calculation of k_(q,lam) > do i=1,nmlon > do k=-2,nlev > if (zpotenm(i,k) < z0) zpotenm(i,k) = z0 > zpotenm3d(i,j,k) = zpotenm(i,k) > enddo ! k=-2,nlev > if (zpotenm(i,nlevp1) < z0) zpotenm(i,nlevp1) = z0 > zpotenm3d(i,j,nlevp1) = zpotenm(i,nlevp1) > enddo ! i=1,nmlon > enddo maglat_loop ! j=2,nmlat-1 Main magnetic latitude loop > ! > ! Polar values for Z: > do k=-2,nlevp1 > zpotenm3d(1,1,k) = > | (4.*sddot(nmlon,unitvm,zpotenm3d(1,2,k))- > | sddot(nmlon,unitvm,zpotenm3d(1,3,k)))/(3.*float(nmlon)) > zpotenm3d(1,nmlat,k) = > | (4.*sddot(nmlon,unitvm,zpotenm3d(1,nmlat-1,k))- > | sddot(nmlon,unitvm,zpotenm3d(1,nmlat-2,k)))/ > | (3.*float(nmlon)) > ! > ! Extend Z over longitude: > do i=1,nmlon > zpotenm3d(i,1,k) = zpotenm3d(1,1,k) > zpotenm3d(i,nmlat,k) = zpotenm3d(1,nmlat,k) > enddo ! i=1,nmlon > ! > ! Periodic points: > do j=1,nmlat > zpotenm3d(nmlonp1,j,k) = zpotenm3d(1,j,k) > enddo ! j=1,nmlat > enddo ! k=-2,nlevp1 > ! > ! Save 3d potential on magnetic grid to secondary history: > ! do j=1,nmlat > ! call addfsech_ik('ZPOTEN3D','ZPOTEN3D','[cm]',zpotenm3d(:,j,:), > ! | 1,nmlonp1,nmlev,nmlev,j) > ! enddo ! j=1,nmlat > ! > end subroutine transf_dyn0 > !----------------------------------------------------------------------- > subroutine prep_dynamo_dyn0(z,lev0,lev1,lon0,lon1,lat0,lat1) > ! am_11/03: if dynamo = 0 need zpoten3d for subroutine threed > ! Prepare geographic-grid fields for input to the dynamo, and gather them > ! to the root task. This is executed by all tasks, and is called from > ! advance before the dynamo itself (which is executed by master task only). > ! > use cons_module,only: gask,grav > ! > ! Args: > integer :: lev0,lev1,lon0,lon1,lat0,lat1 > real :: fmin,fmax > ! > ! Input fields at task subdomains. These are from current (not updated) > ! fields (itp): > real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),intent(in):: > | z ! geopotential height > ! > ! Local: > integer :: k,i,lat,nk,lonbeg,lonend > ! > nk = lev1-lev0+1 > lonbeg = lon0 > if (lon0==1) lonbeg = 3 > lonend = lon1 > if (lon1==nlonp4) lonend = nlonp4-1 > ! > ! Define subdomain part of dynamo input fields (geographic): > ! Define subdomain part of dynamo input fields (geographic): > ! Also transform from (k,i,lat) to (i,lat,k). The latter is used > ! in the dynamo code, whereas the former is used in the rest > ! of the model. > ! (In earlier versions, this was in lamdas, where 1-73 <= 3-75. > ! Here, 3-75 <= 3-75) > ! > do lat=lat0,lat1 > do i=lonbeg,lonend > do k=lev0,lev1-1 > zpoten (i,lat,k) = z (k,i,lat) > enddo ! k=lev0,lev1 > zpoten(i,lat,lev1) = z(lev1,i,lat) > enddo ! i=lon0,lon1 > ! > enddo ! lat=lat0,lat1 > > #ifdef MPI > ! > ! Gather dynamo input fields to the root task, defining module data > ! above at the global domain on the root task: > ! > call mp_dynamo_gather > #endif > end subroutine prep_dynamo_dyn0 > !----------------------------------------------------------- ======================================================================== Diff of /home/tgcm/tiegcm1/dynamo_old.F and /home/foster/tiegcm/amie/modsrc_oct03/dynamo_old.F 104a105,111 > ! The dynamo is symmetric about the magnetic equator, but the high latitude > ! is anti-symmetric in both hemispheres. However, since Mudpack uses the > ! NH potential pattern, then the SH potential pattern must be added > ! back into the 2-D phim before the call threed, and before it is > ! transformed to geographic coordinates. > ! jn index for NH part of potential (nmlat down to ~nmlat0) > ! jp index for NH pfrac (nmlat0 down to 1), the fraction of the dynamo 106a114 > integer :: jn,jp 146a155,162 > ! > ! For dot products: > real,parameter :: unitvm(nmlon)=1., unitv(nlon)=1. > ! > ! Electric potential from heelis or weimer: > real,dimension(nmlonp1,nmlat0) :: > | pfrac, ! NH fraction of potential > | phihm ! potential in magnetic 157c173 < use cons_module,only: pi_dyn,dlonm,dlatm --- > use cons_module,only: pi_dyn,dlonm,dlatm,ylatm,rtd 159a176 > use input_module,only: potential_model 167,171d183 < ! < ! Heelis potential, and related: < real,dimension(nmlonp1,nmlat0) :: colatc ! colatitudes < real,dimension(nmlon0 ,nmlat0) :: pfrac ! fraction of potential < real,dimension(nmlonp1,nmlat) :: phihm ! potential in magnetic 275,278d286 < ! Get heelis colatitudes, pfrac, and potential phihm: < ! < call heelis(colatc,pfrac,phihm,nmlat0,nmlon0,nmlat,nmlon,nmlonp1) < ! 280,283c288,291 < do j=1,nmlat < do i=1,nmlonp1 < phihm_plt(i,:) = phihm(i,j) < enddo ! i=1,nmlonp1 --- > ! do j=1,nmlat > ! do i=1,nmlonp1 > ! phihm_plt(i,:) = phihm(i,j) > ! enddo ! i=1,nmlonp1 286c294 < enddo ! j=1,nmlat --- > ! enddo ! j=1,nmlat 289c297,305 < ! high latitude. --- > ! high latitude. The SH high lat potential will be added back later. > ! pfrac = fraction of dynamo in solution in the NH. = 1 low lat, = 0 hi lat > ! cons_module: crit(1)=15, crit(2)=30 deg colats, or hi-lat > 75 deg, > ! dynamo < 60 deg, and combination between 60-75 mag lat. > ! The dynamo is symmetric about the magnetic equator, but the high latitude > ! is anti-symmetric in both hemispheres. However, since Mudpack uses the > ! NH potential pattern, then the SH potential pattern must be added > ! back into the 2-D phim before the call threed, and before it is > ! transformed to geographic coordinates. 308c324,333 < do j=1,nmlat --- > ! Correct the SH potential for the anti-symmetric imposed NH high lat poten > ! pfrac = fraction of dynamo in solution in the NH. = 1 low lat, = 0 hi lat > ! cons_module: crit(1)=15, crit(2)=30 deg colats, or hi-lat > 75 deg, > ! dynamo < 60 deg, and combination between 60-75 mag lat. > ! ie, need (1.-pfrac(i,nmlat0)) at phim(i,1), etc > ! jn index for NH part of potential (nmlat down to ~nmlat0) > ! jp index for NH pfrac (nmlat0 down to 1) > do j=1,nmlat0 > jn = nmlat - j + 1 > jp = nmlat0 - j + 1 309a335,343 > phim(i,j)=rim(i,j,1)+(1.-pfrac(i,jp))*(phihm(i,j)-phihm(i,jn)) > enddo ! i=1,nmlonp1 > ! write(6,"('dynamo: j=',3i3,' rim phim s,n,frac=',/,(5e12.4))") > ! | j,jn,jp,(rim(i,j,1),phim(i,j),phihm(i,j),phihm(i,jn), > ! | pfrac(i,jp),i=1,nmlonp1,20) > enddo ! j=1,nmlat0 > do j=nmlat0+1,nmlat > jp = j - nmlat0 > do i=1,nmlonp1 312,313c346,347 < ! write(6,"('dynamo: j=',i3,' phim(:,j)=',/,(6e12.4))") < ! | j,phim(:,j) --- > ! write(6,"('dynamo: j=',2i3,' phim n,frac=',/,(6e12.4))") > ! | j,jp,(phim(i,j),phihm(i,j),pfrac(i,jp),i=1,nmlonp1,20) 316c350 < ! Save single-level potential to secondary histories (mag grid): --- > ! Save correct single-level potential to secondary histories (mag grid): 319c353 < ! phim_plt(i,:) = phim(i,j) --- > ! phim_plt(i,:) = phimsn(i,j) 325,326c359,361 < ! Generate 3-d potential array phim3d in geomagnetic coordinates from < ! 2-d mudpack output phim. phim3d is in fields.F. --- > ! Call threed to generate 3-d potential array in geomagnetic coordinates > ! from 2-d solver output phim, corrected for the SH potential. > ! phim3d(nmlonp1,nmlat,-2:nlevp1) is in fields.F. 869c904 < real,parameter :: unitv(nmlon)=1. --- > real :: fmin,fmax 1070a1106,1113 > ! write(6,"(/,'k=',i3,' nlon=',i3,' ssigma1(1,1,k)=',e12.4, > ! | ' (1,2,k)=',e12.4)") k,nlon,ssigma1(1,1,k),ssigma1(1,2,k) > ! write(6,"('k=',i3,' nlat=',i3,' ssigma1(1,nlat,k)=',e12.4, > ! | ' (1,nlat-1,k)=',e12.4)") k,nlat, > ! | ssigma1(1,nlat,k),ssigma1(1,nlat-1,k) > ! write(6,"('ssigma1(1,0,k)=',e12.4,' ssigma1(1,nlatp1,k)=', > ! | e12.4)") ssigma1(1,0,k),ssigma1(1,nlatp1,k) > 1279,1289c1322,1332 < ! do k=-2,nlev < ! sigma_pedm (nmlonp1,k) = sigma_pedm (1,k) < ! sigma_hallm(nmlonp1,k) = sigma_hallm(1,k) < ! wnvelm (nmlonp1,k) = wnvelm (1,k) < ! vxbm (nmlonp1,k) = vxbm (1,k) < ! zpotenm (nmlonp1,k) = zpotenm (1,k) < ! adotvm (nmlonp1,k,1) = adotvm (1,k,1) < ! adotvm (nmlonp1,k,2) = adotvm (1,k,2) < ! axvm (nmlonp1,k,1) = axvm (1,k,1) < ! axvm (nmlonp1,k,2) = axvm (1,k,2) < ! enddo ! k=-2,nlev --- > do k=-2,nlev > sigma_pedm (nmlonp1,k) = sigma_pedm (1,k) > sigma_hallm(nmlonp1,k) = sigma_hallm(1,k) > wnvelm (nmlonp1,k) = wnvelm (1,k) > vxbm (nmlonp1,k) = vxbm (1,k) > zpotenm (nmlonp1,k) = zpotenm (1,k) > adotvm (nmlonp1,k,1) = adotvm (1,k,1) > adotvm (nmlonp1,k,2) = adotvm (1,k,2) > axvm (nmlonp1,k,1) = axvm (1,k,1) > axvm (nmlonp1,k,2) = axvm (1,k,2) > enddo ! k=-2,nlev 1291,1301d1333 < ! call addfsech_ik('SIGMA1M' ,' ',' ',sigma_pedm, < ! | 1,nmlonp1,nmlev,nmlev-1,j) < ! call addfsech_ik('SIGMA2M' ,' ',' ',sigma_hallm, < ! | 1,nmlonp1,nmlev,nmlev-1,j) < ! call addfsech_ik('WNVELM' ,' ',' ',wnvelm, < ! | 1,nmlonp1,nmlev,nmlev-1,j) < ! call addfsech_ik('VXBM' ,' ',' ',vxbm, < ! | 1,nmlonp1,nmlev,nmlev-1,j) < ! call addfsech_ik('ZPOTENM' ,' ',' ',zpotenm, < ! | 1,nmlonp1,nmlev,nmlev-1,j) < ! 1320,1321c1352,1353 < | (4.*sddot(nmlon,unitv,zpotenm3d(1,2,k))- < | sddot(nmlon,unitv,zpotenm3d(1,3,k)))/(3.*float(nmlon)) --- > | (4.*sddot(nmlon,unitvm,zpotenm3d(1,2,k))- > | sddot(nmlon,unitvm,zpotenm3d(1,3,k)))/(3.*float(nmlon)) 1323,1324c1355,1356 < | (4.*sddot(nmlon,unitv,zpotenm3d(1,nmlat-1,k))- < | sddot(nmlon,unitv,zpotenm3d(1,nmlat-2,k)))/ --- > | (4.*sddot(nmlon,unitvm,zpotenm3d(1,nmlat-1,k))- > | sddot(nmlon,unitvm,zpotenm3d(1,nmlat-2,k)))/ 1361,1368c1393,1400 < zigm11(1,nmlat) = (4.*sddot(nmlon,unitv,zigm11(1,nmlat-1))- < | sddot(nmlon,unitv,zigm11(1,nmlat-2)))/(3.*float(nmlon)) < zigmc (1,nmlat) = (4.*sddot(nmlon,unitv,zigmc (1,nmlat-1))- < | sddot(nmlon,unitv,zigmc (1,nmlat-2)))/(3.*float(nmlon)) < zigm2 (1,nmlat) = (4.*sddot(nmlon,unitv,zigm2 (1,nmlat-1))- < | sddot(nmlon,unitv,zigm2 (1,nmlat-2)))/(3.*float(nmlon)) < zigm22(1,nmlat) = (4.*sddot(nmlon,unitv,zigm22(1,nmlat-1))- < | sddot(nmlon,unitv,zigm22(1,nmlat-2)))/(3.*float(nmlon)) --- > zigm11(1,nmlat) = (4.*sddot(nmlon,unitvm,zigm11(1,nmlat-1))- > | sddot(nmlon,unitvm,zigm11(1,nmlat-2)))/(3.*float(nmlon)) > zigmc (1,nmlat) = (4.*sddot(nmlon,unitvm,zigmc (1,nmlat-1))- > | sddot(nmlon,unitvm,zigmc (1,nmlat-2)))/(3.*float(nmlon)) > zigm2 (1,nmlat) = (4.*sddot(nmlon,unitvm,zigm2 (1,nmlat-1))- > | sddot(nmlon,unitvm,zigm2 (1,nmlat-2)))/(3.*float(nmlon)) > zigm22(1,nmlat) = (4.*sddot(nmlon,unitvm,zigm22(1,nmlat-1))- > | sddot(nmlon,unitvm,zigm22(1,nmlat-2)))/(3.*float(nmlon)) 1465c1497 < | sddot(nmlon,unitv,rim(1,nmlat-1,2))/tint1(nmlat-1) --- > | sddot(nmlon,unitvm,rim(1,nmlat-1,2))/tint1(nmlat-1) 1864d1895 < real,parameter :: unitv(nmlon)=1. 1905,1906c1936,1937 < phims = sddot(nmlon,unitv,phim(1,1))/float(nmlon) < phimn = sddot(nmlon,unitv,phim(1,nmlat))/float(nmlon) --- > phims = sddot(nmlon,unitvm,phim(1,1))/float(nmlon) > phimn = sddot(nmlon,unitvm,phim(1,nmlat))/float(nmlon) ======================================================================== >>> WARNING: Cannot find source file /home/tgcm/tiegcm1/dyndiag.F ======================================================================== Diff of /home/tgcm/tiegcm1/elden.F and /home/foster/tiegcm/amie/modsrc_oct03/elden.F 3c3 < | nplus,n2p,nop,o2p,electrons,lev0,lev1,lon0,lon1,lat) --- > | z,nplus,n2p,nop,o2p,electrons,lev0,lev1,lon0,lon1,lat) 7a8 > use params_module,only: nlevp1 12a14 > use dyndiag_module,only: tec_sec,tec 30c32,33 < | xiop2d ! from oplus --- > | xiop2d, ! from oplus > | z ! geopotenital height 49c52 < --- > ! 218a222,230 > C Get ht integral of electron density and save in tec of /sight/ > tec(:,lat) = 0. > do i=lon0,lon1 > do k=lev0,lev1-1 > tec_sec(i,lat,k) = electrons(k,i) > tec(i,lat) = tec(i,lat) +(z(k+1,i)-z(k,i))*electrons(k,i) > enddo > tec_sec(i,lat,nlevp1) = tec(i,lat) > enddo 220,221c232,233 < ! call addfsech('NE_ELDEN',' ',' ',electrons(:,i0:i1), < ! | i0,i1,nk,nk,lat) --- > call addfsech('NE_ELDEN',' ',' ',electrons(:,i0:i1), > | i0,i1,nk,nk,lat) ======================================================================== >>> WARNING: Cannot find source file /home/tgcm/tiegcm1/extraoutput.F ======================================================================== Diff of /home/tgcm/tiegcm1/filter.F and /home/foster/tiegcm/amie/modsrc_oct03/filter.F 8c8 < ! Setfft calls set99, which returns trigs and ifix. --- > ! Setfft calls set99, which returns trigs and ifax. 17d16 < ! 82c81 < integer,parameter :: nlond4 = nlon/4 --- > integer,parameter :: nlond4 = nlon/4, nlond4_5 = 72/4 84c83 < integer :: nn_5(72/4) ! for 5 deg in lon --- > integer :: nn_5(nlond4_5) ! for 5 deg in lon 89a89 > ! nn_5=(/90,40,22,14,10,8,6,4,2,1,1,1,1,1,1,1,1,1/) ! 5 deg grid 94,95c94,97 < elseif (dlon==2.5) then ! 2.5 deg res < do i=1,nlond4 --- > ! > ! For 2.5 degree grid: > elseif (dlon==2.5) then > do i=1,nlond4_5 99,100d100 < ! write(6,"('filter2: nlon=',i3,' nlon/4=',i3,' nn=',/, < ! | (12i4))") nlon,nlon/4,nn 105d104 < ======================================================================== Diff of /home/tgcm/tiegcm1/getms.F and /home/foster/tiegcm/amie/modsrc_oct03/getms.F 82c82 < write(6,"('getms: MPI master task: mspath = ',a)") trim(mspath) --- > ! write(6,"('getms: MPI master task: mspath = ',a)") trim(mspath) 114c114 < return --- > goto 100 ======================================================================== Diff of /home/tgcm/tiegcm1/gswm.F and /home/foster/tiegcm/amie/modsrc_oct03/gswm.F 3c3,4 < use params_module,only: nlon,nlat,nlonp4 --- > use params_module,only: nlon,nlat,nlonp4,nlonp2 > use mpi_module,only: lon0,lon1,lat0,lat1 12a14,17 > ! am 11/27/02: for using the double TIEGCM resolution grid together with the > ! GSWM input, the GSWM grid is set to the double resolution grid and for > ! single TIEGCM resolution only every second GSWM grid point is used > ! 26c31,35 < | mxgswmhour = 24 ! maximum number of hours of gswm data --- > | mxgswmhour = 24, ! maximum number of hours of gswm data > ! | mxgswmlon = 144,! maximum number of longitudinal grid points > ! | mxgswmlat = 72 ! maximum number of latitudinal grid points > | mxgswmlon = 72,! maximum number of longitudinal grid points > | mxgswmlat = 36 ! maximum number of latitudinal grid points 29a39 > ! 46a57 > ! 68c79 < integer :: istat,ncid,i,j,ier --- > integer :: istat,ncid,i,j,ier,i2,j2 70c81,87 < | idv_un,idv_vn,ngswmlon,ngswmlat --- > | idv_un,idv_vn,ngswmlon,ngswmlat,lonbeg,lonend,latbeg, > | lonbeg_dhz,latbeg_dhz > real :: > | gswm_z(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour), ! diurnal tides data org.grid > | gswm_tn(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour), ! diurnal tides data org.grid > | gswm_un(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour), ! diurnal tides data org.grid > | gswm_vn(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour) ! diurnal tides data org.grid 116c133 < if (ngswmlon /= nlon) then --- > if (ngswmlon /= mxgswmlon) then 118,119c135,136 < | 'nlon(GSWM)=',i8,' nlon(TIEGCM)=',i8)") ngswmlon, < | nlon --- > | 'nlon(GSWM)=',i8,' mxgswmlon=',i8)") ngswmlon, > | mxgswmlon 128c145 < if (ngswmlat /= nlat) then --- > if (ngswmlat /= mxgswmlat) then 130,131c147,148 < | 'nlat(GSWM)=',i8,' nlat(TIEGCM)=',i8)") ngswmlat, < | nlat --- > | 'nlat(GSWM)=',i8,' mxgswmlat=',i8)") ngswmlat, > | mxgswmlat 134a152 > ! read whole domain 139c157 < istat = nf_get_var_double(ncid,idv_z,gswm_di_z) --- > istat = nf_get_var_double(ncid,idv_z,gswm_z) 142c160 < gswm_di_z = gswm_di_z*100. ! convert to cm --- > gswm_z = gswm_z*100. ! convert to cm 148c166 < istat = nf_get_var_double(ncid,idv_tn,gswm_di_tn) --- > istat = nf_get_var_double(ncid,idv_tn,gswm_tn) 156c174 < istat = nf_get_var_double(ncid,idv_un,gswm_di_un) --- > istat = nf_get_var_double(ncid,idv_un,gswm_un) 159c177 < gswm_di_un = gswm_di_un*100. ! convert to cm/s --- > gswm_un = gswm_un*100. ! convert to cm/s 165c183 < istat = nf_get_var_double(ncid,idv_vn,gswm_di_vn) --- > istat = nf_get_var_double(ncid,idv_vn,gswm_vn) 168c186 < gswm_di_vn = gswm_di_vn*100. ! convert to cm/s --- > gswm_vn = gswm_vn*100. ! convert to cm/s 173a192,227 > ! > ! check for the TIEGCM dimension of the grid > ! if necessary copy data to right sized arrays > ! copy into subdomain > ! > lonbeg = lon0 > if (lon0==1) lonbeg = 3 > lonend = lon1 > if (lon1==nlonp4) lonend = lon1-2 > > if(nlat == mxgswmlat) then ! same dimension nlon=144,nlat=72 > do i = lonbeg,lonend > do j = lat0,lat1 > gswm_di_z(i,j,:,:) = gswm_z(i-2,j,:,:) > gswm_di_tn(i,j,:,:) = gswm_tn(i-2,j,:,:) > gswm_di_un(i,j,:,:) = gswm_un(i-2,j,:,:) > gswm_di_vn(i,j,:,:) = gswm_vn(i-2,j,:,:) > enddo > enddo > else ! half dimension nlon=72,nlat=36 > > lonbeg_dhz = lonbeg*2-1 > latbeg_dhz = lat0*2-1 > i2 = lonbeg_dhz > do i = lonbeg,lonend > j2 = latbeg_dhz > do j = lat0,lat1 > gswm_di_z(i,j,:,:) = gswm_z(i2-4,j2,:,:) > gswm_di_tn(i,j,:,:) = gswm_tn(i2-4,j2,:,:) > gswm_di_un(i,j,:,:) = gswm_un(i2-4,j2,:,:) > gswm_di_vn(i,j,:,:) = gswm_vn(i2-4,j2,:,:) > j2 = j2+2 > enddo > i2 = i2+2 > enddo > endif 188c242 < integer :: istat,ncid,i,j,ier --- > integer :: istat,ncid,i,j,ier,i2,j2 190c244,250 < | idv_un,idv_vn,ngswmlon,ngswmlat --- > | idv_un,idv_vn,ngswmlon,ngswmlat,lonbeg,lonend,latbeg, > | lonbeg_dhz,latbeg_dhz > real :: > | gswm_z(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour), ! semidiurnal tides data org.grid > | gswm_tn(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour), ! semidiurnal tides data org.grid > | gswm_un(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour), ! semidiurnal tides data org.grid > | gswm_vn(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour) ! semidiurnal tides data org.grid 236c296 < if (ngswmlon /= nlon) then --- > if (ngswmlon /= mxgswmlon) then 238,239c298,299 < | 'nlon(GSWM)=',i8,' nlon(TIEGCM)=',i8)") ngswmlon, < | nlon --- > | 'nlon(GSWM)=',i8,' mxgswmlon=',i8)") ngswmlon, > | mxgswmlon 248c308 < if (ngswmlat /= nlat) then --- > if (ngswmlat /= mxgswmlat) then 250,251c310,311 < | 'nlat(GSWM)=',i8,' nlat(TIEGCM)=',i8)") ngswmlat, < | nlat --- > | 'nlat(GSWM)=',i8,' mxgswmlat=',i8)") ngswmlat, > | mxgswmlat 254a315 > ! read whole domain 259c320 < istat = nf_get_var_double(ncid,idv_z,gswm_sdi_z) --- > istat = nf_get_var_double(ncid,idv_z,gswm_z) 262c323 < gswm_sdi_z = gswm_sdi_z*100. ! convert to cm --- > gswm_z = gswm_z*100. ! convert to cm 268c329 < istat = nf_get_var_double(ncid,idv_tn,gswm_sdi_tn) --- > istat = nf_get_var_double(ncid,idv_tn,gswm_tn) 276c337 < istat = nf_get_var_double(ncid,idv_un,gswm_sdi_un) --- > istat = nf_get_var_double(ncid,idv_un,gswm_un) 279c340 < gswm_sdi_un = gswm_sdi_un*100. ! convert to cm/s --- > gswm_un = gswm_un*100. ! convert to cm/s 285c346 < istat = nf_get_var_double(ncid,idv_vn,gswm_sdi_vn) --- > istat = nf_get_var_double(ncid,idv_vn,gswm_vn) 288c349 < gswm_sdi_vn = gswm_sdi_vn*100. ! convert to cm/s --- > gswm_vn = gswm_vn*100. ! convert to cm/s 293a355,390 > > ! check for the TIEGCM dimension of the grid > ! if necessary copy data to right sized arrays > ! copy into subdomain > ! > lonbeg = lon0 > if (lon0==1) lonbeg = 3 > lonend = lon1 > if (lon1==nlonp4) lonend = lon1-2 > > if(nlat == mxgswmlat) then ! same dimension nlon=144,nlat=72 > do i = lonbeg,lonend > do j = lat0,lat1 > gswm_sdi_z(i,j,:,:) = gswm_z(i-2,j,:,:) > gswm_sdi_tn(i,j,:,:) = gswm_tn(i-2,j,:,:) > gswm_sdi_un(i,j,:,:) = gswm_un(i-2,j,:,:) > gswm_sdi_vn(i,j,:,:) = gswm_vn(i-2,j,:,:) > enddo > enddo > else ! half dimension nlon=72,nlat=36 > > lonbeg_dhz = lonbeg*2-1 > latbeg_dhz = lat0*2-1 > i2 = lonbeg_dhz > do i = lonbeg,lonend > j2 = latbeg_dhz > do j = lat0,lat1 > gswm_sdi_z(i,j,:,:) = gswm_z(i2-4,j2,:,:) > gswm_sdi_tn(i,j,:,:) = gswm_tn(i2-4,j2,:,:) > gswm_sdi_un(i,j,:,:) = gswm_un(i2-4,j2,:,:) > gswm_sdi_vn(i,j,:,:) = gswm_vn(i2-4,j2,:,:) > j2 = j2+2 > enddo > i2 = i2+2 > enddo > endif 308c405 < integer :: istat,ncid,i,j,ier --- > integer :: istat,ncid,i,j,ier,i2,j2 310c407,413 < | idv_un,idv_vn,ngswmlon,ngswmlat --- > | idv_un,idv_vn,ngswmlon,ngswmlat,lonbeg,lonend,latbeg, > | lonbeg_dhz,latbeg_dhz > real :: ! org. GSWM grid > | gswm_z(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour), ! nonmigrating diurnal tides data > | gswm_tn(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour), ! nonmigrating diurnal tides data > | gswm_un(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour), ! nonmigrating diurnal tides data > | gswm_vn(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour) ! nonmigrating diurnal tides data 357c460 < if (ngswmlon /= nlon) then --- > if (ngswmlon /= mxgswmlon) then 359,360c462,463 < | 'nlon(GSWM)=',i8,' nlon(TIEGCM)=',i8)") ngswmlon, < | nlon --- > | 'nlon(GSWM)=',i8,' mxgswmlon=',i8)") ngswmlon, > | mxgswmlon 369c472 < if (ngswmlat /= nlat) then --- > if (ngswmlat /= mxgswmlat) then 371,372c474,475 < | 'nlat(GSWM)=',i8,' nlat(TIEGCM)=',i8)") ngswmlat, < | nlat --- > | 'nlat(GSWM)=',i8,' mxgswmlat=',i8)") ngswmlat, > | mxgswmlat 375a479 > ! read whole domain 380c484 < istat = nf_get_var_double(ncid,idv_z,gswm_nmidi_z) --- > istat = nf_get_var_double(ncid,idv_z,gswm_z) 383c487 < gswm_nmidi_z = gswm_nmidi_z*100. ! convert to cm --- > gswm_z = gswm_z*100. ! convert to cm 389c493 < istat = nf_get_var_double(ncid,idv_tn,gswm_nmidi_tn) --- > istat = nf_get_var_double(ncid,idv_tn,gswm_tn) 397c501 < istat = nf_get_var_double(ncid,idv_un,gswm_nmidi_un) --- > istat = nf_get_var_double(ncid,idv_un,gswm_un) 400c504 < gswm_nmidi_un = gswm_nmidi_un*100. ! convert to cm/s --- > gswm_un = gswm_un*100. ! convert to cm/s 406c510 < istat = nf_get_var_double(ncid,idv_vn,gswm_nmidi_vn) --- > istat = nf_get_var_double(ncid,idv_vn,gswm_vn) 409c513 < gswm_nmidi_vn = gswm_nmidi_vn*100. ! convert to cm/s --- > gswm_vn = gswm_vn*100. ! convert to cm/s 414a519,555 > > ! check for the TIEGCM dimension of the grid > ! if necessary copy data to right sized arrays > ! copy into subdomain > ! > lonbeg = lon0 > if (lon0==1) lonbeg = 3 > lonend = lon1 > if (lon1==nlonp4) lonend = lon1-2 > > if(nlat == mxgswmlat) then ! same dimension nlon=144,nlat=72 > do i = lonbeg,lonend > do j = lat0,lat1 > gswm_nmidi_z(i,j,:,:) = gswm_z(i-2,j,:,:) > gswm_nmidi_tn(i,j,:,:) = gswm_tn(i-2,j,:,:) > gswm_nmidi_un(i,j,:,:) = gswm_un(i-2,j,:,:) > gswm_nmidi_vn(i,j,:,:) = gswm_vn(i-2,j,:,:) > enddo > enddo > else ! half dimension nlon=72,nlat=36 > > lonbeg_dhz = lonbeg*2-1 > latbeg_dhz = lat0*2-1 > i2 = lonbeg_dhz > > do i = lon0,lon1 > j2 = latbeg_dhz > do j = lat0,lat1 > gswm_nmidi_z(i,j,:,:) = gswm_z(i2-4,j2,:,:) > gswm_nmidi_tn(i,j,:,:) = gswm_tn(i2-4,j2,:,:) > gswm_nmidi_un(i,j,:,:) = gswm_un(i2-4,j2,:,:) > gswm_nmidi_vn(i,j,:,:) = gswm_vn(i2-4,j2,:,:) > j2 = j2+2 > enddo > i2 = i2+2 > enddo > endif 430c571 < integer :: istat,ncid,i,j,ier --- > integer :: istat,ncid,i,j,ier,i2,j2 432c573,579 < | idv_un,idv_vn,ngswmlon,ngswmlat --- > | idv_un,idv_vn,ngswmlon,ngswmlat,lonbeg,lonend,latbeg, > | lonbeg_dhz,latbeg_dhz > real :: > | gswm_z(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour), ! nonmigrating semidiurnal tides data > | gswm_tn(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour),! nonmigrating semidiurnal tides data > | gswm_un(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour),! nonmigrating semidiurnal tides data > | gswm_vn(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour)! nonmigrating semidiurnal tides data 479c626 < if (ngswmlon /= nlon) then --- > if (ngswmlon /= mxgswmlon) then 481,482c628,629 < | 'fit: nlon(GSWM)=',i8,' nlon(TIEGCM)=',i8)") ngswmlon, < | nlon --- > | 'fit: nlon(GSWM)=',i8,' mxgswmlon=',i8)") ngswmlon, > | mxgswmlon 491c638 < if (ngswmlat /= nlat) then --- > if (ngswmlat /= mxgswmlat) then 493,494c640,641 < | 'nlat(GSWM)=',i8,' nlat(TIEGCM)=',i8)") ngswmlat, < | nlat --- > | 'nlat(GSWM)=',i8,' mxgswmlat=',i8)") ngswmlat, > | mxgswmlat 497a645 > ! read whole domain 502c650 < istat = nf_get_var_double(ncid,idv_z,gswm_nmisdi_z) --- > istat = nf_get_var_double(ncid,idv_z,gswm_z) 505c653 < gswm_nmisdi_z = gswm_nmisdi_z*100. ! convert to cm --- > gswm_z = gswm_z*100. ! convert to cm 511c659 < istat = nf_get_var_double(ncid,idv_tn,gswm_nmisdi_tn) --- > istat = nf_get_var_double(ncid,idv_tn,gswm_tn) 519c667 < istat = nf_get_var_double(ncid,idv_un,gswm_nmisdi_un) --- > istat = nf_get_var_double(ncid,idv_un,gswm_un) 522c670 < gswm_nmisdi_un = gswm_nmisdi_un*100. ! convert to cm/s --- > gswm_un = gswm_un*100. ! convert to cm/s 528c676 < istat = nf_get_var_double(ncid,idv_vn,gswm_nmisdi_vn) --- > istat = nf_get_var_double(ncid,idv_vn,gswm_vn) 531c679 < gswm_nmisdi_vn = gswm_nmisdi_vn*100. ! convert to cm/s --- > gswm_vn = gswm_vn*100. ! convert to cm/s 536a685,719 > > ! check for the TIEGCM dimension of the grid > ! if necessary copy data to right sized arrays > ! > lonbeg = lon0 > if (lon0==1) lonbeg = 3 > lonend = lon1 > if (lon1==nlonp4) lonend = lon1-2 > > if(nlat == mxgswmlat) then ! same dimension nlon=144,nlat=72 > do i = lonbeg,lonend > do j = lat0,lat1 > gswm_nmisdi_z(i,j,:,:) = gswm_z(i-2,j,:,:) > gswm_nmisdi_tn(i,j,:,:) = gswm_tn(i-2,j,:,:) > gswm_nmisdi_un(i,j,:,:) = gswm_un(i-2,j,:,:) > gswm_nmisdi_vn(i,j,:,:) = gswm_vn(i-2,j,:,:) > enddo > enddo > else ! half dimension nlon=72,nlat=36 > > lonbeg_dhz = lonbeg*2-1 > latbeg_dhz = lat0*2-1 > i2 = lonbeg_dhz > do i = lonbeg,lonend > j2 = latbeg_dhz > do j = lat0,lat1 > gswm_nmisdi_z(i,j,:,:) = gswm_z(i2-4,j2,:,:) > gswm_nmisdi_tn(i,j,:,:) = gswm_tn(i2-4,j2,:,:) > gswm_nmisdi_un(i,j,:,:) = gswm_un(i2-4,j2,:,:) > gswm_nmisdi_vn(i,j,:,:) = gswm_vn(i2-4,j2,:,:) > j2 = j2+2 > enddo > i2 = i2+2 > enddo > endif 551a735,737 > #ifdef MPI > use mpi_module,only: mp_periodic_f2d > #endif 558c744,745 < | mon_cur,mon_nxt,ihr_cur,ihr_nxt,nointp=0 --- > | mon_cur,mon_nxt,ihr_cur,ihr_nxt,nointp=0, > | lonbeg,lonend 560,565c747,758 < | z_curmo(nlon,nlat)=0.,z_nxtmo(nlon,nlat)=0., ! UT interpolation < | tn_curmo(nlon,nlat)=0.,tn_nxtmo(nlon,nlat)=0., ! UT interpolation < | un_curmo(nlon,nlat)=0.,un_nxtmo(nlon,nlat)=0., ! UT interpolation < | vn_curmo(nlon,nlat)=0.,vn_nxtmo(nlon,nlat)=0., ! UT interpolation < | z_tmp(nlon,nlat)=0.,tn_tmp(nlon,nlat)=0., ! tmp < | un_tmp(nlon,nlat)=0.,vn_tmp(nlon,nlat)=0. ! tmp --- > | z_curmo(lon0:lon1,lat0:lat1), > | z_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation > | tn_curmo(lon0:lon1,lat0:lat1), > | tn_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation > | un_curmo(lon0:lon1,lat0:lat1), > | un_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation > | vn_curmo(lon0:lon1,lat0:lat1), > | vn_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation > | z_tmp(lon0:lon1,lat0:lat1), > | tn_tmp(lon0:lon1,lat0:lat1), ! tmp > | un_tmp(lon0:lon1,lat0:lat1), > | vn_tmp(lon0:lon1,lat0:lat1) ! tmp 616,623c809,816 < z_curmo(:,:) = gswm_di_z(:,:,mon_cur,ihr_cur+1) ! current month < tn_curmo(:,:) = gswm_di_tn(:,:,mon_cur,ihr_cur+1) < un_curmo(:,:) = gswm_di_un(:,:,mon_cur,ihr_cur+1) < vn_curmo(:,:) = gswm_di_vn(:,:,mon_cur,ihr_cur+1) < z_nxtmo(:,:) = gswm_di_z(:,:,mon_nxt,ihr_cur+1) ! next month < tn_nxtmo(:,:) = gswm_di_tn(:,:,mon_nxt,ihr_cur+1) < un_nxtmo(:,:) = gswm_di_un(:,:,mon_nxt,ihr_cur+1) < vn_nxtmo(:,:) = gswm_di_vn(:,:,mon_nxt,ihr_cur+1) --- > z_curmo(:,:) =gswm_di_z (lon0:lon1,lat0:lat1,mon_cur,ihr_cur+1) ! current month > tn_curmo(:,:)=gswm_di_tn(lon0:lon1,lat0:lat1,mon_cur,ihr_cur+1) > un_curmo(:,:)=gswm_di_un(lon0:lon1,lat0:lat1,mon_cur,ihr_cur+1) > vn_curmo(:,:)=gswm_di_vn(lon0:lon1,lat0:lat1,mon_cur,ihr_cur+1) > z_nxtmo(:,:) =gswm_di_z (lon0:lon1,lat0:lat1,mon_nxt,ihr_cur+1) ! next month > tn_nxtmo(:,:)=gswm_di_tn(lon0:lon1,lat0:lat1,mon_nxt,ihr_cur+1) > un_nxtmo(:,:)=gswm_di_un(lon0:lon1,lat0:lat1,mon_nxt,ihr_cur+1) > vn_nxtmo(:,:)=gswm_di_vn(lon0:lon1,lat0:lat1,mon_nxt,ihr_cur+1) 628,651c821,844 < call timeliper_array(gswm_di_z(:,:,mon_cur,ihr_cur+1), ! cur. month < | gswm_di_z(:,:,mon_cur,ihr_nxt+1),difsec, < | secint,z_curmo,1) < call timeliper_array(gswm_di_tn(:,:,mon_cur,ihr_cur+1), ! cur. month < | gswm_di_tn(:,:,mon_cur,ihr_nxt+1),difsec, < | secint,tn_curmo,1) < call timeliper_array(gswm_di_un(:,:,mon_cur,ihr_cur+1), ! cur. month < | gswm_di_un(:,:,mon_cur,ihr_nxt+1),difsec, < | secint,un_curmo,1) < call timeliper_array(gswm_di_vn(:,:,mon_cur,ihr_cur+1), ! cur. month < | gswm_di_vn(:,:,mon_cur,ihr_nxt+1),difsec, < | secint,vn_curmo,1) < call timeliper_array(gswm_di_z(:,:,mon_nxt,ihr_cur+1), ! next month < | gswm_di_z(:,:,mon_nxt,ihr_nxt+1),difsec, < | secint,z_nxtmo,1) < call timeliper_array(gswm_di_tn(:,:,mon_nxt,ihr_cur+1), ! next month < | gswm_di_tn(:,:,mon_nxt,ihr_nxt+1),difsec, < | secint,tn_nxtmo,1) < call timeliper_array(gswm_di_un(:,:,mon_nxt,ihr_cur+1), ! next month < | gswm_di_un(:,:,mon_nxt,ihr_nxt+1),difsec, < | secint,un_nxtmo,1) < call timeliper_array(gswm_di_vn(:,:,mon_nxt,ihr_cur+1), ! next month < | gswm_di_vn(:,:,mon_nxt,ihr_nxt+1),difsec, < | secint,vn_nxtmo,1) --- > call timeliper_array(gswm_di_z(lon0:lon1,lat0:lat1,mon_cur, ! cur. month > | ihr_cur+1),gswm_di_z(lon0:lon1,lat0:lat1,mon_cur,ihr_nxt+1), > | difsec,secint,z_curmo,1) > call timeliper_array(gswm_di_tn(lon0:lon1,lat0:lat1,mon_cur, ! cur. month > | ihr_cur+1),gswm_di_tn(lon0:lon1,lat0:lat1,mon_cur,ihr_nxt+1), > | difsec,secint,tn_curmo,1) > call timeliper_array(gswm_di_un(lon0:lon1,lat0:lat1,mon_cur, ! cur. month > | ihr_cur+1),gswm_di_un(lon0:lon1,lat0:lat1,mon_cur,ihr_nxt+1), > | difsec,secint,un_curmo,1) > call timeliper_array(gswm_di_vn(lon0:lon1,lat0:lat1,mon_cur, ! cur. month > | ihr_cur+1),gswm_di_vn(lon0:lon1,lat0:lat1,mon_cur,ihr_nxt+1), > | difsec,secint,vn_curmo,1) > call timeliper_array(gswm_di_z(lon0:lon1,lat0:lat1,mon_nxt, ! next month > | ihr_cur+1),gswm_di_z(lon0:lon1,lat0:lat1,mon_nxt,ihr_nxt+1), > | difsec,secint,z_nxtmo,1) > call timeliper_array(gswm_di_tn(lon0:lon1,lat0:lat1,mon_nxt, ! next month > | ihr_cur+1),gswm_di_tn(lon0:lon1,lat0:lat1,mon_nxt,ihr_nxt+1), > | difsec,secint,tn_nxtmo,1) > call timeliper_array(gswm_di_un(lon0:lon1,lat0:lat1,mon_nxt, ! next month > | ihr_cur+1),gswm_di_un(lon0:lon1,lat0:lat1,mon_nxt,ihr_nxt+1), > | difsec,secint,un_nxtmo,1) > call timeliper_array(gswm_di_vn(lon0:lon1,lat0:lat1,mon_nxt, ! next month > | ihr_cur+1),gswm_di_vn(lon0:lon1,lat0:lat1,mon_nxt,ihr_nxt+1), > | difsec,secint,vn_nxtmo,1) 692,695c885,888 < zdi_gswm(3:74,:) = z_curmo < tndi_gswm(3:74,:) = tn_curmo < undi_gswm(3:74,:) = un_curmo < vndi_gswm(3:74,:) = vn_curmo --- > zdi_gswm (lon0:lon1,lat0:lat1) = z_curmo > tndi_gswm(lon0:lon1,lat0:lat1) = tn_curmo > undi_gswm(lon0:lon1,lat0:lat1) = un_curmo > vndi_gswm(lon0:lon1,lat0:lat1) = vn_curmo 698,701c891,894 < zdi_gswm(3:74,:) = z_nxtmo < tndi_gswm(3:74,:) = tn_nxtmo < undi_gswm(3:74,:) = un_nxtmo < vndi_gswm(3:74,:) = vn_nxtmo --- > zdi_gswm (lon0:lon1,lat0:lat1) = z_nxtmo > tndi_gswm(lon0:lon1,lat0:lat1) = tn_nxtmo > undi_gswm(lon0:lon1,lat0:lat1) = un_nxtmo > vndi_gswm(lon0:lon1,lat0:lat1) = vn_nxtmo 705,717c898,902 < zdi_gswm(3:74,:) = z_tmp < tndi_gswm(3:74,:) = tn_tmp < undi_gswm(3:74,:) = un_tmp < vndi_gswm(3:74,:) = vn_tmp < ! periodic points < 30 zdi_gswm(1:2,:) = z_tmp(71:72,:) < tndi_gswm(1:2,:) = tn_tmp(71:72,:) < undi_gswm(1:2,:) = un_tmp(71:72,:) < vndi_gswm(1:2,:) = vn_tmp(71:72,:) < zdi_gswm(75:76,:) = z_tmp(1:2,:) < tndi_gswm(75:76,:) = tn_tmp(1:2,:) < undi_gswm(75:76,:) = un_tmp(1:2,:) < vndi_gswm(75:76,:) = vn_tmp(1:2,:) --- > zdi_gswm (lon0:lon1,lat0:lat1) = z_tmp > tndi_gswm(lon0:lon1,lat0:lat1) = tn_tmp > undi_gswm(lon0:lon1,lat0:lat1) = un_tmp > vndi_gswm(lon0:lon1,lat0:lat1) = vn_tmp > 30 continue 718a904,924 > ! do mpi periodic points exchange for gswm with f2d(:) > ! Moved here from sub getgswm* because mpi necessary when f2d data > ! is allocated only for task-local subdomain block. > ! lons 1,2 <- nlonp4-3,nlonp4-2 and nlonp4-1,nlonp4 <- 3,4 > > #ifdef MPI > call mp_periodic_f2d(zdi_gswm(lon0:lon1,lat0:lat1), > | lon0,lon1,lat0,lat1) > call mp_periodic_f2d(tndi_gswm(lon0:lon1,lat0:lat1), > | lon0,lon1,lat0,lat1) > call mp_periodic_f2d(undi_gswm(lon0:lon1,lat0:lat1), > | lon0,lon1,lat0,lat1) > call mp_periodic_f2d(vndi_gswm(lon0:lon1,lat0:lat1), > | lon0,lon1,lat0,lat1) > # else > set_periodic_f2d(zdi_gswm) > set_periodic_f2d(tndi_gswm) > set_periodic_f2d(undi_gswm) > set_periodic_f2d(vndi_gswm) > #endif > 736a943,945 > #ifdef MPI > use mpi_module,only: mp_periodic_f2d > #endif 743c952,953 < | mon_cur,mon_nxt,ihr_cur,ihr_nxt,nointp=0 --- > | mon_cur,mon_nxt,ihr_cur,ihr_nxt,nointp=0, > | lonbeg,lonend 745,750c955,966 < | z_curmo(nlon,nlat)=0.,z_nxtmo(nlon,nlat)=0., ! UT interpolation < | tn_curmo(nlon,nlat)=0.,tn_nxtmo(nlon,nlat)=0., ! UT interpolation < | un_curmo(nlon,nlat)=0.,un_nxtmo(nlon,nlat)=0., ! UT interpolation < | vn_curmo(nlon,nlat)=0.,vn_nxtmo(nlon,nlat)=0., ! UT interpolation < | z_tmp(nlon,nlat)=0.,tn_tmp(nlon,nlat)=0., ! tmp < | un_tmp(nlon,nlat)=0.,vn_tmp(nlon,nlat)=0. ! tmp --- > | z_curmo(lon0:lon1,lat0:lat1), > | z_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation > | tn_curmo(lon0:lon1,lat0:lat1), > | tn_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation > | un_curmo(lon0:lon1,lat0:lat1), > | un_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation > | vn_curmo(lon0:lon1,lat0:lat1), > | vn_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation > | z_tmp(lon0:lon1,lat0:lat1), > | tn_tmp(lon0:lon1,lat0:lat1), ! tmp > | un_tmp(lon0:lon1,lat0:lat1), > | vn_tmp(lon0:lon1,lat0:lat1) ! tmp 800,807c1016,1031 < z_curmo(:,:) = gswm_sdi_z(:,:,mon_cur,ihr_cur+1) ! current month < tn_curmo(:,:) = gswm_sdi_tn(:,:,mon_cur,ihr_cur+1) < un_curmo(:,:) = gswm_sdi_un(:,:,mon_cur,ihr_cur+1) < vn_curmo(:,:) = gswm_sdi_vn(:,:,mon_cur,ihr_cur+1) < z_nxtmo(:,:) = gswm_sdi_z(:,:,mon_nxt,ihr_cur+1) ! next month < tn_nxtmo(:,:) = gswm_sdi_tn(:,:,mon_nxt,ihr_cur+1) < un_nxtmo(:,:) = gswm_sdi_un(:,:,mon_nxt,ihr_cur+1) < vn_nxtmo(:,:) = gswm_sdi_vn(:,:,mon_nxt,ihr_cur+1) --- > z_curmo(:,:) =gswm_sdi_z (lon0:lon1,lat0:lat1, ! current month > | mon_cur,ihr_cur+1) > tn_curmo(:,:)=gswm_sdi_tn(lon0:lon1,lat0:lat1, > | mon_cur,ihr_cur+1) > un_curmo(:,:)=gswm_sdi_un(lon0:lon1,lat0:lat1, > | mon_cur,ihr_cur+1) > vn_curmo(:,:)=gswm_sdi_vn(lon0:lon1,lat0:lat1, > | mon_cur,ihr_cur+1) > z_nxtmo(:,:) =gswm_sdi_z (lon0:lon1,lat0:lat1, ! next month > | mon_nxt,ihr_cur+1) > tn_nxtmo(:,:)=gswm_sdi_tn(lon0:lon1,lat0:lat1, > | mon_nxt,ihr_cur+1) > un_nxtmo(:,:)=gswm_sdi_un(lon0:lon1,lat0:lat1, > | mon_nxt,ihr_cur+1) > vn_nxtmo(:,:)=gswm_sdi_vn(lon0:lon1,lat0:lat1, > | mon_nxt,ihr_cur+1) 812,835c1036,1059 < call timeliper_array(gswm_sdi_z(:,:,mon_cur,ihr_cur+1), ! cur. month < | gswm_sdi_z(:,:,mon_cur,ihr_nxt+1),difsec, < | secint,z_curmo,1) < call timeliper_array(gswm_sdi_tn(:,:,mon_cur,ihr_cur+1), ! cur. month < | gswm_sdi_tn(:,:,mon_cur,ihr_nxt+1),difsec, < | secint,tn_curmo,1) < call timeliper_array(gswm_sdi_un(:,:,mon_cur,ihr_cur+1), ! cur. month < | gswm_sdi_un(:,:,mon_cur,ihr_nxt+1),difsec, < | secint,un_curmo,1) < call timeliper_array(gswm_sdi_vn(:,:,mon_cur,ihr_cur+1), ! cur. month < | gswm_sdi_vn(:,:,mon_cur,ihr_nxt+1),difsec, < | secint,vn_curmo,1) < call timeliper_array(gswm_sdi_z(:,:,mon_nxt,ihr_cur+1), ! next month < | gswm_sdi_z(:,:,mon_nxt,ihr_nxt+1),difsec, < | secint,z_nxtmo,1) < call timeliper_array(gswm_sdi_tn(:,:,mon_nxt,ihr_cur+1), ! next month < | gswm_sdi_tn(:,:,mon_nxt,ihr_nxt+1),difsec, < | secint,tn_nxtmo,1) < call timeliper_array(gswm_sdi_un(:,:,mon_nxt,ihr_cur+1), ! next month < | gswm_sdi_un(:,:,mon_nxt,ihr_nxt+1),difsec, < | secint,un_nxtmo,1) < call timeliper_array(gswm_sdi_vn(:,:,mon_nxt,ihr_cur+1), ! next month < | gswm_sdi_vn(:,:,mon_nxt,ihr_nxt+1),difsec, < | secint,vn_nxtmo,1) --- > call timeliper_array(gswm_sdi_z(lon0:lon1,lat0:lat1,mon_cur, ! cur. month > | ihr_cur+1),gswm_sdi_z(lon0:lon1,lat0:lat1,mon_cur,ihr_nxt+1), > | difsec,secint,z_curmo,1) > call timeliper_array(gswm_sdi_tn(lon0:lon1,lat0:lat1,mon_cur, ! cur. month > | ihr_cur+1),gswm_sdi_tn(lon0:lon1,lat0:lat1,mon_cur,ihr_nxt+1), > | difsec,secint,tn_curmo,1) > call timeliper_array(gswm_sdi_un(lon0:lon1,lat0:lat1,mon_cur, ! cur. month > | ihr_cur+1),gswm_sdi_un(lon0:lon1,lat0:lat1,mon_cur,ihr_nxt+1), > | difsec,secint,un_curmo,1) > call timeliper_array(gswm_sdi_vn(lon0:lon1,lat0:lat1,mon_cur, ! cur. month > | ihr_cur+1),gswm_sdi_vn(lon0:lon1,lat0:lat1,mon_cur,ihr_nxt+1), > | difsec,secint,vn_curmo,1) > call timeliper_array(gswm_sdi_z(lon0:lon1,lat0:lat1,mon_nxt, ! next month > | ihr_cur+1),gswm_sdi_z(lon0:lon1,lat0:lat1,mon_nxt,ihr_nxt+1), > | difsec,secint,z_nxtmo,1) > call timeliper_array(gswm_sdi_tn(lon0:lon1,lat0:lat1,mon_nxt, ! next month > | ihr_cur+1),gswm_sdi_tn(lon0:lon1,lat0:lat1,mon_nxt,ihr_nxt+1), > | difsec,secint,tn_nxtmo,1) > call timeliper_array(gswm_sdi_un(lon0:lon1,lat0:lat1,mon_nxt, ! next month > | ihr_cur+1),gswm_sdi_un(lon0:lon1,lat0:lat1,mon_nxt,ihr_nxt+1), > | difsec,secint,un_nxtmo,1) > call timeliper_array(gswm_sdi_vn(lon0:lon1,lat0:lat1,mon_nxt, ! next month > | ihr_cur+1),gswm_sdi_vn(lon0:lon1,lat0:lat1,mon_nxt,ihr_nxt+1), > | difsec,secint,vn_nxtmo,1) 875,878c1099,1102 < zsdi_gswm(3:74,:) = z_curmo < tnsdi_gswm(3:74,:) = tn_curmo < unsdi_gswm(3:74,:) = un_curmo < vnsdi_gswm(3:74,:) = vn_curmo --- > zsdi_gswm (lon0:lon1,lat0:lat1) = z_curmo > tnsdi_gswm(lon0:lon1,lat0:lat1) = tn_curmo > unsdi_gswm(lon0:lon1,lat0:lat1) = un_curmo > vnsdi_gswm(lon0:lon1,lat0:lat1) = vn_curmo 881,884c1105,1108 < zsdi_gswm(3:74,:) = z_nxtmo < tnsdi_gswm(3:74,:) = tn_nxtmo < unsdi_gswm(3:74,:) = un_nxtmo < vnsdi_gswm(3:74,:) = vn_nxtmo --- > zsdi_gswm (lon0:lon1,lat0:lat1) = z_nxtmo > tnsdi_gswm(lon0:lon1,lat0:lat1) = tn_nxtmo > unsdi_gswm(lon0:lon1,lat0:lat1) = un_nxtmo > vnsdi_gswm(lon0:lon1,lat0:lat1) = vn_nxtmo 888,901c1112,1116 < zsdi_gswm(3:74,:) = z_tmp < tnsdi_gswm(3:74,:) = tn_tmp < unsdi_gswm(3:74,:) = un_tmp < vnsdi_gswm(3:74,:) = vn_tmp < ! < ! periodic points < 30 zsdi_gswm(1:2,:) = z_tmp(71:72,:) < tnsdi_gswm(1:2,:) = tn_tmp(71:72,:) < unsdi_gswm(1:2,:) = un_tmp(71:72,:) < vnsdi_gswm(1:2,:) = vn_tmp(71:72,:) < zsdi_gswm(75:76,:) = z_tmp(1:2,:) < tnsdi_gswm(75:76,:) = tn_tmp(1:2,:) < unsdi_gswm(75:76,:) = un_tmp(1:2,:) < vnsdi_gswm(75:76,:) = vn_tmp(1:2,:) --- > zsdi_gswm (lon0:lon1,lat0:lat1) = z_tmp > tnsdi_gswm(lon0:lon1,lat0:lat1) = tn_tmp > unsdi_gswm(lon0:lon1,lat0:lat1) = un_tmp > vnsdi_gswm(lon0:lon1,lat0:lat1) = vn_tmp > 30 continue 902a1118,1138 > ! do mpi periodic points exchange for gswm with f2d(:) > ! Moved here from sub getgswm* because mpi necessary when f2d data > ! is allocated only for task-local subdomain block. > ! lons 1,2 <- nlonp4-3,nlonp4-2 and nlonp4-1,nlonp4 <- 3,4 > > #ifdef MPI > call mp_periodic_f2d(zsdi_gswm(lon0:lon1,lat0:lat1), > | lon0,lon1,lat0,lat1) > call mp_periodic_f2d(tnsdi_gswm(lon0:lon1,lat0:lat1), > | lon0,lon1,lat0,lat1) > call mp_periodic_f2d(unsdi_gswm(lon0:lon1,lat0:lat1), > | lon0,lon1,lat0,lat1) > call mp_periodic_f2d(vnsdi_gswm(lon0:lon1,lat0:lat1), > | lon0,lon1,lat0,lat1) > # else > set_periodic_f2d(zsdi_gswm) > set_periodic_f2d(tnsdi_gswm) > set_periodic_f2d(unsdi_gswm) > set_periodic_f2d(vnsdi_gswm) > #endif > ! 919a1156,1158 > #ifdef MPI > use mpi_module,only: mp_periodic_f2d > #endif 926c1165,1166 < | mon_cur,mon_nxt,ihr_cur,ihr_nxt,nointp=0 --- > | mon_cur,mon_nxt,ihr_cur,ihr_nxt,nointp=0, > | lonbeg,lonend 928,933c1168,1179 < | z_curmo(nlon,nlat)=0.,z_nxtmo(nlon,nlat)=0., ! UT interpolation < | tn_curmo(nlon,nlat)=0.,tn_nxtmo(nlon,nlat)=0., ! UT interpolation < | un_curmo(nlon,nlat)=0.,un_nxtmo(nlon,nlat)=0., ! UT interpolation < | vn_curmo(nlon,nlat)=0.,vn_nxtmo(nlon,nlat)=0., ! UT interpolation < | z_tmp(nlon,nlat)=0.,tn_tmp(nlon,nlat)=0., ! tmp < | un_tmp(nlon,nlat)=0.,vn_tmp(nlon,nlat)=0. ! tmp --- > | z_curmo(lon0:lon1,lat0:lat1), > | z_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation > | tn_curmo(lon0:lon1,lat0:lat1), > | tn_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation > | un_curmo(lon0:lon1,lat0:lat1), > | un_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation > | vn_curmo(lon0:lon1,lat0:lat1), > | vn_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation > | z_tmp(lon0:lon1,lat0:lat1), > | tn_tmp(lon0:lon1,lat0:lat1), ! tmp > | un_tmp(lon0:lon1,lat0:lat1), > | vn_tmp(lon0:lon1,lat0:lat1) ! tmp 984,991c1230,1245 < z_curmo(:,:) = gswm_nmidi_z(:,:,mon_cur,ihr_cur+1) ! current month < tn_curmo(:,:) = gswm_nmidi_tn(:,:,mon_cur,ihr_cur+1) < un_curmo(:,:) = gswm_nmidi_un(:,:,mon_cur,ihr_cur+1) < vn_curmo(:,:) = gswm_nmidi_vn(:,:,mon_cur,ihr_cur+1) < z_nxtmo(:,:) = gswm_nmidi_z(:,:,mon_nxt,ihr_cur+1) ! next month < tn_nxtmo(:,:) = gswm_nmidi_tn(:,:,mon_nxt,ihr_cur+1) < un_nxtmo(:,:) = gswm_nmidi_un(:,:,mon_nxt,ihr_cur+1) < vn_nxtmo(:,:) = gswm_nmidi_vn(:,:,mon_nxt,ihr_cur+1) --- > z_curmo(:,:) = gswm_nmidi_z (lon0:lon1,lat0:lat1, ! current month > | mon_cur,ihr_cur+1) > tn_curmo(:,:) = gswm_nmidi_tn(lon0:lon1,lat0:lat1, > | mon_cur,ihr_cur+1) > un_curmo(:,:) = gswm_nmidi_un(lon0:lon1,lat0:lat1, > | mon_cur,ihr_cur+1) > vn_curmo(:,:) = gswm_nmidi_vn(lon0:lon1,lat0:lat1, > | mon_cur,ihr_cur+1) > z_nxtmo(:,:) = gswm_nmidi_z (lon0:lon1,lat0:lat1, ! next month > | mon_nxt,ihr_cur+1) > tn_nxtmo(:,:) = gswm_nmidi_tn(lon0:lon1,lat0:lat1, > | mon_nxt,ihr_cur+1) > un_nxtmo(:,:) = gswm_nmidi_un(lon0:lon1,lat0:lat1, > | mon_nxt,ihr_cur+1) > vn_nxtmo(:,:) = gswm_nmidi_vn(lon0:lon1,lat0:lat1, > | mon_nxt,ihr_cur+1) 996,1019c1250,1273 < call timeliper_array(gswm_nmidi_z(:,:,mon_cur,ihr_cur+1), ! cur. month < | gswm_nmidi_z(:,:,mon_cur,ihr_nxt+1),difsec, < | secint,z_curmo,1) < call timeliper_array(gswm_nmidi_tn(:,:,mon_cur,ihr_cur+1), ! cur. month < | gswm_nmidi_tn(:,:,mon_cur,ihr_nxt+1),difsec, < | secint,tn_curmo,1) < call timeliper_array(gswm_nmidi_un(:,:,mon_cur,ihr_cur+1), ! cur. month < | gswm_nmidi_un(:,:,mon_cur,ihr_nxt+1),difsec, < | secint,un_curmo,1) < call timeliper_array(gswm_nmidi_vn(:,:,mon_cur,ihr_cur+1), ! cur. month < | gswm_nmidi_vn(:,:,mon_cur,ihr_nxt+1),difsec, < | secint,vn_curmo,1) < call timeliper_array(gswm_nmidi_z(:,:,mon_nxt,ihr_cur+1), ! next month < | gswm_nmidi_z(:,:,mon_nxt,ihr_nxt+1),difsec, < | secint,z_nxtmo,1) < call timeliper_array(gswm_nmidi_tn(:,:,mon_nxt,ihr_cur+1), ! next month < | gswm_nmidi_tn(:,:,mon_nxt,ihr_nxt+1),difsec, < | secint,tn_nxtmo,1) < call timeliper_array(gswm_nmidi_un(:,:,mon_nxt,ihr_cur+1), ! next month < | gswm_nmidi_un(:,:,mon_nxt,ihr_nxt+1),difsec, < | secint,un_nxtmo,1) < call timeliper_array(gswm_nmidi_vn(:,:,mon_nxt,ihr_cur+1), ! next month < | gswm_nmidi_vn(:,:,mon_nxt,ihr_nxt+1),difsec, < | secint,vn_nxtmo,1) --- > call timeliper_array(gswm_nmidi_z(lon0:lon1,lat0:lat1, ! cur. month > | mon_cur,ihr_cur+1),gswm_nmidi_z(lon0:lon1,lat0:lat1,mon_cur, > | ihr_nxt+1),difsec,secint,z_curmo,1) > call timeliper_array(gswm_nmidi_tn(lon0:lon1,lat0:lat1, ! cur. month > | mon_cur,ihr_cur+1),gswm_nmidi_tn(lon0:lon1,lat0:lat1,mon_cur, > | ihr_nxt+1),difsec,secint,tn_curmo,1) > call timeliper_array(gswm_nmidi_un(lon0:lon1,lat0:lat1, ! cur. month > | mon_cur,ihr_cur+1),gswm_nmidi_un(lon0:lon1,lat0:lat1,mon_cur, > | ihr_nxt+1),difsec,secint,un_curmo,1) > call timeliper_array(gswm_nmidi_vn(lon0:lon1,lat0:lat1, ! cur. month > | mon_cur,ihr_cur+1),gswm_nmidi_vn(lon0:lon1,lat0:lat1,mon_cur, > | ihr_nxt+1),difsec,secint,vn_curmo,1) > call timeliper_array(gswm_nmidi_z(lon0:lon1,lat0:lat1, ! next month > | mon_nxt,ihr_cur+1),gswm_nmidi_z(lon0:lon1,lat0:lat1,mon_nxt, > | ihr_nxt+1),difsec,secint,z_nxtmo,1) > call timeliper_array(gswm_nmidi_tn(lon0:lon1,lat0:lat1, ! next month > | mon_nxt,ihr_cur+1),gswm_nmidi_tn(lon0:lon1,lat0:lat1,mon_nxt, > | ihr_nxt+1),difsec,secint,tn_nxtmo,1) > call timeliper_array(gswm_nmidi_un(lon0:lon1,lat0:lat1, ! next month > | mon_nxt,ihr_cur+1),gswm_nmidi_un(lon0:lon1,lat0:lat1,mon_nxt, > | ihr_nxt+1),difsec,secint,un_nxtmo,1) > call timeliper_array(gswm_nmidi_vn(lon0:lon1,lat0:lat1, ! next month > | mon_nxt,ihr_cur+1),gswm_nmidi_vn(lon0:lon1,lat0:lat1,mon_nxt, > | ihr_nxt+1),difsec,secint,vn_nxtmo,1) 1060,1063c1314,1317 < znmidi_gswm(3:74,:) = z_curmo < tnnmidi_gswm(3:74,:) = tn_curmo < unnmidi_gswm(3:74,:) = un_curmo < vnnmidi_gswm(3:74,:) = vn_curmo --- > znmidi_gswm (lon0:lon1,lat0:lat1) = z_curmo > tnnmidi_gswm(lon0:lon1,lat0:lat1) = tn_curmo > unnmidi_gswm(lon0:lon1,lat0:lat1) = un_curmo > vnnmidi_gswm(lon0:lon1,lat0:lat1) = vn_curmo 1066,1069c1320,1323 < znmidi_gswm(3:74,:) = z_nxtmo < tnnmidi_gswm(3:74,:) = tn_nxtmo < unnmidi_gswm(3:74,:) = un_nxtmo < vnnmidi_gswm(3:74,:) = vn_nxtmo --- > znmidi_gswm (lon0:lon1,lat0:lat1) = z_nxtmo > tnnmidi_gswm(lon0:lon1,lat0:lat1) = tn_nxtmo > unnmidi_gswm(lon0:lon1,lat0:lat1) = un_nxtmo > vnnmidi_gswm(lon0:lon1,lat0:lat1) = vn_nxtmo 1073,1085c1327,1331 < znmidi_gswm(3:74,:) = z_tmp < tnnmidi_gswm(3:74,:) = tn_tmp < unnmidi_gswm(3:74,:) = un_tmp < vnnmidi_gswm(3:74,:) = vn_tmp < ! periodic points < 30 znmidi_gswm(1:2,:) = z_tmp(71:72,:) < tnnmidi_gswm(1:2,:) = tn_tmp(71:72,:) < unnmidi_gswm(1:2,:) = un_tmp(71:72,:) < vnnmidi_gswm(1:2,:) = vn_tmp(71:72,:) < znmidi_gswm(75:76,:) = z_tmp(1:2,:) < tnnmidi_gswm(75:76,:)= tn_tmp(1:2,:) < unnmidi_gswm(75:76,:)= un_tmp(1:2,:) < vnnmidi_gswm(75:76,:)= vn_tmp(1:2,:) --- > znmidi_gswm (lon0:lon1,lat0:lat1) = z_tmp > tnnmidi_gswm(lon0:lon1,lat0:lat1) = tn_tmp > unnmidi_gswm(lon0:lon1,lat0:lat1) = un_tmp > vnnmidi_gswm(lon0:lon1,lat0:lat1) = vn_tmp > 30 continue 1086a1333,1353 > ! do mpi periodic points exchange for gswm with f2d(:) > ! Moved here from sub getgswm* because mpi necessary when f2d data > ! is allocated only for task-local subdomain block. > ! lons 1,2 <- nlonp4-3,nlonp4-2 and nlonp4-1,nlonp4 <- 3,4 > > #ifdef MPI > call mp_periodic_f2d(znmidi_gswm(lon0:lon1,lat0:lat1), > | lon0,lon1,lat0,lat1) > call mp_periodic_f2d(tnnmidi_gswm(lon0:lon1,lat0:lat1), > | lon0,lon1,lat0,lat1) > call mp_periodic_f2d(unnmidi_gswm(lon0:lon1,lat0:lat1), > | lon0,lon1,lat0,lat1) > call mp_periodic_f2d(vnnmidi_gswm(lon0:lon1,lat0:lat1), > | lon0,lon1,lat0,lat1) > # else > set_periodic_f2d(znmidi_gswm) > set_periodic_f2d(tnnmidi_gswm) > set_periodic_f2d(unnmidi_gswm) > set_periodic_f2d(vnnmidi_gswm) > #endif > ! 1104a1372,1374 > #ifdef MPI > use mpi_module,only: mp_periodic_f2d > #endif 1111c1381,1382 < | mon_cur,mon_nxt,ihr_cur,ihr_nxt,nointp=0 --- > | mon_cur,mon_nxt,ihr_cur,ihr_nxt,nointp=0, > | lonbeg,lonend 1113,1118c1384,1395 < | z_curmo(nlon,nlat)=0.,z_nxtmo(nlon,nlat)=0., ! UT interpolation < | tn_curmo(nlon,nlat)=0.,tn_nxtmo(nlon,nlat)=0., ! UT interpolation < | un_curmo(nlon,nlat)=0.,un_nxtmo(nlon,nlat)=0., ! UT interpolation < | vn_curmo(nlon,nlat)=0.,vn_nxtmo(nlon,nlat)=0., ! UT interpolation < | z_tmp(nlon,nlat)=0.,tn_tmp(nlon,nlat)=0., ! tmp < | un_tmp(nlon,nlat)=0.,vn_tmp(nlon,nlat)=0. ! tmp --- > | z_curmo(lon0:lon1,lat0:lat1), > | z_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation > | tn_curmo(lon0:lon1,lat0:lat1), > | tn_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation > | un_curmo(lon0:lon1,lat0:lat1), > | un_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation > | vn_curmo(lon0:lon1,lat0:lat1), > | vn_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation > | z_tmp(lon0:lon1,lat0:lat1), > | tn_tmp(lon0:lon1,lat0:lat1), ! tmp > | un_tmp(lon0:lon1,lat0:lat1), > | vn_tmp(lon0:lon1,lat0:lat1) ! tmp 1169,1176c1446,1461 < z_curmo(:,:) = gswm_nmisdi_z(:,:,mon_cur,ihr_cur+1) ! current month < tn_curmo(:,:) = gswm_nmisdi_tn(:,:,mon_cur,ihr_cur+1) < un_curmo(:,:) = gswm_nmisdi_un(:,:,mon_cur,ihr_cur+1) < vn_curmo(:,:) = gswm_nmisdi_vn(:,:,mon_cur,ihr_cur+1) < z_nxtmo(:,:) = gswm_nmisdi_z(:,:,mon_nxt,ihr_cur+1) ! next month < tn_nxtmo(:,:) = gswm_nmisdi_tn(:,:,mon_nxt,ihr_cur+1) < un_nxtmo(:,:) = gswm_nmisdi_un(:,:,mon_nxt,ihr_cur+1) < vn_nxtmo(:,:) = gswm_nmisdi_vn(:,:,mon_nxt,ihr_cur+1) --- > z_curmo(:,:) = gswm_nmisdi_z (lon0:lon1,lat0:lat1, ! current month > | mon_cur,ihr_cur+1) > tn_curmo(:,:) = gswm_nmisdi_tn(lon0:lon1,lat0:lat1, > | mon_cur,ihr_cur+1) > un_curmo(:,:) = gswm_nmisdi_un(lon0:lon1,lat0:lat1, > | mon_cur,ihr_cur+1) > vn_curmo(:,:) = gswm_nmisdi_vn(lon0:lon1,lat0:lat1, > | mon_cur,ihr_cur+1) > z_nxtmo(:,:) = gswm_nmisdi_z (lon0:lon1,lat0:lat1, ! next month > | mon_nxt,ihr_cur+1) > tn_nxtmo(:,:) = gswm_nmisdi_tn(lon0:lon1,lat0:lat1, > | mon_nxt,ihr_cur+1) > un_nxtmo(:,:) = gswm_nmisdi_un(lon0:lon1,lat0:lat1, > | mon_nxt,ihr_cur+1) > vn_nxtmo(:,:) = gswm_nmisdi_vn(lon0:lon1,lat0:lat1, > | mon_nxt,ihr_cur+1) 1181,1204c1466,1489 < call timeliper_array(gswm_nmisdi_z(:,:,mon_cur,ihr_cur+1), ! cur. month < | gswm_nmisdi_z(:,:,mon_cur,ihr_nxt+1),difsec, < | secint,z_curmo,1) < call timeliper_array(gswm_nmisdi_tn(:,:,mon_cur,ihr_cur+1), ! cur. month < | gswm_nmisdi_tn(:,:,mon_cur,ihr_nxt+1),difsec, < | secint,tn_curmo,1) < call timeliper_array(gswm_nmisdi_un(:,:,mon_cur,ihr_cur+1), ! cur. month < | gswm_nmisdi_un(:,:,mon_cur,ihr_nxt+1),difsec, < | secint,un_curmo,1) < call timeliper_array(gswm_nmisdi_vn(:,:,mon_cur,ihr_cur+1), ! cur. month < | gswm_nmisdi_vn(:,:,mon_cur,ihr_nxt+1),difsec, < | secint,vn_curmo,1) < call timeliper_array(gswm_nmisdi_z(:,:,mon_nxt,ihr_cur+1), ! next month < | gswm_nmisdi_z(:,:,mon_nxt,ihr_nxt+1),difsec, < | secint,z_nxtmo,1) < call timeliper_array(gswm_nmisdi_tn(:,:,mon_nxt,ihr_cur+1), ! next month < | gswm_nmisdi_tn(:,:,mon_nxt,ihr_nxt+1),difsec, < | secint,tn_nxtmo,1) < call timeliper_array(gswm_nmisdi_un(:,:,mon_nxt,ihr_cur+1), ! next month < | gswm_nmisdi_un(:,:,mon_nxt,ihr_nxt+1),difsec, < | secint,un_nxtmo,1) < call timeliper_array(gswm_nmisdi_vn(:,:,mon_nxt,ihr_cur+1), ! next month < | gswm_nmisdi_vn(:,:,mon_nxt,ihr_nxt+1),difsec, < | secint,vn_nxtmo,1) --- > call timeliper_array(gswm_nmisdi_z(lon0:lon1,lat0:lat1, ! cur. month > | mon_cur,ihr_cur+1),gswm_nmisdi_z(lon0:lon1,lat0:lat1,mon_cur, > | ihr_nxt+1),difsec,secint,z_curmo,1) > call timeliper_array(gswm_nmisdi_tn(lon0:lon1,lat0:lat1, ! cur. month > | mon_cur,ihr_cur+1),gswm_nmisdi_tn(lon0:lon1,lat0:lat1,mon_cur, > | ihr_nxt+1),difsec,secint,tn_curmo,1) > call timeliper_array(gswm_nmisdi_un(lon0:lon1,lat0:lat1, ! cur. month > | mon_cur,ihr_cur+1),gswm_nmisdi_un(lon0:lon1,lat0:lat1,mon_cur, > | ihr_nxt+1),difsec,secint,un_curmo,1) > call timeliper_array(gswm_nmisdi_vn(lon0:lon1,lat0:lat1, ! cur. month > | mon_cur,ihr_cur+1),gswm_nmisdi_vn(lon0:lon1,lat0:lat1,mon_cur, > | ihr_nxt+1),difsec,secint,vn_curmo,1) > call timeliper_array(gswm_nmisdi_z(lon0:lon1,lat0:lat1, ! next month > | mon_nxt,ihr_cur+1),gswm_nmisdi_z(lon0:lon1,lat0:lat1,mon_nxt, > | ihr_nxt+1),difsec,secint,z_nxtmo,1) > call timeliper_array(gswm_nmisdi_tn(lon0:lon1,lat0:lat1, ! next month > | mon_nxt,ihr_cur+1),gswm_nmisdi_tn(lon0:lon1,lat0:lat1,mon_nxt, > | ihr_nxt+1),difsec,secint,tn_nxtmo,1) > call timeliper_array(gswm_nmisdi_un(lon0:lon1,lat0:lat1, ! next month > | mon_nxt,ihr_cur+1),gswm_nmisdi_un(lon0:lon1,lat0:lat1,mon_nxt, > | ihr_nxt+1),difsec,secint,un_nxtmo,1) > call timeliper_array(gswm_nmisdi_vn(lon0:lon1,lat0:lat1, ! next month > | mon_nxt,ihr_cur+1),gswm_nmisdi_vn(lon0:lon1,lat0:lat1,mon_nxt, > | ihr_nxt+1),difsec,secint,vn_nxtmo,1) 1244,1247c1529,1532 < znmisdi_gswm(3:74,:) = z_curmo < tnnmisdi_gswm(3:74,:) = tn_curmo < unnmisdi_gswm(3:74,:) = un_curmo < vnnmisdi_gswm(3:74,:) = vn_curmo --- > znmisdi_gswm (lon0:lon1,lat0:lat1) = z_curmo > tnnmisdi_gswm(lon0:lon1,lat0:lat1) = tn_curmo > unnmisdi_gswm(lon0:lon1,lat0:lat1) = un_curmo > vnnmisdi_gswm(lon0:lon1,lat0:lat1) = vn_curmo 1250,1253c1535,1538 < znmisdi_gswm(3:74,:) = z_nxtmo < tnnmisdi_gswm(3:74,:) = tn_nxtmo < unnmisdi_gswm(3:74,:) = un_nxtmo < vnnmisdi_gswm(3:74,:) = vn_nxtmo --- > znmisdi_gswm (lon0:lon1,lat0:lat1) = z_nxtmo > tnnmisdi_gswm(lon0:lon1,lat0:lat1) = tn_nxtmo > unnmisdi_gswm(lon0:lon1,lat0:lat1) = un_nxtmo > vnnmisdi_gswm(lon0:lon1,lat0:lat1) = vn_nxtmo 1257,1269c1542,1546 < znmisdi_gswm(3:74,:) = z_tmp < tnnmisdi_gswm(3:74,:) = tn_tmp < unnmisdi_gswm(3:74,:) = un_tmp < vnnmisdi_gswm(3:74,:) = vn_tmp < ! periodic points < 30 znmisdi_gswm(1:2,:) = z_tmp(71:72,:) < tnnmisdi_gswm(1:2,:) = tn_tmp(71:72,:) < unnmisdi_gswm(1:2,:) = un_tmp(71:72,:) < vnnmisdi_gswm(1:2,:) = vn_tmp(71:72,:) < znmisdi_gswm(75:76,:) = z_tmp(1:2,:) < tnnmisdi_gswm(75:76,:)= tn_tmp(1:2,:) < unnmisdi_gswm(75:76,:)= un_tmp(1:2,:) < vnnmisdi_gswm(75:76,:)= vn_tmp(1:2,:) --- > znmisdi_gswm (lon0:lon1,lat0:lat1) = z_tmp > tnnmisdi_gswm(lon0:lon1,lat0:lat1) = tn_tmp > unnmisdi_gswm(lon0:lon1,lat0:lat1) = un_tmp > vnnmisdi_gswm(lon0:lon1,lat0:lat1) = vn_tmp > 30 continue 1270a1548,1568 > ! do mpi periodic points exchange for gswm with f2d(:) > ! Moved here from sub getgswm* because mpi necessary when f2d data > ! is allocated only for task-local subdomain block. > ! lons 1,2 <- nlonp4-3,nlonp4-2 and nlonp4-1,nlonp4 <- 3,4 > > #ifdef MPI > call mp_periodic_f2d(znmisdi_gswm(lon0:lon1,lat0:lat1), > | lon0,lon1,lat0,lat1) > call mp_periodic_f2d(tnnmisdi_gswm(lon0:lon1,lat0:lat1), > | lon0,lon1,lat0,lat1) > call mp_periodic_f2d(unnmisdi_gswm(lon0:lon1,lat0:lat1), > | lon0,lon1,lat0,lat1) > call mp_periodic_f2d(vnnmisdi_gswm(lon0:lon1,lat0:lat1), > | lon0,lon1,lat0,lat1) > # else > set_periodic_f2d(znmisdi_gswm) > set_periodic_f2d(tnnmisdi_gswm) > set_periodic_f2d(unnmisdi_gswm) > set_periodic_f2d(vnnmisdi_gswm) > #endif > ! 1291,1292c1589,1591 < real,intent(in) :: difd1d2,difint,d1(nlon,nlat),d2(nlon,nlat) < real,intent(out) :: fout(nlon,nlat) --- > real,intent(in) :: difd1d2,difint,d1(lon0:lon1,lat0:lat1), > | d2(lon0:lon1,lat0:lat1) > real,intent(out) :: fout(lon0:lon1,lat0:lat1) 1310,1311c1609,1610 < do i = 1,nlon < do j = 1,nlat --- > do i = lon0,lon1 > do j = lat0,lat1 1316a1616,1625 > !----------------------------------------------------------------------- > subroutine set_periodic_f2d(f) > ! > ! Set periodic points for all f2d fields (serial or non-mpi only): > ! > real,intent(inout) :: f(nlonp4,nlat) > > ! lons 1,2 <- nlonp4-3,nlonp4-2 and nlonp4-1,nlonp4 <- 3,4 > f(1:2,:) = f(nlonp4-3:nlonp4-2,:) > f(nlonp4-1:nlonp4,:) = f(3:4,:) 1317a1627,1628 > end subroutine set_periodic_f2d > ======================================================================== Diff of /home/tgcm/tiegcm1/heelis.F and /home/foster/tiegcm/amie/modsrc_oct03/heelis.F 2,3c2 < subroutine heelis(colatc,pfrac,phihm,nmlat0,nmlon0,nmlat,nmlon, < | nmlonp1) --- > module heelis_module 4a4,13 > ! Module used to calculcate the Heelis model potential in both hemispheres > ! Byimf, Ctpoten and Power at a minimum using paramaters from aurora_cons > ! > use params_module,only: nmlat,nmlonp1,nmlon,nmlonp1 > use dynamo_module,only: nmlat0,pfrac,phihm > implicit none > contains > !----------------------------------------------------------------------- > subroutine heelis > ! 6c15 < ! These routines return colatc, pfrac, and phihm to the dynamo. --- > ! These routines return pfrac and phihm to the dynamo. 8,9d16 < ! This is not a module because of difficulties with circular dependencies < ! with the dynamo. Instead, fields are passed through argument lists. 14d20 < ! colatc: Colatitudes of northern geomagnetic grid points: 16c22 < ! colatitudes crit(2). --- > ! convection colatitudes crit(2). 19,22c25,26 < integer,intent(in) :: nmlat0,nmlon0,nmlon,nmlat,nmlonp1 < real,dimension(nmlonp1,nmlat0),intent(out) :: colatc < real,dimension(nmlon0,nmlat0),intent(out) :: pfrac < real,dimension(nmlonp1,nmlat),intent(out) :: phihm --- > ! Calculate pfrac fractional presence of dynamo equation using critical > ! convection colatitudes crit(2). (crit is in cons module) 24,27c28 < ! Define colatc with northern auroral colatitudes corresponding to < ! each northern geomagnetic grid point. Also calculate fractional < ! presence of dynamo equation given critical colatitudes crit(2). < ! (crit is in cons module) --- > call colath 29,30d29 < call colath(colatc,pfrac,nmlat0,nmlon0,nmlonp1) < ! 34c33 < call potm(phihm,nmlon,nmlonp1,nmlat) --- > call potm 38c37 < subroutine colath(colatc,pfrac,nmlat0,nmlon0,nmlonp1) --- > subroutine colath 40,42c39,40 < ! Define colatc with northern auroral circle coordinates, corresponding < ! to each northern geomagnetic grid point. Also calculate fractional < ! presence of dynamo equation given critical colatitudes crit(2). --- > ! Calculate pfrac fractional presence of dynamo equation using critical > ! convection colatitudes crit(2). (crit is in cons module) 52,54c50,53 < integer,intent(in) :: nmlat0,nmlon0,nmlonp1 < real,intent(out) :: colatc(nmlonp1,nmlat0) < real,intent(out) :: pfrac(nmlon0,nmlat0) --- > ! integer,intent(in) :: nmlat0,nmlonp1 > ! real,intent(out) :: pfrac(nmlonp1,nmlat0) > > real,dimension(nmlonp1,nmlat0) :: colatc 64a64,67 > ! TEMP > ! write (6,"(1x,'COLATH: crit1,2 dskofc offc deg=',6e12.4)") > ! | crit(1)*rtd,crit(2)*rtd,dskofc(1)*rtd,offc(1)*rtd, > ! | dskofc(2)*rtd,offc(2)*rtd 66c69 < ! Define colatc with colatitude values: --- > ! Define colatc with northern convection circle coordinates 78c81 < ! | j,colatc(:,j) --- > ! | j,colatc(:,j)*rtd 82c85 < ! hemisphere geomagnetic grid point. Output in pfrac(nmlon0,nmlat0) --- > ! hemisphere geomagnetic grid point. Output in pfrac(nmlonp1,nmlat0) 84c87 < do i=1,nmlon0 --- > do i=1,nmlonp1 88c91 < enddo ! i=1,nmlon0 --- > enddo ! i=1,nmlonp1 97c100 < subroutine potm(phihm,nmlon,nmlonp1,nmlat) --- > subroutine potm 101a105 > use amie_module,only: tiepot ! (nmlonp1,nmlat) 106,109d109 < ! Args: < integer,intent(in) :: nmlon,nmlonp1,nmlat < real,intent(out) :: phihm(nmlonp1,nmlat) < ! 304a305,306 > !----------------------------------------------------------------------- > end module heelis_module ======================================================================== Diff of /home/tgcm/tiegcm1/init.F and /home/foster/tiegcm/amie/modsrc_oct03/init.F 77c77,78 < | mxhist_prim,mxhist_sech,output,secout,mkhvols,nmc,source_start --- > | mxhist_prim,mxhist_sech,output,secout,mkhvols,nmc,source_start, > | amienh,amiesh,aurora_proton 86a88,89 > use amie_module,only: rdamie_nh,rdamie_sh > use low_proton_module,only: rdproton,getproton 189c192 < ! Initialize amie, and get amie file if necessary: --- > ! Read amie data files if necessary: 190a194,205 > if (len_trim(aurora_proton) > 0) then > write(6,"('Reading AURORA PROTON file ',a)") trim(aurora_proton) > call rdproton > endif > if (len_trim(amienh) > 0) then > write(6,"('Reading AMIENH file ',a)") trim(amienh) > call rdamie_nh > endif > if (len_trim(amiesh) > 0) then > write(6,"('Reading AMIESH file ',a)") trim(amiesh) > call rdamie_sh > endif ======================================================================== Diff of /home/tgcm/tiegcm1/input.F and /home/foster/tiegcm/amie/modsrc_oct03/input.F 5c5 < | spval,ispval,mxday,nlonp4 --- > | spval,ispval,mxday,nlonp4,mxind_time 42c42,44 < | amievol ! file or mss path of amie data file (optional) --- > | aurora_proton, ! file or mss path of aurora proton file (optional) > | amiesh, ! file or mss path of amie SH data file (optional) > | amienh ! file or mss path of amie NH data file (optional) 61a64,66 > | ilow_proton, ! 0/1 flag for auroral proton data (not in namelist read) > | iamie, ! 0/1 flag for AMIE data (not in namelist read) > | iamie_bkg, ! 0/1/2 flag for read real, 1st, or 24-hr averaged data 74a80,81 > | power_time(4,mxind_time), ! time-dependent hemispheric power input > | ctpoten_time(4,mxind_time), ! time-dependent cross tail potential input 75a83,87 > | bzimf, ! Bz component of IMF in nT > | swvel, ! Solar wind velocity in km/s > | swden, ! Solar wind density in #/cm3 > | AL, ! AL lower magnetic auroral activity index in nT > ! if present, ALUSE=true; if absent, AL=-20, ALUSE=false 76a89 > logical :: aluse ! logical to use AL in Weimer model or not 77a91,101 > | potential_model, ! electric potential model used > ! Values can be 'HEELIS', 'WEIMER', or 'NONE' > ! If absent, the default value is set to 'HEELIS' > | weimer_ncfile, ! mss path or file path to netcdf weimer coef file > ! see comments in weimer_mod.f > | weimer_ncfiledef, ! default path to weimer coef data file > ! hrindices_ncfile NOT MADE YET! 11/15/02: > | hrindices_ncfile, ! mss path or file path to netcdf hourly geophysical > ! data indices file. see comments in hrindices_mod.f > | hrindices_ncfiledef, ! default path to hrly indices data file > ! 94a119,121 > integer :: > | ntimes_ctpoten, ! number of times provided in ctpoten_time > | ntimes_power ! number of times provided in power_time 133,135c160,163 < | label, tempdir, magvol, amievol, gpi_ncfile, < | gswm_di_ncfile, gswm_sdi_ncfile, gswm_nmidi_ncfile, < | gswm_nmisdi_ncfile --- > | label, tempdir, magvol, amiesh,amienh, gpi_ncfile, > | gswm_di_ncfile, gswm_sdi_ncfile, gswm_nmidi_ncfile, > | gswm_nmisdi_ncfile, potential_model, weimer_ncfile, > | hrindices_ncfile,aurora_proton 137c165,166 < | step,start_day,start_year,calendar_advance, --- > | step,start_day,start_year,calendar_advance, ilow_proton, > | iamie, iamie_bkg, 142c171,172 < | f107,f107a, power,ctpoten,byimf,colfac --- > | f107,f107a, power,ctpoten,byimf,colfac,bzimf,swvel,swden,al, > | power_time(4,mxind_time),ctpoten_time(4,mxind_time) 179c209 < | label,tempdir,magvol,amievol,date,calday,step,dispose, --- > | label,tempdir,magvol,date,calday,step,dispose, 182c212 < | secfmag,secfgeo2d,secfmag2d,secfmagphr, --- > | secfmag,secfgeo2d,secfmag2d,secfmagphr,potential_model, 184,187c214,219 < | power,ctpoten,byimf,colfac,tideann,aurora,magphr,gpi_ncfile, < | gswm_di_ncfile,gswm_sdi_ncfile,gswm_nmidi_ncfile, < | gswm_nmisdi_ncfile,mxhist_prim,mxhist_sech,msreten,ntask_lat, < | ntask_lon,start_day,start_year,calendar_advance --- > | power,ctpoten,byimf,bzimf,swvel,swden,al,colfac,tideann,aurora, > | magphr,gpi_ncfile,gswm_di_ncfile,gswm_sdi_ncfile, > | gswm_nmidi_ncfile,gswm_nmisdi_ncfile,mxhist_prim,mxhist_sech, > | msreten,ntask_lat,ntask_lon,start_day,start_year, > | calendar_advance,ctpoten_time,power_time,amiesh,amienh, > | iamie_bkg,aurora_proton,ilow_proton 234c266,268 < amievol = ' ' --- > amiesh = ' ' > amienh = ' ' > aurora_proton = ' ' 239c273,275 < write(gpi_ncfiledef,"('/TGCM/data/gpi_1979001-2001031.nc')") --- > ! write(gpi_ncfiledef,"('/TGCM/data/gpi_1979001-2001031.nc')") > ! write(gpi_ncfiledef,"('/TGCM/data/gpi_1979001-2002181.nc')") > write(gpi_ncfiledef,"('/TGCM/data/gpi_1979001-2002365.nc')") 241a278,286 > weimer_ncfile = ' ' ! user input weimer coef file path > weimer_ncfiledef = ' ' ! default weimer coef file path > write(weimer_ncfiledef,"('/TGCM/data/weimer2001_coeffs.nc')") > ! > hrindices_ncfile = ' ' ! user input hrly indices data file path > hrindices_ncfiledef = ' ' ! default hrly indices data file path > ! write(hrindices_ncfiledef, > ! | "('/TGCM/data/hrindices_1988001-2002181.nc')") > ! 265a311,313 > ilow_proton = 0 ! will be set 1 if the auroral proton file is given > iamie = 0 ! will be set 1 if amienh or amiesh files are given > iamie_bkg = 0 ! default is read in real time amie data 272a321,322 > power_time(:,:) = spval > ctpoten_time(:,:) = spval 273a324,327 > bzimf = spval > swvel = spval > swden = spval > al = spval 275a330 > potential_model = ' ' 308a364 > real :: rday,rhour,rmin,rval 544,545d599 < ! c(62) is replaced by f107a from input_mod. < ! c(62) = f107a ! this may get reset later by getgpi or tail 547c601,620 < ! Cross-tail potential: --- > ! There are 3 input options for ctpoten and power: > ! 1) Constants are provided by the user (ctpoten, power) > ! 2) Time series are provided by the user (ctpoten_time, power_time) > ! 3) GPI database (neither constans nor time series were provided) > ! > ! Cannot provide both constant and time series: > if (ctpoten /= spval .and. any(ctpoten_time /= spval)) then > write(6,"(/,'>>> INPUT: you cannot provide both a constant', > | '(CTPOTEN) and a time-series (CTPOTEN_TIME) for cross-', > | 'tail potential.')") > stop 'CTPOTEN' > endif > if (power /= spval .and. any(power_time /= spval)) then > write(6,"(/,'>>> INPUT: you cannot provide both a constant', > | ' (POWER) and a time-series (POWER_TIME) for hemispheric', > | ' power.')") > stop 'POWER' > endif > ! > ! Validate cross-tail potential constant, if provided: 554c627 < ! Hemispheric power: --- > ! Validate hemispheric power constant, if provided: 560d632 < ! hpower = power ! hpower is in /ingpi/ 561a634,640 > ! Validate times and values in user-provided time series for > ! crosstail potential and hemispheric power: > call validate_indices_time(ctpoten_time,ntimes_ctpoten, > | 'CTPOTEN_TIME') > call validate_indices_time(power_time,ntimes_power, > | 'POWER_TIME') > ! 571a651,652 > if (any(power_time /= spval)) gpi_vars(3) = 0. ! dummy non-spval > if (any(ctpoten_time /= spval)) gpi_vars(4) = 0. ! dummy non-spval 619c700,705 < write(6,"(' ',a,' = ',f10.2)") gpi_names(i),gpi_vars(i) --- > if ((trim(gpi_names(i))=="ctpoten".and.ntimes_ctpoten > 0).or. > | (trim(gpi_names(i))=="ctpoten".and.ntimes_ctpoten > 0)) then > write(6,"(' ',a,' = (time series)')") gpi_names(i) > else > write(6,"(' ',a,' = ',f10.2)") gpi_names(i),gpi_vars(i) > endif 634a721,748 > ! Check to see if have hrindices_ncfile (when avail! 11/15/02) > ! > ! Check electric potential model: > if (potential_model == 'WEIMER') then > write (6,"(4x,'Will use the Weimer 2001 electric potential', > | ' model')") > if (len_trim(weimer_ncfile)==0) then > write(6,"(4x,'Since Weimer coef data file WEIMER_NCFILE was ', > | 'NOT provided,',/,4x,'I will use the default ', > | 'WEIMER_NCFILE: ',a)") trim(weimer_ncfiledef) > weimer_ncfile = weimer_ncfiledef > else > write(6,"(4x,'Will use user-provided path to Weimer coefs', > | ' Weimer_NCFILE = ',a)") trim(weimer_ncfile) > endif > elseif (potential_model == 'HEELIS') then > write (6,"(4x,'Will use the Heelis param electric potential', > | ' model')") > elseif (potential_model == 'NONE') then > write (6,"(4x,'Will use NONE (ie zero) electric potential', > | ' model')") > endif > if (len_trim(potential_model)==0) then > write (6,"(4x,'Will use default Heelis elecric potential', > | ' model')") > potential_model='HEELIS' > endif > ! 643,646c757,769 < if (byimf > 7.) then < write(6,"(/'>>> INPUT: byimf too big (must be <= 7.):', < | ' byimf=',e12.4)") byimf < stop 'BYIMF' --- > ! > ! Limits on By NOT necessary for Weimer (or for NONE) > if (potential_model == 'HEELIS') then > if (byimf > 7.) then > write(6,"(/'>>> INPUT: byimf too big (must be <= 7.):', > | ' byimf=',e12.4)") byimf > stop 'BYIMF' > endif > if (byimf < -11.) then > write(6,"(/'>>> INPUT: byimf too small (must be >= -11.):', > | ' byimf=',e12.4)") byimf > stop 'BYIMF' > endif 648,651c771,778 < if (byimf < -11.) then < write(6,"(/'>>> INPUT: byimf too small (must be >= -11.):', < | ' byimf=',e12.4)") byimf < stop 'BYIMF' --- > ! > ! Bz component of solar IMF magnetic field: > if (bzimf==spval) then > write(6,"(/'>>> INPUT: need 1 value for bzimf', > | ' e.g.: BZIMF=0.')") > write(6,"('NOTE: multiple (time-dependent) values for BZIMF', > | ' not yet supported')") > stop 'BZIMF' 653a781,808 > ! Solar wind velocity: > if (swvel==spval) then > write(6,"(/'>>> INPUT: Use default value for swvel=400 km/s')") > write(6,"('NOTE: multiple (time-dependent) values for SWVEL', > | ' not yet supported')") > swvel = 400. > endif > ! > ! Solar wind density: > if (swden==spval) then > write(6,"(/'>>> INPUT: Use default value for swden=4 #/cm3')") > write(6,"('NOTE: multiple (time-dependent) values for SWDEN', > | ' not yet supported')") > swden = 4. > endif > ! > ! AL, lower auroral magnetic activity index: > if (AL==spval) then > write(6,"(/'>>> INPUT: AL missing, set AL=-20. ALUSE=FALSE')") > write(6,"('NOTE: multiple (time-dependent) values for AL', > | ' not yet supported')") > AL = -20. > ALUSE = .FALSE. > else > write(6,"(/'Input: AL given, setting ALUSE=TRUE')") > ALUSE = .TRUE. > endif > ! 692a848,863 > ! > ! Set ilow_proton flag to 1 if user has provided the FUV proton data: > ! > if (len_trim(aurora_proton) > 0) then > write(6,"('Will read auroral proton data inputs')") > ilow_proton = 1 > endif > ! > ! Set iamie flag to 1 if user has provided either amiesh or amienh > ! data files: > ! > if (len_trim(amiesh) + len_trim(amienh) > 0) then > write(6,"('Will read AMIE data inputs')") > iamie = 1 > potential_model = 'AMIE' > endif 760c931,932 < | nsec_hist (mxseries_sech),nsec_save(mxseries_sech) --- > | nsec_hist (mxseries_sech),nsec_save(mxseries_sech), > | nsec_index,nsec0 776,777c948,949 < elseif (dispose /= 0 .and. dispose /= 1) then < write(6,"(/,'>>> INPUT: DISPOSE flag must be 0 or 1:', --- > elseif (dispose /= 0 .and. dispose /= 1 .and. dispose /= 2) then > write(6,"(/,'>>> INPUT: DISPOSE flag must be 0, 1, or 2:', 778a951,955 > write(6,"('DISPOSE = 0 -> do not dispose to mss.')") > write(6,"('DISPOSE = 1 -> dispose to mss during model', > | ' execution.')") > write(6,"('DISPOSE = 2 -> dispose to mss after model', > | ' execution.')") 1067a1245,1299 > ! Validate that times given by ctpoten_time and/or power_time are > ! within start and stop times (the times themselves were validated > ! in inp_model). Note user must provide a value for each start time. > ! > if (ntimes_ctpoten > 0) then > ! > ! Must provide first ctpoten_time(:,1) same as START time. > nsec0 = mtime_to_nsec(int(ctpoten_time(1:3,1))) > if (nsec0 /= nsec_start(1)) then > write(6,"(/,'>>> INPUT: Please provide first CTPOTEN_TIME', > | ' time equal to first START time.')") > stop 'starting ctpoten_time' > endif > ! > ! All ctpoten times must be within start and stop times: > do i=1,ntimes_ctpoten > nsec_index = mtime_to_nsec(int(ctpoten_time(1:3,i))) > if (nsec_index < nsec_start(1) .or. > | nsec_index > nsec_stop(1)) then > write(6,"(/,'>>> INPUT: CTPOTEN_TIME ',3i4,' is outside ', > | 'model START/STOP times: START=',3i4,' STOP=',3i4)") > | int(ctpoten_time(1:3,i)),start(:,1),stop(:,1) > stop 'ctpoten_time' > endif > enddo > ! > ! Init: > ctpoten = ctpoten_time(4,1) > endif ! ntimes_ctpoten > if (ntimes_power > 0) then > ! > ! Must provide first power_time(:,1) same as START time. > nsec0 = mtime_to_nsec(int(power_time(1:3,1))) > if (nsec0 /= nsec_start(1)) then > write(6,"(/,'>>> INPUT: Please provide first POWER_TIME', > | ' time equal to first START time.')") > stop 'starting power_time' > endif > ! > ! All power times must be within start and stop times: > do i=1,ntimes_power > nsec_index = mtime_to_nsec(int(power_time(1:3,i))) > if (nsec_index < nsec_start(1) .or. > | nsec_index > nsec_stop(1)) then > write(6,"(/,'>>> INPUT: POWER_TIME ',3i4,' is outside ', > | 'model START/STOP times: START=',3i4,' STOP=',3i4)") > | int(power_time(1:3,i)),start(:,1),stop(:,1) > stop 'power_time' > endif > enddo > ! > ! Init: > power = power_time(4,1) > endif ! ntimes_power > ! 1656c1888,1893 < inp%amievol = amievol --- > inp%amiesh = aurora_proton > inp%amiesh = amiesh > inp%amienh = amienh > inp%potential_model = potential_model > inp%hrindices_ncfile = hrindices_ncfile > inp%weimer_ncfile = weimer_ncfile 1678a1916,1918 > inp%ilow_proton = ilow_proton > inp%iamie = iamie > inp%iamie_bkg = iamie_bkg 1682a1923,1924 > inp%power_time(:,:) = power_time(:,:) > inp%ctpoten_time(:,:) = ctpoten_time(:,:) 1683a1926,1929 > inp%bzimf = bzimf > inp%swvel = swvel > inp%swden = swden > inp%AL = AL 1732a1979,1995 > if (len_trim(potential_model) > 0) > | write(6,"(' high-lat electric potential model: ', > | 'potential_model = ',a)") trim(potential_model) > if (len_trim(hrindices_ncfile) > 0) > | write(6,"(' hrly indices: hrindices_ncfile = ',a)") > | trim(hrindices_ncfile) > if (len_trim(weimer_ncfile) > 0) > | write(6,"(' weimer coefs: weimer_ncfile = ',a)") > | trim(weimer_ncfile) > if (len_trim(amiesh) > 0) > | write(6,"(' amiesh = ',a,/,4x, > | '(file or mss path containing amie SH data)')") > | trim(inp%amiesh) > if (len_trim(amienh) > 0) > | write(6,"(' amienh = ',a,/,4x, > | '(file or mss path containing amie NH data)')") > | trim(inp%amienh) 1747,1749d2009 < if (len_trim(amievol) > 0) < | write(6,"(' amievol = ',a,/,4x, < | '(file or mss path containing amie data)')") trim(inp%amievol) 1765c2025 < write(6,"(' dispose =',i4,' (mss dispose flag 0/1)')") --- > write(6,"(' dispose =',i4,' (mss dispose flag 0/1/2)')") 1915,1916c2175,2181 < write(6,"(' power = ',f9.3,' (hemispheric power (gw))')") < | inp%power --- > if (ntimes_power==0) then > write(6,"(' power = ',f9.3,' (hemispheric power (gw)')") > | inp%power > else > write(6,"(' power = ',f9.3,' (hemispheric power (gw)', > | ' (user provided time series)')") inp%power > endif 1923,1924c2188,2194 < write(6,"(' ctpoten = ',f9.3,' (cross-cap potential ', < | '(volts))')") inp%ctpoten --- > if (ntimes_ctpoten==0) then > write(6,"(' ctpoten = ',f9.3,' (cross-cap potential ', > | '(volts)')") inp%ctpoten > else > write(6,"(' ctpoten = ',f9.3,' (cross-cap potential ', > | '(volts) (user provided time series)')") inp%ctpoten > endif 1932c2202,2209 < --- > write(6,"(' bzimf = ',f9.3,' (Bz component of IMF)')") > | inp%bzimf > write(6,"(' swvel = ',f9.3,' (solar wind velocity)')") > | inp%swvel > write(6,"(' swden = ',f9.3,' (solar wind density)')") > | inp%swden > write(6,"(' AL = ',f9.3,' (AL, lower auroral mag index)')") > | inp%AL 2135a2413,2477 > !----------------------------------------------------------------------- > subroutine validate_indices_time(rindex,ntimes,name) > ! > ! Validate user provided time series for solar index "name". Rindex will > ! be either ctpoten_time or power_time, from namelist input. > ! > ! Args: > real,intent(in) :: rindex(4,mxind_time) > integer,intent(out) :: ntimes > character(len=*),intent(in) :: name > ! > ! Local: > integer :: i,ii,n > integer(kind=8) :: nsec0,nsec1 > ! > ! External: > integer(kind=8),external :: mtime_to_nsec > ! > ! Validate times and values in user provided time series: > ntimes = 0 > if (any(rindex /= spval)) then > ! > ! Time series must be provided in groups of 4 (day,hr,min,value): > n = 0 > do i=1,mxind_time > do ii=1,4 > if (rindex(ii,i) /= spval) n = n+1 > enddo > enddo > if (mod(n,4) /= 0) then > write(6,"('>>> INPUT: must provide ',a,' in groups', > | ' of 4: n=',i5,' mod(n,4)=',i5)") name,n,mod(n,4) > stop 'validate_indices_time' > endif > do i=1,mxind_time > if (any(rindex(:,i) /= spval)) then > call validate_mtime(int(rindex(1:3,i)),mxday,name) > if (rindex(4,i) <= 0.) then > write(6,"(/,'>>> INPUT: ',a,' values must be ', > | 'positive: i=',i3,' rindex(4,i)=',e12.4)") > | name,i,rindex(4,i) > stop 'validate_indices_time' > endif > ntimes = ntimes+1 > endif > enddo ! i=1,mxind_time > ! > ! Times must be increasing (assume non-spvals are from 1 to ntimes, > ! i.e., not interspersed with spvals): > if (ntimes > 1) then > nsec0 = mtime_to_nsec(int(rindex(1:3,1))) > do i=2,ntimes > nsec1 = mtime_to_nsec(int(rindex(1:3,i))) > if (nsec0 >= nsec1) then > write(6,"(/,'>>> INPUT: ',a,' times must increase', > | '. Check time at i=',i3,' day,hr,min=',3i4,' and the', > | ' time previous to i.')") name,i,int(rindex(1:3,i)) > stop 'validate_indices_time' > endif > nsec0 = nsec1 > enddo ! i=1,ntimes > endif ! ntimes > 1 > endif ! ctpoten_time > end subroutine validate_indices_time > !----------------------------------------------------------------------- ======================================================================== Diff of /home/tgcm/tiegcm1/lamdas.F and /home/foster/tiegcm/amie/modsrc_oct03/lamdas.F /home/tgcm/tiegcm1/lamdas.F and /home/foster/tiegcm/amie/modsrc_oct03/lamdas.F are identical ======================================================================== >>> WARNING: Cannot find source file /home/tgcm/tiegcm1/low_proton.F ======================================================================== Diff of /home/tgcm/tiegcm1/magfield.F and /home/foster/tiegcm/amie/modsrc_oct03/magfield.F 15a16,18 > ! > ! Use-association of nc_open, nc_close, and handle_ncerr is not > ! done in order to avoid circular dependencies. 20d22 < use nchist_module,only: nc_open,nc_close,handle_ncerr 199,200c201,202 < call nc_open(ncid,dskfile,'OLD','READ') < if (ncid <= 0) then --- > istat = nf_open(dskfile,NF_NOWRITE,ncid) > if (istat /= NF_NOERR) then 362c364 < call nc_close(ncid) --- > istat = nf_close(ncid) 696a699,714 > subroutine handle_ncerr(istat,msg) > ! > ! Handle a netcdf lib error: > ! > integer,intent(in) :: istat > character(len=*),intent(in) :: msg > ! > write(6,"(/72('-'))") > write(6,"('>>> Error from netcdf library:')") > write(6,"(a)") trim(msg) > write(6,"('istat=',i5)") istat > write(6,"(a)") nf_strerror(istat) > write(6,"(72('-')/)") > return > end subroutine handle_ncerr > !----------------------------------------------------------------------- ======================================================================== Diff of /home/tgcm/tiegcm1/mpi.F and /home/foster/tiegcm/amie/modsrc_oct03/mpi.F 76a77 > | rtc0_periodic_f2d , rtcmp_periodic_f2d =0., 2244a2246,2354 > subroutine mp_periodic_f2d(f,lon0,lon1,lat0,lat1) > ! > ! Define periodic points for field f(lon0:lon1,lat0:lat1): > ! lons 1,2 <- nlonp4-3,nlonp4-2 and nlonp4-1,nlonp4 <- 3,4 > ! > ! Args: > integer,intent(in) :: lon0,lon1,lat0,lat1 > real,intent(inout) :: f(lon0:lon1,lat0:lat1) > ! > ! Local: > integer :: j,idest,isrc,isend,irecv,ier,len,nlons,nlats > real,allocatable :: sndbuf(:,:),rcvbuf(:,:) ! (2,nlats) > ! > ! The restriction that every task must carry at least 4 longitudes > ! means that tasks with mytidi==0 will have lons 1->4 and tasks > ! with mytidi==ntaski-1 will have lons nlonp4-3->nlonp4. > ! > if (mytidi /= 0 .and. mytidi /= ntaski-1) return > #ifdef VT > ! call vtsymdef(109, 'mp_periodic_f2d','Communication',ier) > call vtbegin(109,ier) > #endif > if (do_rtc_mpi) call timer(rtc0_periodic_f2d,tsec,'begin') > ! > ! Allocate send and receive buffers: > nlons = lon1-lon0+1 > nlats = lat1-lat0+1 > allocate(sndbuf(2,nlats),stat=ier) > if (ier /= 0) > | write(6,"('>>> mp_periodic_f2d: error allocating sndbuf:', > | ' nlats=',i3,' ier=',i4)") nlats,ier > allocate(rcvbuf(2,nlats),stat=ier) > if (ier /= 0) > | write(6,"('>>> mp_periodic_f2d: error allocating rcvbuf:', > | ' nlats=',i3,' ier=',i4)") nlats,ier > len = 2*nlats > ! > ! Send lons 3,4 to lons nlonp4-1,nlonp4 at task with mytidi==ntaski-1, > ! and receive lons 1,2 from the same task. > if (mytidi==0) then > idest = itask_table(ntaski-1,mytidj) > > do j=lat0,lat1 > sndbuf(1:2,j-lat0+1) = f(3:4,j) > enddo ! j=lat0,lat1 > > call mpi_isend(sndbuf,len,MPI_REAL8,idest,1,MPI_COMM_WORLD, > | isend,ier) > if (ier /= 0) > | call handle_mpi_err(ier,'mp_periodic_f2d send to idest') > isrc = idest > call mpi_irecv(rcvbuf,len,MPI_REAL8,isrc,1,MPI_COMM_WORLD, > | irecv,ier) > if (ier /= 0) > | call handle_mpi_err(ier,'mp_periodic_f2d recv fm isrc') > ! > ! Send lons nlonp4-3,nlonp4-2 to lons 1,2 at task with mytidi==0, > ! and receive lons nlonp4-1,nlonp4 from same task: > elseif (mytidi==ntaski-1) then > idest = itask_table(0,mytidj) > > do j=lat0,lat1 > sndbuf(1:2,j-lat0+1) = f(nlonp4-3:nlonp4-2,j) > enddo ! j=lat0,lat1 > > call mpi_isend(sndbuf,len,MPI_REAL8,idest,1,MPI_COMM_WORLD, > | isend,ier) > if (ier /= 0) > | call handle_mpi_err(ier,'mp_periodic_f2d send to idest') > isrc = idest > call mpi_irecv(rcvbuf,len,MPI_REAL8,isrc,1,MPI_COMM_WORLD, > | irecv,ier) > if (ier /= 0) > | call handle_mpi_err(ier,'mp_periodic_f2d recv fm isrc') > endif > ! > ! Wait for completions: > call mpi_wait(isend,irstat,ier) > if (ier /= 0) > | call handle_mpi_err(ier,'mp_periodic_f2d wait for send') > call mpi_wait(irecv,irstat,ier) > if (ier /= 0) > | call handle_mpi_err(ier,'mp_periodic_f2d wait for recv') > ! > ! Copy from receive buffer: > if (mytidi==0) then > do j=lat0,lat1 > f(1:2,j) = rcvbuf(1:2,j-lat0+1) > enddo ! j=lat0,lat1 > elseif (mytidi==ntaski-1) then > do j=lat0,lat1 > f(nlonp4-1:nlonp4,j) = rcvbuf(1:2,j-lat0+1) > enddo ! j=lat0,lat1 > endif > ! > ! Free local buffer space: > deallocate(sndbuf) > deallocate(rcvbuf) > ! > if (do_rtc_mpi) then > call timer(rtc0_periodic_f2d,tsec,'end') > rtcmp_periodic_f2d = rtcmp_periodic_f2d+tsec > endif > #ifdef VT > ! call vtsymdef(109, 'mp_periodic_f2d','Communication',ier) > call vtend(109,ier) > #endif > end subroutine mp_periodic_f2d > !----------------------------------------------------------------------- ======================================================================== Diff of /home/tgcm/tiegcm1/mud.F and /home/foster/tiegcm/amie/modsrc_oct03/mud.F 18c18 < subroutine mud(pe,jntl,isolve) --- > subroutine mud(pe,jntl,isolve,ier) 24c24 < integer jntl --- > integer jntl,ier ! output: not converged ier < 0 195a196,198 > ier = ierror ! < 0 not converged > if(ier < 0 ) goto 108 > 210,218c213,215 < C ITRANS = 0 < C CALL EZCNTR(PE(1,JMX0),IMX0,JMX0) < C ITRANS = 1 < C CALL SET(.05,.95,.05,.95,-1.,1.,-1.,1.,1) < C CALL CONREC(PE(1,JMX0),IMX0,IMX0,JMX0,0.,0.,0.,1,0,-1430B) < C CALL FRAME < C ITRANS = 0 < C CALL EZCNTR(PE(1,JMX0),IMX0,JMX0) < C ITRANS = 1 --- > > 108 continue > 427c424,427 < if (relmax.gt.tolmax) ierror = -1 ! flag convergenc failure --- > if (relmax.gt.tolmax) then > write(6,*) "no convergence with mud" > ierror = -1 ! flag convergenc failure > end if ======================================================================== Diff of /home/tgcm/tiegcm1/mudmod.F and /home/foster/tiegcm/amie/modsrc_oct03/mudmod.F 25c25 < subroutine mudmod(pe,jntl) --- > subroutine mudmod(pe,jntl,ier) 28c28 < integer jntl --- > integer jntl,ier ! output: not converged ier < 0 191a192,193 > ier = ierror ! ier < 0 not converged > if(ier < 0 ) goto 108 207a210 > 108 continue ======================================================================== Diff of /home/tgcm/tiegcm1/nchist.F and /home/foster/tiegcm/amie/modsrc_oct03/nchist.F 656a657 > use input_module,only: potential_model 1305a1307,1315 > ! > ! Potential model used (weimer, heelis, or none): > istat = nf_put_att_text(ncid,NF_GLOBAL,"potential_model", > | len_trim(potential_model),trim(potential_model)) > if (istat /= NF_NOERR) then > write(char120,"('Error return from nf_put_att_text ', > | 'for potential_model global attribute')") > call handle_ncerr(istat,char120) > endif ======================================================================== Diff of /home/tgcm/tiegcm1/oplus.F and /home/foster/tiegcm/amie/modsrc_oct03/oplus.F 932a933,936 > ! | phid = 5.0e8, > ! | phin = -2.0e8, > ! | phid = -3.0e8, > ! | phin = -3.0e8, 943a948 > ! if (abs(rlatm(i,lat))-pi/12.>=0.) then 946a952 > ! a(i)=.5*(1.+sin(pi*(abs(rlatm(i,lat))-pi/24.)/(pi/12.))) ======================================================================== Diff of /home/tgcm/tiegcm1/output.F and /home/foster/tiegcm/amie/modsrc_oct03/output.F 139a140,145 > ! > ! If writing to mss after model execution (dispose==2), update > ! csh dispose script after every disk write: > elseif (dispose==2) then > call savefile(ncid,diskfile,output(ioutfile),tempdir, > | dispose,.true.,iprint) 140a147 > 231a239,244 > ! > ! If writing to mss after model execution (dispose==2), update > ! csh dispose script after every disk write: > elseif (dispose==2) then > call savefile(ncidsech,diskfile,secout(isecout),tempdir, > | dispose,.true.,iprint) 237a251 > use dispose_module,only: add_dispose 272a287 > ! idispose==1 -> dispose to mss from model now. 277a293,300 > ! > ! idispose==2 -> add mswrite line to dispose script, which is executed > ! after model execution: > elseif (idispose==2) then > call add_dispose(outpath,diskfile,tempdir,' ','NCARTGCM', > | fileinfo,msreten) > ! > ! idispose==0 -> do not dispose to mss, but do link to tempdir. ======================================================================== Diff of /home/tgcm/tiegcm1/params.F and /home/foster/tiegcm/amie/modsrc_oct03/params.F 61c61,62 < | mxfsech = 50 ! max number of fields on secondary histories --- > | mxfsech = 50, ! max number of fields on secondary histories > | mxind_time = 100 ! max number of time-dependent solar index points ======================================================================== Diff of /home/tgcm/tiegcm1/putms.F and /home/foster/tiegcm/amie/modsrc_oct03/putms.F 115a116,125 > ! 12/9/02 bf: Currently, the only OSF system being used is lemieux.psc.edu. > ! Lemieux's frontend archive storage is a machine named golem.psc.edu. The > ! Below links and rcp attempt do not work from the lemieux compute nodes, > ! and may be using unnecessary wallclock time, so am not using them. Storage > ! of history output to golem will be via rcp in the job script after model > ! execution (e.g, see job script tiegcm1lx.job) > ! > write(6,"('PUTMS returning because this is an OSF job.')") > return > ! 233a244,251 > ! > ! MSS not set -- try rcp. For now, assume target is golem.psc.edu > ! (for job running on lemieux.psc.edu) > ! > #else > write(6,"('putms calling rcpfile: dskfile=',a)") > | trim(dskfile) > call rcpfile(dskfile,'golem.psc.edu:'//trim(dskfile)) ======================================================================== Diff of /home/tgcm/tiegcm1/qjoule.F and /home/foster/tiegcm/amie/modsrc_oct03/qjoule.F /home/tgcm/tiegcm1/qjoule.F and /home/foster/tiegcm/amie/modsrc_oct03/qjoule.F are identical ======================================================================== Diff of /home/tgcm/tiegcm1/qrj.F and /home/foster/tiegcm/amie/modsrc_oct03/qrj.F 252,254d251 < < if (debug) write(6,"('qrj: l=',i3,' lat=',i3)") l,lat < 397c394 < do l=1,l1-1 --- > do l=1,l1-1 ! 1->15 524a522,536 > ! > ! write(6,"(' n wave1 wave2 sflux ', > ! | 'euvflx fsrc feuv')") > ! do n=1,lmax ! 1->59 > ! if (n < l1) then ! 1->15 > ! write(6,"(i3,2f8.2,e12.4,12x,e12.4,e12.4)") > ! | n,wave1(n),wave2(n),sflux(n),fsrc(n),feuv(n) > ! elseif (n < l1+neuv-1) then ! 16->52 > ! write(6,"(i3,2f8.2,e12.4,e12.4,12x,e12.4)") > ! | n,wave1(n),wave2(n),sflux(n),euvflx(n-l1+1),feuv(n) > ! else ! 53->59 > ! write(6,"(i3,2f8.2,e12.4,12x,12x,e12.4)") > ! | n,wave1(n),wave2(n),sflux(n),feuv(n) > ! endif > ! enddo 526c538 < !--------------------------------------------------------------- --- > !----------------------------------------------------------------------- 874c886 < real,intent(out) :: euvflx(37) --- > real,intent(out) :: euvflx(neuv) 878c890,891 < real :: AFAC(37),F74113(37),flxfac --- > real :: AFAC(neuv),F74113(neuv),flxfac > real :: afac_tmp(neuv),f74113_tmp(neuv) 900a914,922 > ! 6/20/03 bf: Reverse order of F74113 and AFAC to correspond to descending > ! wavelength bins (as in wleuv1,2 and WAVEL,WAVES): > do i=1,neuv > afac_tmp(i) = afac(neuv-i+1) > f74113_tmp(i) = f74113(neuv-i+1) > enddo > afac = afac_tmp ! whole array > f74113 = f74113_tmp > ! 905c927 < DO 50 I=1,37 --- > DO 50 I=1,neuv ======================================================================== >>> WARNING: Cannot find source file /home/tgcm/tiegcm1/rgrd1.F ======================================================================== >>> WARNING: Cannot find source file /home/tgcm/tiegcm1/rgrd2.F ======================================================================== >>> WARNING: Cannot find source file /home/tgcm/tiegcm1/rgrd3.F ======================================================================== Diff of /home/tgcm/tiegcm1/tgcm.F and /home/foster/tiegcm/amie/modsrc_oct03/tgcm.F 6c6,7 < use input_module,only: input,dynamo,step --- > use input_module,only: input,dynamo,step,dispose > use dispose_module,only: init_dispose 40a42 > if (dispose==2.and.mytid==0) call init_dispose 43a46 > if (dispose==2) call init_dispose ======================================================================== Diff of /home/tgcm/tiegcm1/util.F and /home/foster/tiegcm/amie/modsrc_oct03/util.F 17c17 < if (lureq(lu)<=0) then ! do not use previously returned lu --- > ! if (lureq(lu)<=0) then ! do not use previously returned lu 21c21 < endif --- > ! endif 243a244,257 > real function finterp_bigints(f0,f1,isec0,isec1,isec) > ! > ! Do linear interpolation between f0 (which is at isec0) and > ! f1 (which is at isec1) to isec. Same as finterp except integer > ! parameters are 8-byte. > ! > ! Args: > real,intent(in) :: f0,f1 > integer(kind=8),intent(in) :: isec0,isec1,isec > ! > finterp_bigints = > | f0+(f1-f0)*float(isec-isec0)/float(isec1-isec0) > end function finterp_bigints > !------------------------------------------------------------------- 1094a1109,1144 > !----------------------------------------------------------------------- > subroutine rcpfile(source,target) > ! > ! Do remote copy of source to target. > ! > implicit none > ! > ! Args: > character(len=*),intent(in) :: source,target > ! > ! Local: > integer,parameter :: maxlen=1024 > character(len=maxlen) :: command > integer :: len,ier > ! > ! External: > integer,external :: isystem > ! > len = len_trim(source)+len_trim(target)+5 > if (len > maxlen) then > write(6,"(/,'>>> WARNING rcpfile: len(source)+len(target)+5=', > | i6,' is > maxlen=',i6)") len,maxlen > write(6,"('Please increase maxlen in sub rcpfile to at least ', > | i6,' -- no rcp done.')") len > return > endif > write(command,"('rcp ',a,' ',a)") trim(source),trim(target) > ier = isystem(command) > if (ier==0) then > write(6,"('rcpfile: successful rcp of ',a,' to ',a)") > | trim(source),trim(target) > else > write(6,"('>>> rcpfile: Error from rcp of ',a,' to ',a)") > | trim(source),trim(target) > endif > end subroutine rcpfile ======================================================================== >>> WARNING: Cannot find source file /home/tgcm/tiegcm1/wei01gcm.F