module time3d ! ! Driver module for TIME3D model, called from module dpie_coupling. ! ! TIME3D model: ! 3D Theoretical Ionospheric Model of Earth in IGGCAS ! Developed by Zhipeng Ren ! Email: zpren@mail.iggcas.ac.cn ! use shr_kind_mod ,only: r8 => shr_kind_r8 use cam_logfile ,only: iulog use time_manager ,only: get_nstep use spmd_utils ,only: masterproc use cam_abortutils ,only: endrun use savefield_edyn,only: savefld_edyn ! save fields on edynamo history use savefield_t3d ,only: savefld_t3d ! save global time3d fields to netcdf file(s) use edyn_init ,only: nstep_savefld_t3d ! interval to save time3d fields use interpolate_data,only: bilin ! bilinear interpolation routine implicit none save integer :: nstep character(len=1024) :: filename character(len=16) :: action private public time3d_drv,prep_time3d_input contains !----------------------------------------------------------------------- subroutine time3d_drv(waccm_op,waccm_z,i0,i1,j0,j1,plev) ! ! Driver for time3d model. ! use spmd_utils ,only: masterproc use time_manager ,only: get_curr_date,get_curr_calday,get_step_size use time3d_gcmsim ,only: pregrid,gcmpara,EUVFluxCrossSection,prealt,gcmsim use time3d_grid ,only: time3dparagrid use time3d_main ,only: time3d2gcm,time3dmodule use time3d_init ,only: time3dInitialCondition use time3d_flds ,only: iyd,iyr,imo,ida,f107,f107a,ap,uth,timestep,hindex use edyn_maggrid ,only: nmlat,nmlonp1 use time3d_geogrid,only: nlat,nlon,nlev use time3d_flds ,only: & ! time3d inputs: all dimensioned (nlat,nlon,nlev) Z,Tn,cO,cO2,cN2,cH,cHe,cN,Qop,QHep,OpL,Qep,QHp,HepL,HpL,cO2p, & cNOp,cN2p,Te,Ti,Ws,We,Wu use time3d_flds ,only: & ! time3d gcmsim fields Ws_msis,We_msis,Wu_msis,Tn_sim,cO_sim,cO2_sim,cN2_sim,cN_sim, & cO2p_sim,cNOp_sim,cN2p_sim,Qop_sim,Qep_sim,OpL_sim use edynamo ,only: ed1_glb,ed2_glb ! efield gathered to global on root task ! (nmlat,nmlonp1) for input to time3d use edyn_init,only: & use_time3d_integ, & ! if true, use time3d conductance in edynamo use_time3d_gcmsim, & ! if true, call time3d_gcmsim, otherwise use waccm ! fields for time3d input ! ! if use_edyn_geogrid=true, then time3d will use edyn geogrid (does not work) ! if use_edyn_maggrid=true, then time3d will use edyn magrid (does work, should be true) ! use_edyn_geogrid=>time3d_use_edyn_geogrid, & use_edyn_maggrid=>time3d_use_edyn_maggrid, & use_time3d_output ! if true, pass time3d output to waccm (Op,Te,Ti) ! ! Time3d output fields on time3d_geogrid (nlat,nlon,nlev): ! use time3d_flds ,only: Z,Te,Ti,Vi,Oppc,Tn, & Ped,Hall,Ped_t3d,Hall_t3d,Ped_waccm,Hall_waccm, & Op,Op_t3d,Op_waccm,Qep,Qop use pmgrid, only: plevp ! ! Args: ! waccm_op is overwritten only if use_time3d_output=true integer,intent(in) :: i0,i1,j0,j1,plev real(r8),dimension(i0:i1,j0:j1,plev),intent(inout) :: waccm_op ! waccm O+ real(r8),dimension(i0:i1,j0:j1,plev),intent(in) :: waccm_z ! waccm geometric ht (m) ! ! Local: ! ! Time3d global 2d outputs on mag grid (defined by sub TIME3D2Dynamo (time3d_main.f)): ! (data is on regular gmlat grid from sub pregrid) ! real(r8),dimension(nmlat,nmlonp1),save :: zigm11,zigm22,zigm2,zigmc,rim1,rim2 integer :: k,iyear,imon,iday,tod,calday real(r8) :: hb ! output by time3dparagrid character(len=80) :: string80 character(len=16) :: string16 ! ! Exec: ! nstep = get_nstep() ! ! Everybody call pregrid to set time3d mag grid (regular in latitude) ! call pregrid(use_edyn_geogrid,use_edyn_maggrid) if (.not.masterproc) goto 100 ! only masterproc calls time3d call get_curr_date(iyear,imon,iday,tod) ! tod is integer time-of-day in seconds calday = get_curr_calday() ! day of year call get_curr_date(iyear,imon,iday,tod) ! tod is integer time-of-day in seconds calday = get_curr_calday() ! day of year ! ! Set time3d GCM simulation parameters in module time3d_flds: ! f107 = 210._r8 ! F10.7 daily f107a = 210._r8 ! F10.7 average ap = 4._r8 ! AP index uth = real(tod)/3600._r8 ! UT (hours) iyd = calday ! day of year iyr = iyear ! year imo = imon ! month ida = iday ! day of month ! timestep = 300._r8 ! time step (seconds) timestep = get_step_size() ! time step (seconds) hindex(1) = 0.0_r8 hindex(2) = -0.5_r8 hindex(3) = 300._r8 hindex(4) = 10._r8 hindex(5) = 10._r8 hindex(6) = 20._r8 write(iulog,"('Enter time3d_drv: nstep=',i6,' yyyy mm ddd=',3i5,' calday=',i4,' uth=',f9.3)") & nstep,iyear,imon,iday,calday,uth ! time3d is commented for debug purposes: ! Note pregrid was called above. if (nstep==1) then call gcmpara ! time3d_gcmsim.f call time3dparagrid(iyear,calday,hb) ! time3d_grid.f call EUVFluxCrossSection(f107,f107a,calday) ! time3d_gcmsim.f call prealt(Z) ! (km) ! time3d_gcmsim.f ! write(iulog,"('time3d_drv after prealt: Z min,max=',2(1pe12.4))") & ! minval(Z),maxval(Z) call time3dInitialCondition(calday,f107,f107a,ap,uth) ! time3d_init.f ! ! Get Op,Te,Ti,Vi for GCM code ! This will also be called from sub time3dmodule at every timestep. ! ! These inits taken from Hanli's time3d/edynamo standalone: Tn(1:nlat,1:nlon,1:Nlev) = 200._r8 Oppc(1:nlat,1:nlon,1:Nlev) = 1.0e-15_r8 ! ! Call time3d2gcm in first timestep (subsequently, it is called from time3dmodule) ! Z,Oppc,Tn are input, all others are output (Z was set by prealt call above): ! call time3d2gcm (Z,Op_t3d,Te,Ti,Vi,Ped_t3d,Hall_t3d,Oppc,Tn) ! time3d_main.f Ped = Ped_t3d Hall = Hall_t3d Op = Op_t3d endif ! nstep==1 ! ! Now in waccm timestep loop. ! ! Call the time3d GCM simulator (empirical models simulate host model): ! (uth and iyd module data were set above, and are used by gcmsim) ! If use_time3d_gcmsim is true, output of this call will be passed ! on to time3dmodule, otherwise waccm inputs will be used in time3dmodule. ! ! SUBROUTINE GCMsim(Z,Tn,cO,cO2,cN2,cH,cHe,cN,Ws,We,Wu,Qop,QHep, ! | OpL,Qep,QHp,HepL,HpL,cO2p,cNOp,cN2p,Op,Te,Ti) ! call gcmsim(Z,Tn_sim,cO_sim,cO2_sim,cN2_sim,cH,cHe,cN_sim,& Ws_msis,We_msis,Wu_msis,Qop_sim,QHep,OpL_sim,Qep_sim,QHp,& HepL,HpL,cO2p_sim,cNOp_sim,cN2p_sim,Op,Te,Ti) string16 = ' ' if (nstep==1.or.use_time3d_gcmsim) then string16 = '(from gcmsim)' Ws = Ws_msis We = We_msis Wu = Wu_msis Tn = Tn_sim cO2 = cO2_sim cO = cO_sim cN2 = cN2_sim cN = cN_sim cO2p= cO2p_sim cNOp= cNOp_sim cN2p= cN2p_sim Qop = Qop_sim Qep = Qep_sim OpL = OpL_sim else string16 = '(from waccm)' ! ! Override Qop,Qep, i.e., use them from gcmsim even when ! use_time3d_gcmsim=false, since using either one of these ! from waccm will cause ion drifts to become -Infinity. ! ! 3/3/15: Setting Qop,Qep to 1.e-10 where they are zero solves ! this problem (see sub prep_time3d_input). ! ! Qop = Qop_sim ! commenting this, i.e., using waccm Qop, will cause drifts to become -Infinity ! Qep = Qep_sim ! commenting this, i.e., using waccm Qep, will cause drifts to become -Infinity endif ! ! If use_time3d_gcmsim=T, then Op,Te,Ti were returned by time3d2gcm ! (time3d2gcm is called by time3dmodule, and above, when nstep=1) ! (these fields are intent(in) in gcmsim, called above) ! ! If use_time3d_gcmsim=F, then Ped_waccm, and Hall_waccm ! were interpolated from waccm fields to global time3d grid by ! sub prep_time3d_input, and are assigned to time3d input arrays here: ! ! if (.not.use_time3d_gcmsim) then ! use fields from waccm for time3d input ! Ped = Ped_waccm ! Hall = Hall_waccm ! endif ! ! TIME3D module ! Input: TimeStep: time step; Z: altitude in meter; Tn: neutral temperature in K; ! cO,cO2,cN2,cH,cHe,cN: O, O2, N2, H, He, and N densities in m-3 ! (in now cersion, [He] is set to be 0 or 1.0D-25) ! Ws,We,Wu: southward, eastward, and upward wind velocities in m/s; ! Qop,QHep,QHp,Qep: production rates of O+,He+,H+, and photoelectron in M-3*s-1 ! (in now cersion, QHep,QHp are set to be 0 or 1.0D-25) ! OpL,HepL,HpL: loss rates of O+,He+,and H+ in s-1 ! (in now cersion, HepL,HpL are set to be 0 or 1.0D-25) ! cO2p,cNOp,cN2p: O2+, NO+, and N2+ density in m-3; (above are output from GCM) ! ed1,ed2: electric field index; (output from dynamo) ! Output: Op: O+ density in m-3; Te,Ti: electron and ion temperature in K; Vi: anti-field velocity of O+; ! (input into GCM) ! Ped,Hall: Pedersen and Hall conductivities in s/m; (input into dynamo and GCM) ! zigm11,zigm22,zigm2,zigmc,rim1,rim2; dynamo index in S or S*m/s*Be3 (input into dynamo); ! the output of [subroutine fieldline_integrals] in dynamo ! ! Ped,Hall not used by edynamo (edynamo uses sigma_ped,sigma_hall from WACCM) ! Conductances zigmxx,rim1,2 are used by edynamo if use_time3d_integ==.true. ! ! write(iulog,"('time3d_drv: time3dmodule inputs ',a,' nstep=',i5)") string16,nstep ! write(iulog,"(' Z min,max=',2e12.4)") minval(Z),maxval(Z) ! write(iulog,"(' Tn min,max=',2e12.4)") minval(Tn),maxval(Tn) ! write(iulog,"(' cO min,max=',2e12.4)") minval(cO),maxval(cO) ! write(iulog,"(' cO2 min,max=',2e12.4)") minval(cO2),maxval(cO2) ! write(iulog,"(' cN2 min,max=',2e12.4)") minval(cN2),maxval(cN2) ! write(iulog,"(' cH min,max=',2e12.4)") minval(cH),maxval(cH) ! write(iulog,"(' cHe min,max=',2e12.4)") minval(cHe),maxval(cHe) ! write(iulog,"(' cN min,max=',2e12.4)") minval(cN),maxval(cN) ! write(iulog,"(' Ws min,max=',2e12.4)") minval(Ws),maxval(Ws) ! write(iulog,"(' We min,max=',2e12.4)") minval(We),maxval(We) ! write(iulog,"(' Wu min,max=',2e12.4)") minval(Wu),maxval(Wu) ! msis Wu is all zero ! write(iulog,"(' Qop min,max=',2e12.4)") minval(Qop),maxval(Qop) ! write(iulog,"(' QHep min,max=',2e12.4)") minval(QHep),maxval(QHep) ! write(iulog,"(' OpL min,max=',2e12.4)") minval(OpL),maxval(OpL) ! write(iulog,"(' QeP min,max=',2e12.4)") minval(QeP),maxval(QeP) ! write(iulog,"(' QHP min,max=',2e12.4)") minval(QHP),maxval(QHP) ! write(iulog,"(' HepL min,max=',2e12.4)") minval(HepL),maxval(HepL) ! write(iulog,"(' HpL min,max=',2e12.4)") minval(HpL),maxval(HpL) ! write(iulog,"(' cO2p min,max=',2e12.4)") minval(cO2p),maxval(cO2p) ! write(iulog,"(' cNOp min,max=',2e12.4)") minval(cNOp),maxval(cNOp) ! write(iulog,"(' cN2p min,max=',2e12.4)") minval(cN2p),maxval(cN2p) ! write(iulog,"(' Op min,max=',2e12.4)") minval(Op),maxval(Op) ! write(iulog,"(' Ped min,max=',2e12.4)") minval(Ped),maxval(Ped) ! write(iulog,"(' Hall min,max=',2e12.4)") minval(Hall),maxval(Hall) ! write(iulog,"(' ed1_glb min,max=',2e12.4)") minval(ed1_glb),maxval(ed1_glb) ! write(iulog,"(' ed2_glb min,max=',2e12.4)") minval(ed2_glb),maxval(ed2_glb) ! ! Save time3d inputs to netcdf file: ! if (time2write(nstep,nstep_savefld_t3d)) then filename = 'time3d_waccm_input.nc' if (use_time3d_gcmsim) filename = 'time3d_gcmsim_input.nc' action = 'append' if (nstep==1) action = 'create' call savefld_t3d(filename,'Op','O+ '//string16,'m-3',& Op,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) action = 'append' ! call savefld_t3d(filename,'Tn','Tn '//string16,'deg K',& ! Tn,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) ! call savefld_t3d(filename,'Ws','Ws '//string16,'m/s',& ! Ws,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) ! call savefld_t3d(filename,'We','We '//string16,'m/s',& ! We,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) ! call savefld_t3d(filename,'Wu','Wu '//string16,'m/s',& ! Wu,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) ! call savefld_t3d(filename,'cO','O '//string16,'m-3',& ! cO,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) ! call savefld_t3d(filename,'cO2','O2 '//string16,'m-3',& ! cO2,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) ! call savefld_t3d(filename,'cN2','N2 '//string16,'m-3',& ! cN2,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) ! call savefld_t3d(filename,'cN','N '//string16,'m-3',& ! cN,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) ! call savefld_t3d(filename,'cO2p','O2+ '//string16,'m-3',& ! cO2p,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) ! call savefld_t3d(filename,'cNOp','NO+ '//string16,'m-3',& ! cNOp,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) ! call savefld_t3d(filename,'cN2p','N2+ '//string16,'m-3',& ! cN2p,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) ! call savefld_t3d(filename,'Op','O+ '//string16,'m-3',& ! Op,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) call savefld_t3d(filename,'Qop','O+ Production Rate'//string16,'m-3/s',& Qop,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) call savefld_t3d(filename,'Qep','PE Production Rate'//string16,'m-3/s',& Qep,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) call savefld_t3d(filename,'OpL','O+ Loss Rate'//string16,'m-3/s',& OpL,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) ! call savefld_t3d(filename,'Ped','Pedersen '//string16,' ',& ! Ped,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) ! call savefld_t3d(filename,'Hall','Hall '//string16,' ',& ! Hall,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) ! call savefld_t3d(filename,'Z','Z geopotential','m',& ! Z,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) endif ! ! Save efield (input to time3d from edynamo (will be zero in first step)) ! ! if (time2write(nstep,nstep_savefld_t3d)) then ! filename = 'time3d_edyn_maginput.nc' ! action = 'append' ! if (nstep==1) action = 'create' ! call savefld_t3d(filename,'ed1','ed1 time3d input from edynamo','()',& ! ed1_glb,(/nmlat,nmlonp1,1/),(/'nlat','nlon','nlev'/),action,nstep) ! action = 'append' ! call savefld_t3d(filename,'ed2','ed2 time3d input from edynamo','()',& ! ed2_glb,(/nmlat,nmlonp1,1/),(/'nlat','nlon','nlev'/),action,nstep) ! endif ! ! Call main time3d driver: ! call time3dmodule(timestep,Z,Tn,cO,cO2,cN2,cH,cHe,cN,Ws,We,Wu, & Qop,QHep,OpL,Qep,QHp,HepL,HpL,cO2p,cNOp,cN2p,ed1_glb,ed2_glb,& Op,Te,Ti,Vi,Ped,Hall,zigm11,zigm22,zigm2,zigmc,rim1,rim2) ! ! Save time3d output on global 3d geo grid: ! ! write(iulog,"('time3d_drv: time3dmodule output: nstep=',i5)") nstep ! write(iulog,"(' Op min,max=',2e12.4)") minval(Op),maxval(Op) ! write(iulog,"(' Te min,max=',2e12.4)") minval(Te),maxval(Te) ! write(iulog,"(' Ti min,max=',2e12.4)") minval(Ti),maxval(Ti) ! ! Save time3d output on global 2d mag grid: ! ! write(iulog,"(' zigm11 min,max=',2e12.4)") minval(zigm11),maxval(zigm11) ! write(iulog,"(' zigm22 min,max=',2e12.4)") minval(zigm22),maxval(zigm22) ! write(iulog,"(' zigm2 min,max=',2e12.4)") minval(zigm2) ,maxval(zigm2) ! write(iulog,"(' zigmc min,max=',2e12.4)") minval(zigmc) ,maxval(zigmc) ! write(iulog,"(' rim1 min,max=',2e12.4)") minval(rim1) ,maxval(rim1) ! write(iulog,"(' rim2 min,max=',2e12.4)") minval(rim2) ,maxval(rim2) ! ! Save time3d output to netcdf file: ! if (time2write(nstep,nstep_savefld_t3d)) then filename = 'time3d_waccm_output.nc' if (use_time3d_gcmsim) filename = 'time3d_gcmsim_output.nc' action = 'append' if (nstep==1) action = 'create' call savefld_t3d(filename,'Op','O+ (time3d output) ','m-3',& Op,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) action = 'append' call savefld_t3d(filename,'Te','Te (time3d output)','deg K',& Te,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) call savefld_t3d(filename,'Ti','Ti (time3d output)','deg K',& Ti,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) endif 100 continue ! if use_time3d_integ is set, then slave tasks will ! wait in sub prep_conductances for the root task to catch up ! ! If using conductances from time3d in edynamo, sub prep_conductances will ! transform the global arrays from (nmlat,nmlonp1) to (nmlonp1,nmlat), and ! scatter them to the mag subdomains so edynamo can call complete_integrals. ! Edynamo will then gather them back to the global domain for the potential ! solver (in this case edynamo will not call fieldline_integrals). ! ! If use_edyn_maggrid=false, then time3d used its own mag grid (regular in ! latitude, 0->360 in longitude). In this case, sub prep_conductances will ! interpolate to the irregular edynamo mag lats and phase shift to -180to180. ! if (use_time3d_integ) then call prep_conductances(zigm11,zigm22,zigm2,zigmc,rim1,rim2,nmlat,nmlonp1) endif ! ! Write global time3d conductances to nc file. ! If .not.use_edyn_maggrid, then these have been interpolated to the edyn ! maggrid, otherwise they remain on the time3d maggrid. ! ! if (masterproc) then if (time2write(nstep,nstep_savefld_t3d)) then filename = 'time3d_waccm_magoutput.nc' if (use_time3d_gcmsim) filename = 'time3d_gcmsim_magoutput.nc' action = 'append' if (nstep==1) action = 'create' call savefld_t3d(filename,'zigm11','time3d zigm11','()',& zigm11,(/nmlat,nmlonp1,1/),(/'nlat','nlon','nlev'/),action,nstep) action = 'append' call savefld_t3d(filename,'zigm22','time3d zigm22','()',& zigm22,(/nmlat,nmlonp1,1/),(/'nlat','nlon','nlev'/),action,nstep) call savefld_t3d(filename,'zigm2','time3d zigm2','()',& zigm2,(/nmlat,nmlonp1,1/),(/'nlat','nlon','nlev'/),action,nstep) call savefld_t3d(filename,'zigmc','time3d zigmc','()',& zigmc,(/nmlat,nmlonp1,1/),(/'nlat','nlon','nlev'/),action,nstep) call savefld_t3d(filename,'rim1','time3d rim1','()',& rim1,(/nmlat,nmlonp1,1/),(/'nlat','nlon','nlev'/),action,nstep) call savefld_t3d(filename,'rim2','time3d rim2','()',& rim2,(/nmlat,nmlonp1,1/),(/'nlat','nlon','nlev'/),action,nstep) endif ! time to write to nc file endif ! masterproc ! ! If use_time3d_output is true, we pass time3d output O+,Te,Ti to waccm ! This could be with use_time3d_gcmsim T or F, but is more consistent ! with use_time3d_gcmsim=false, i.e., time3d using waccm input, esp Op. ! Also note that Op is intent(inout) in time3d, and Te,Ti are intent(out). ! ! Op,Te,Ti are in time3d_flds (nlat,nlon,nlev) (available only on masterproc) ! Sub prep_time3d_out will interpolate to waccm global geo grid, and scatter ! to waccm block decomposition, defining waccm arrays on subdomains. ! if (use_time3d_output) then call prep_time3d_output(Op,Te,Ti,nlat,nlon,nlev, & waccm_op,waccm_z,i0,i1,j0,j1) endif end subroutine time3d_drv !----------------------------------------------------------------------- subroutine prep_time3d_output(op_t3d,te_t3d,ti_t3d,nlat,nlon,nlev,& waccm_op,waccm_z,i0,i1,j0,j1) use pmgrid ,only: & plat,plon,plev,plevp ! waccm global grid dimensions use time3d_geogrid,only: & t3d_londeg=>londeg, & ! time3d global longitudes (deg) (nlon+1 includes periodic point) t3d_latdeg_n2s=>latdeg, & ! time3d global latitudes (deg) (north to south) alt ! time3d heights (nlev) use commap ,only: & ! waccm grid info waccm_latdeg=>latdeg, & ! global latitude coordinate (plat) waccm_londeg=>londeg ! global longitude coordinate (plon,plat) use edyn_mpi ,only: mp_scatter_2d ! scatter solution to slave tasks use cam_history ,only: outfld ! CAM routine to output fields to history files use htinterp ,only: interp_column ! column interpolator ! ! Args: ! ! Input fields on global time3d geographic grid: integer,intent(in) :: nlat,nlon,nlev ! global time3d geo grid real(r8),dimension(nlat,nlon,nlev),intent(in) :: & op_t3d,te_t3d,ti_t3d ! ! Input waccm 3d: real(r8),dimension(i0:i1,j0:j1,plev),intent(in) :: waccm_z ! (m) ! ! Output fields on distributed waccm geographic grid: ! (intent(inout) because it is needed at levels below the time3d bottom boundary) ! integer,intent(in) :: i0,i1,j0,j1 ! waccm subdomain boundaries real(r8),dimension(i0:i1,j0:j1,plev),intent(inout) :: & waccm_op ! ! Local: integer :: i,j,k,n,logint,kout,k90_waccm real(r8) :: zmax_waccm integer :: koutplon(plon) integer :: plon_out(plon) integer,parameter :: nf=3 ! op,te,ti character(len=8) :: fnames(nf) = (/'op','te','ti'/) real(r8) :: r,aveop(plev),avez(plev) ! ! Local on time3d grid: real(r8) :: t3d_latdeg(nlat) real(r8) :: t3d_ikj(nlon,nlev,nlat) ! ! Local on waccm global grid: real(r8) :: waccm_glb_t3dk (plon,plat,nlev) real(r8) :: waccm_glb_ijk (plon,plat,plev) real(r8) :: waccm_glb_ijk_z (plon,plat,plev) real(r8) :: waccm_glb_ijk_op(plon,plat,plev) ! ! Local on waccm subdomains: real(r8) :: waccm_blk_ij (i0:i1,j0:j1) real(r8) :: waccm_blk_ik (i0:i1,plev) real(r8) :: waccm_blk_ijk(i0:i1,j0:j1,plev) real(r8) :: waccm_blk_ijk_op(i0:i1,j0:j1,plev) real(r8) :: diag_ik(i0:i1,plev) ! ! Save waccm input field for testing,debug. This will be restored ! to the original input below for rest of the run, until this routine ! is ready for testing. ! ! write(iulog,"('prep_time3d_out: waccm_op (input on waccm subdomains): min,max =',2e12.4)") & ! minval(waccm_op),maxval(waccm_op) waccm_blk_ijk_op = waccm_op ! save for testing,debug ! ! Gather waccm geometric height to global: ! (waccm_glb_ijk_z(plon,plat,plev) is output). ! call gather_field(waccm_z,i0,i1,j0,j1,waccm_glb_ijk_z,plon,plat,plev,.true.,0) ! ! Gather waccm inout fields as they will be needed globally below ! the time3d bottom boundary. ! call gather_field(waccm_blk_ijk_op,i0,i1,j0,j1,waccm_glb_ijk_op,plon,plat,plev,.true.,0) if (.not.masterproc) goto 100 ! write(iulog,"('Enter prep_time3d_output: nstep=',i4)") nstep ! write(iulog,"(' op_t3d min,max =',2e12.4)") minval(op_t3d),maxval(op_t3d) ! write(iulog,"(' te_t3d min,max =',2e12.4)") minval(te_t3d),maxval(te_t3d) ! write(iulog,"(' ti_t3d min,max =',2e12.4)") minval(ti_t3d),maxval(ti_t3d) ! write(iulog,"(' waccm_glb_ijk_z min,max=',2e12.4)") & ! minval(waccm_glb_ijk_z),maxval(waccm_glb_ijk_z) ! ! Set up files for post-plotting: ! ! time3d output (to be interpolated to waccm grid for input to waccm): if (time2write(nstep,nstep_savefld_t3d)) then filename = 'Begin_prep_time3d_out.nc' action = 'append' if (nstep==1) action = 'create' call savefld_t3d(filename,'op_t3d','O+ out of time3d (N2S,bot2top)',& '()',op_t3d,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) action = 'append' ! call savefld_t3d(filename,'op_t3d','time3d te (N2S) (begin prep_time3d_out)',& ! '()',op_te,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) ! call savefld_t3d(filename,'op_t3d','time3d ti (N2S) (begin prep_time3d_out)',& ! '()',op_ti,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) endif ! ! Initial gathered waccm O+ and Z (these will be interpolated to time3d grid): ! if (time2write(nstep,nstep_savefld_t3d)) then ! filename = 'waccm_prep_time3d_out.nc' ! action = 'append' ! if (nstep==1) action = 'create' ! call savefld_t3d(filename,'waccm_z','waccm Z (top2bot) (begin prep_time3d_out)',& ! '()',waccm_glb_ijk_z,(/plon,plat,plev/),(/'nlon','nlat','nlev'/),action,nstep) ! action = 'append' ! call savefld_t3d(filename,'waccm_op','waccm O+ (top2bot) (begin prep_time3d_out)',& ! '()',waccm_glb_ijk_op,(/plon,plat,plev/),(/'nlon','nlat','nlev'/),action,nstep) ! endif ! ! Reverse vertical order of waccm Z so bottom boundary is at k==1, top is at k=plev. ! Also convert Z from m to km. ! waccm_glb_ijk_z(:,:,:) = waccm_glb_ijk_z(:,:,plev:1:-1)*1.e-3 ! ! Find waccm level on which waccm Z is globally < 90 km (bottom time3d level) k90_waccm = 0 do k=1,plev if (all(waccm_glb_ijk_z(:,:,k) > 90._r8)) then k90_waccm = k exit endif enddo if (k90_waccm==0) call endrun('prep_time3d_out could not find k90_waccm') ! write(iulog,"('prep_time3d_out: k90_waccm=',i4,' (plev=',i4,')')") k90_waccm,plev ! ! Find max top boundary of waccm (km): zmax_waccm = maxval(waccm_glb_ijk_z(:,:,plev)) ! write(iulog,"('prep_time3d_out: zmax_waccm=',f12.4)") zmax_waccm ! ! Look at global average waccm O+ and Z at each waccm level (note k=1,plev is "top2bot"): ! r = 1._r8/(plon*plat) ! do k=1,plev ! aveop(k) = 0._r8 ! avez(k) = 0._r8 ! do j=1,plat ! do i=1,plon ! aveop(k) = aveop(k) + waccm_glb_ijk_op(i,j,k) * r ! avez(k) = avez(k) + waccm_glb_ijk_z (i,j,k) * r ! enddo ! enddo ! write(iulog,"('k=',i3,' aveop(k)=',es12.4,' avez(k)=',es12.4)") & ! k,aveop(k),avez(k) ! enddo ! ! Reorder time3d latitude coord from north->south (N2S) to south->north (S2N): do j=1,nlat t3d_latdeg(j) = t3d_latdeg_n2s(nlat-j+1) enddo ! write(iulog,"('prep_time3d_out: nlat=',i3,' t3d_latdeg (S2N)=',/,(8f9.2))") & ! nlat,t3d_latdeg plon_out(:) = plon ! ! Prepare time3d inputs for bilin interpolation: ! (reorder data in latitude coord from N2S to S2N) ! 100 continue do n=1,nf if (.not.masterproc) goto 101 select case(trim(fnames(n))) case('op') do k=1,nlev do j=1,nlat t3d_ikj(:,k,j) = op_t3d(nlat-j+1,:,k) enddo enddo logint = 1 ! if (time2write(nstep,nstep_savefld_t3d)) then ! filename = 't3d_ikj_prep_time3d_out.nc' ! action = 'append' ! if (nstep==1) action = 'create' ! call savefld_t3d(filename,'t3d_ikj_op','time3d O+ (S2N) (prep_time3d_out)',& ! '()',t3d_ikj,(/nlon,nlev,nlat/),(/'nlon','nlev','nlat'/),action,nstep) ! endif case('te') do k=1,nlev do j=1,nlat t3d_ikj(:,k,j) = te_t3d(nlat-j+1,:,k) enddo enddo logint = 0 case('ti') do k=1,nlev do j=1,nlat t3d_ikj(:,k,j) = ti_t3d(nlat-j+1,:,k) enddo enddo logint = 0 case default write(iulog,"('>>> prep_time3d_output: unknown fname=',a)") trim(fnames(n)) call endrun end select ! ! Horizontal interpolation from time3d geo grid to waccm geo grid, at nlev levels: ! (no height or pressure interpolation yet) ! Bilin input: t3d_ikj(nlon,nlev,nlat) ! Bilin output: waccm_glb_t3dk (plon,plat,nlev) ! (bot2top t3d O+ interpolated to waccm horizontal grid) ! waccm_glb_t3dk = 0._r8 ! (plon,nlev,plat) do k=1,nlev call bilin(t3d_ikj(:,k,:),t3d_londeg,t3d_latdeg,nlon,nlon,1,1,nlat, & waccm_glb_t3dk(:,:,k),waccm_londeg,waccm_latdeg,plon,plon_out,1,plat) enddo ! k=1,nlev ! if (time2write(nstep,nstep_savefld_t3d).and.trim(fnames(n))=='op') then ! filename = 'waccm_ijt3dk_prep_time3d_out.nc' ! action = 'append' ! if (nstep==1) action = 'create' ! call savefld_t3d(filename,'waccm_t3dk_bilin_op','waccm O+ (bot2top) (after bilin)',& ! '()',waccm_glb_t3dk,(/plon,plat,nlev/),(/'nlon','nlat','nlev'/),action,nstep) ! endif ! ! Now interpolate from time3d heights to waccm pressure levels. ! This will define waccm_glb_ijk from k90_waccm (level at ~90 km) ! upward to waccm upper boundary. ! ! Note waccm_glb_ijk_z is "bottom2top" (see above), and for now ! waccm_glb_ijk is also "bottom2top", but will be inverted to ! "top2bottom" in final prep for input to waccm. ! ! subroutine interp_column(fin,levin,nlevin,fout,levout,nlevout,& ! logint,iprint) ! Interp_column input: waccm_glb_t3dk(i,j,1:nlev) ! Interp_column output: waccm_glb_ijk(i,j,k90_waccm:plev) ! do j=1,plat do i=1,plon call interp_column(waccm_glb_t3dk(i,j,:),alt,nlev,& waccm_glb_ijk(i,j,k90_waccm:plev),waccm_glb_ijk_z(i,j,k90_waccm:plev),& plev-k90_waccm+1,logint,0) enddo ! i=1,plon enddo ! j=1,plat ! ! Assign data at levels from k90_waccm to bottom boundary from original ! waccm input waccm_glb_ijk_xx (note waccm_glb_ijk_xx is still in "top2bot" ! order. Waccm_glb_ijk is in "bot2top" order, but will be reversed in ! the final array, ready to be input to waccm). ! select case(trim(fnames(n))) case('op') do k=1,k90_waccm waccm_glb_ijk(:,:,k) = waccm_glb_ijk_op(:,:,plev-k+1) enddo case('te') case('ti') end select if (time2write(nstep,nstep_savefld_t3d)) then filename = 'End_prep_time3d_out.nc' ! was created above at nstep=1 action = 'append' if (nstep==1.and.n==1) action = 'create' select case(trim(fnames(n))) case('op') call savefld_t3d(filename,'waccm_glb_ijk_op','waccm O+ on time3d grid (S2N,bot2top)',& '()',waccm_glb_ijk,(/plon,plat,plev/),(/'nlon','nlat','nlev'/),action,nstep) action = 'append' case('te') case('ti') end select endif ! ! Reverse order of waccm_glb_ijk from "bot2top" to "top2bot" for waccm: waccm_glb_ijk(:,:,1:plev) = waccm_glb_ijk(:,:,plev:1:-1) ! if (time2write(nstep,nstep_savefld_t3d)) then ! filename = 'End_prep_time3d_out.nc' ! action = 'append' ! if (nstep==1) action = 'create' ! select case(trim(fnames(n))) ! case('op') ! call savefld_t3d(filename,'OP','O+ in to waccm (S2N,top2bot)',& ! '()',waccm_glb_ijk,(/plon,plat,plev/),(/'nlon','nlat','nlev'/),action,nstep) ! action = 'append' ! case('te') ! case('ti') ! end select ! endif ! 101 continue ! all procs ! ! Scatter to waccm subdomains: ! subroutine mp_scatter_2d(fsub,ilon0,ilon1,ilat0,ilat1,fglb,nglblon,nglblat,nf) ! real(r8) :: waccm_blk_ijk(i0:i1,j0:j1,plev) ! do k=1,plev call mp_scatter_2d(waccm_blk_ijk(:,:,k),i0,i1,j0,j1,& waccm_glb_ijk(:,:,k),plon,plat,1) ! ! Assign to output arrays (i0:i1,j0:j1,k): ! real(r8),dimension(i0:i1,j0:j1,plev),intent(inout) :: waccm_op ! select case(trim(fnames(n))) case('op') waccm_op(:,:,k) = waccm_blk_ijk(:,:,k) case('te') ! waccm_te(i,j,k) = waccm_blk_ijk(:,:,k) case('ti') ! waccm_ti(i,j,k) = waccm_blk_ijk(:,:,k) end select enddo ! k=1,plev enddo ! n=1,nf ! write(iulog,"('prep_time3d_output returning: waccm_op (interp from time3d Op) min,max=',2e12.4)") & ! minval(waccm_op),maxval(waccm_op) ! ! Save fields prepared for waccm input to waccm history file: ! do j=j0,j1 do i=i0,i1 diag_ik(i,:) = waccm_op(i,j,:) enddo ! ! See sub add_fields in edyn_init.F90 for names. These names also have ! to be in namelist nl_cam to be written to the history. ! call outfld('T3D_OP_2WACCM',diag_ik,i1-i0+1,j) enddo ! j=j0,j1 ! ! Remove this when ready to test: ! waccm_op = waccm_blk_ijk_op ! save for testing,debug end subroutine prep_time3d_output !----------------------------------------------------------------------- subroutine prep_conductances(zigm11_glbin,zigm22_glbin,zigm2_glbin,& zigmc_glbin,rim1_glbin,rim2_glbin,nmlat,nmlonp1) ! ! Prepare time3d output conductances for input to edynamo (use_time3d_integ=T) ! ! Input args are global mag grid from time3d (on masterproc only). ! On output, edynamo zigmxx,rim1,2 are scattered to mag subdomains, ! ready for input to the edynamo. ! If NOT use_edyn_maggrid, then the conductances have been calculated ! with time3d regular mag lat grid (gmlat_t3d), and gmlon_t3d = 0->360. ! In this case, interpolate to the edynamo mag grid (irregular in mlat), ! and phase shift to -180->+180. ! use edyn_mpi ,only: mp_scatter_2d ! scatter solution to slave tasks use edyn_mpi ,only: mlon0,mlon1,mlat0,mlat1 use edyn_init ,only: lonshift_global ! routine to phase shift global longitude data use edynamo ,only: & ! output on mag subdomains (mlon0:mlon1,mlat0:mlat1) zigm11, & ! sigma11*cos(theta0) zigmc, & ! sigmac zigm2, & ! sigma2 zigm22, & ! sigma22/cos(theta0) rim1,rim2 ! use edyn_init,only: & use_edyn_maggrid=>time3d_use_edyn_maggrid, use_time3d_gcmsim use edyn_maggrid,only: nmlon,gmlat,gmlon use time3d_grid,only: gmlat_t3d,gmlon_t3d ! ! Args (input global conductivities from time3d): integer,intent(in) :: nmlat,nmlonp1 real(r8),dimension(nmlat,nmlonp1),intent(in) :: & ! input zigm11_glbin , & zigm22_glbin , & zigm2_glbin , & zigmc_glbin , & rim1_glbin,rim2_glbin ! ! Local: integer :: i,j,n real(r8),dimension(nmlonp1,nmlat) :: & zigm11_glbinx , & zigm22_glbinx , & zigm2_glbinx , & zigmc_glbinx , & rim1_glbinx,rim2_glbinx real(r8),dimension(nmlonp1,nmlat) :: time3d_glb_ij, edyn_glb_ij real(r8),dimension(nmlonp1,nmlat,6) :: fglb ! global fields real(r8),dimension(mlon0:mlon1,mlat0:mlat1,6) :: fsub ! subdomain fields real(r8) :: gmlon_shifted(nmlon,nmlat) ! edynamo lon coord shifted to 0->360 integer :: nmlon_out(nmlat) ! number of mag lons per mag lat nstep = get_nstep() ! write(iulog,"('Enter prep_conductances: nstep=',i5)") nstep if (masterproc) then ! write(iulog,"(' zigm11_glbin=',2e12.4)") minval(zigm11_glbin),maxval(zigm11_glbin) ! write(iulog,"(' zigm22_glbin=',2e12.4)") minval(zigm22_glbin),maxval(zigm22_glbin) ! write(iulog,"(' zigm2_glbin =',2e12.4)") minval(zigm2_glbin) ,maxval(zigm2_glbin) ! write(iulog,"(' zigmc_glbin =',2e12.4)") minval(zigmc_glbin) ,maxval(zigmc_glbin) ! write(iulog,"(' rim1_glbin =',2e12.4)") minval(rim1_glbin) ,maxval(rim1_glbin) ! write(iulog,"(' rim2_glbin =',2e12.4)") minval(rim2_glbin) ,maxval(rim2_glbin) if (.not.use_edyn_maggrid) then ! ! Prepare lon coords for bilin: do j=1,nmlat gmlon_shifted(:,j) = gmlon(1:nmlon) call lonshift_global(gmlon_shifted(:,j),nmlon,'zeroto360',.true.) nmlon_out(j) = nmlon enddo ! ! Loop over fields: do n=1,6 select case(n) case(1) do j=1,nmlat time3d_glb_ij(:,j) = zigm11_glbin(j,:) enddo case(2) do j=1,nmlat time3d_glb_ij(:,j) = zigm22_glbin(j,:) enddo case(3) do j=1,nmlat time3d_glb_ij(:,j) = zigm2_glbin(j,:) enddo case(4) do j=1,nmlat time3d_glb_ij(:,j) = zigmc_glbin(j,:) enddo case(5) do j=1,nmlat time3d_glb_ij(:,j) = rim1_glbin(j,:) enddo case(6) do j=1,nmlat time3d_glb_ij(:,j) = rim2_glbin(j,:) enddo end select ! ! Interpolate from time3d mag grid (regular in latitude) to edynamo mag grid (irreg in lat): ! (Data and lon coords are at 0->360) ! call bilin(time3d_glb_ij(1:nmlon,:),gmlon_t3d,gmlat_t3d,nmlon,nmlon,1,1,nmlat, & edyn_glb_ij(1:nmlon,:),gmlon_shifted,gmlat,nmlon,nmlon_out,1,nmlat) ! ! Lon shift interpolated data from 0->360 to -180to180 for edynamo: ! do j=1,nmlat call lonshift_global(edyn_glb_ij(1:nmlon,j),nmlon,'-180to180',.false.) edyn_glb_ij(nmlonp1,j) = edyn_glb_ij(1,j) ! periodic point enddo select case(n) case(1) do j=1,nmlat zigm11_glbinx(:,j) = edyn_glb_ij(:,j) enddo case(2) do j=1,nmlat zigm22_glbinx(:,j) = edyn_glb_ij(:,j) enddo case(3) do j=1,nmlat zigm2_glbinx(:,j) = edyn_glb_ij(:,j) enddo case(4) do j=1,nmlat zigmc_glbinx(:,j) = edyn_glb_ij(:,j) enddo case(5) do j=1,nmlat rim1_glbinx(:,j) = edyn_glb_ij(:,j) enddo case(6) do j=1,nmlat rim2_glbinx(:,j) = edyn_glb_ij(:,j) enddo end select enddo ! n=1,6 else ! time3d used edynamo mag grid do j=1,nmlat do i=1,nmlonp1 zigm11_glbinx(i,j) = zigm11_glbin(j,i) zigm22_glbinx(i,j) = zigm22_glbin(j,i) zigmc_glbinx (i,j) = zigmc_glbin (j,i) zigm2_glbinx (i,j) = zigm2_glbin (j,i) rim1_glbinx (i,j) = rim1_glbin (j,i) rim2_glbinx (i,j) = rim2_glbin (j,i) enddo enddo endif ! use_edyn_maggrid or not ! ! Scatter time3d zigmxx and rim1,2 back to subdomains so complete_integrals ! can be called from edynamo: ! fglb(:,:,1) = zigm11_glbinx(:,:) fglb(:,:,2) = zigm22_glbinx(:,:) fglb(:,:,3) = zigm2_glbinx(:,:) fglb(:,:,4) = zigmc_glbinx(:,:) fglb(:,:,5) = rim1_glbinx(:,:) fglb(:,:,6) = rim2_glbinx(:,:) endif ! masterproc ! ! fglb is input (defined on masterproc only), fsub is output (all procs on subdomains): ! call mp_scatter_2d(fsub,mlon0,mlon1,mlat0,mlat1,fglb,nmlonp1,nmlat,6) ! ! Note zigmxx,rim1,2 are allocated with halos. Here we ignore halos. ! zigm11(mlon0:mlon1,mlat0:mlat1) = fsub(:,:,1) zigm22(mlon0:mlon1,mlat0:mlat1) = fsub(:,:,2) zigm2 (mlon0:mlon1,mlat0:mlat1) = fsub(:,:,3) zigmc (mlon0:mlon1,mlat0:mlat1) = fsub(:,:,4) rim1 (mlon0:mlon1,mlat0:mlat1) = fsub(:,:,5) rim2 (mlon0:mlon1,mlat0:mlat1) = fsub(:,:,6) ! write(iulog,"('prep_conductances: time3d integrations on mag subdomains: nstep=',i4)") & ! nstep ! write(iulog,"(' zigm11 min,max=',2e12.4)") minval(zigm11(mlon0:mlon1,mlat0:mlat1)),& ! maxval(zigm11(mlon0:mlon1,mlat0:mlat1)) ! write(iulog,"(' zigm22 min,max=',2e12.4)") minval(zigm22(mlon0:mlon1,mlat0:mlat1)),& ! maxval(zigm22(mlon0:mlon1,mlat0:mlat1)) ! write(iulog,"(' zigm2 min,max=',2e12.4)") minval(zigm2 (mlon0:mlon1,mlat0:mlat1)),& ! maxval(zigm2 (mlon0:mlon1,mlat0:mlat1)) ! write(iulog,"(' zigmc min,max=',2e12.4)") minval(zigmc (mlon0:mlon1,mlat0:mlat1)),& ! maxval(zigmc (mlon0:mlon1,mlat0:mlat1)) ! write(iulog,"(' rim1 min,max=',2e12.4)") minval(rim1 (mlon0:mlon1,mlat0:mlat1)),& ! maxval(rim1 (mlon0:mlon1,mlat0:mlat1)) ! write(iulog,"(' rim2 min,max=',2e12.4)") minval(rim2 (mlon0:mlon1,mlat0:mlat1)),& ! maxval(rim2 (mlon0:mlon1,mlat0:mlat1)) ! ! Save to edynamo.nc history file (exclude halo points): ! subroutine savefld_edyn(name,long_name,units,f,dname1,lb1,ub1,dname2,lb2,ub2,idx) ! (See also savefld_edyn calls in edynamo.F90) ! ! call savefld_edyn('T3D2EDYN_ZIGM11','ZIGM11 at end of prep_conductances',' ',& ! zigm11(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call savefld_edyn('T3D2EDYN_ZIGM22','ZIGM22 at end of prep_conductances',' ',& ! zigm22(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call savefld_edyn('T3D2EDYN_ZIGM2','ZIGM2 at end of prep_conductances',' ',& ! zigm2(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call savefld_edyn('T3D2EDYN_ZIGMC','ZIGMC at end of prep_conductances',' ',& ! zigmc(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call savefld_edyn('T3D2EDYN_RIM1','RIM1 at end of prep_conductances',' ',& ! rim1(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) ! call savefld_edyn('T3D2EDYN_RIM2','RIM2 at end of prep_conductances',' ',& ! rim2(mlon0:mlon1,mlat0:mlat1),'mlon',mlon0,mlon1,'mlat',mlat0,mlat1,0) end subroutine prep_conductances !----------------------------------------------------------------------- subroutine prep_time3d_input(z,t,u,v,w,o2,o1,n2,n1,o2p,nop,n2p,op,& sigma_ped,sigma_hall,qepp,qopp,lop,i0,i1,j0,j1,plev,nflds,fnames) ! ! Prepare waccm and edynamo fields for input to time3d. ! This is called from dpie_coupling before time3d_drv: ! ! call prep_time3d_input(zht,tn,u,v,w,o2,o1,n2,n1,o2p,nop,n2p,op, & ! sigma_ped,sigma_hall,qep,qop,lop,i0,i1,j0,j1,plev,18, & ! (/'z','t','u','v','w','o2','o','n2','n','o2p','nop','n2p','op',& ! 'Ped','Hall','Qep','Qop','OpL'/)) ! ! Neutral fields on geographic grid are prepared from waccm (z,u,v,w,t,etc), ! and electric field on magnetic grid from edynamo (ed1,ed2). ! This involves regridding, interpolation in height, longitude phase-shifting, ! and dimension/coordinate reversal and sign change in some cases. ! use pmgrid, only: & plat, plon ! waccm global grid dimensions use edyn_mpi, only: & mlat0,mlat1,mlon0,mlon1 ! edynamo mag subdomain first,last indices use edyn_maggrid, only: & nmlat,nmlon,nmlonp1, & ! edynamo global mag dimensions gmlat,gmlon ! edynamo mag coordinates use edyn_init,only: & use_edyn_maggrid=>time3d_use_edyn_maggrid use time3d_geogrid,only: & nlat,nlon,nlev, & ! 3d time3d global grid dims alt, & ! time3d heights (nlev) t3d_londeg=>londeg, & ! time3d global longitudes (deg) (nlon+1 includes periodic point) t3d_latdeg_n2s=>latdeg ! time3d global latitudes (deg) (north to south) use time3d_grid,only: & gmlat_t3d,gmlon_t3d ! time3d mag grid coordinates use edynamo, only: & ed1,ed2, & ! electric field on mag subdomains ed1_glb,ed2_glb ! global efield for input to time3d (nmlat,nmlonp1) use edyn_init, only: & lonshift_global, & ! routine to phase shift global longitude data use_time3d_gcmsim ! if true, call time3d_gcmsim, otherwise use waccm fields for time3d input use time3d_flds, only: & ! (time3d inputs: nlat,nlon,nlev) Tn,Wu,Ws,We,cO2,cO,cN2, & cN,cO2p,cNOp,cN2p, & ! ions Op_waccm, & ! O+ interpolated to waccm grid Ped_waccm,Hall_waccm, & ! conductivities interpolated to waccm grid Z_waccm=>Z, & ! geopotential interpolated to waccm grid Qop,Qep, & ! O+ and Photo-Electron production rates OpL ! O+ loss rate use commap, only: & ! waccm grid info latdeg, & ! global latitude coordinate (plat) londeg1=>londeg ! global longitude coordinate (plon,plat) use htinterp,only: & cuthtint ! 2d height interpolation routine ! ! Args: ! integer,intent(in) :: i0,i1,j0,j1,plev ! waccm subdomains real(r8),dimension(i0:i1,j0:j1,plev),intent(in) :: & ! waccm fields on subdomains z,t,u,v,w,o2,o1,n2,n1, & ! dynamics and composition o2p,nop,n2p,op, & ! ions sigma_ped,sigma_hall, & ! conductivities qepp,qopp, & ! photo-electron and O+ production rates lop ! O+ loss rate integer,intent(in) :: nflds character(len=*),intent(in) :: fnames(nflds) ! ! Local: integer :: i,j,k,n real(r8) :: londeg(plon) ! global waccm longitudes (deg) integer :: nlon_out(nlat) ! number of t3d lons per lat for bilin integer :: nmlon_out(nmlat) ! number of mag lons per mag lat real(r8) :: t3d_londeg_out(nlon,nlat) ! 2d t3d lons for bilin call real(r8) :: t3d_latdeg(nlat) ! time3d lats (south to north) for bilin call real(r8),dimension(plon,plat,plev) :: & ! 3d fields on global waccm grid z_glb,t_glb,u_glb,v_glb,w_glb,o2_glb, & o_glb,n2_glb,n_glb,o2p_glb,nop_glb,n2p_glb,op_glb,& Ped_glb,Hall_glb,qep_glb,qop_glb,lop_glb real(r8) :: f3d_waccm_ikj(plon,plev,plat) ! 3d global waccm field for bilin input real(r8) :: f3d_time3d_ikj(nlon,plev,nlat) ! 3d global time3d field for bilin output real(r8) :: f2d_time3d_ij(nlon,nlat) ! temporary holding array for 2d fields real(r8) :: f2d_time3d_ik(nlon,nlev) ! cuthtint output real(r8) :: f2d_time3d_z(nlon,plev) ! horizontally interpolated time3d geopotential integer :: logint,extrp,ier real(r8),dimension(nmlonp1,nmlat,1) :: & ed1_glb_ij,ed2_glb_ij, & ! efield on edynamo mag grid (irregular in lat) ed1_t3d_ij,ed2_t3d_ij ! efield on time3d mag grid (regular in lat) real(r8) :: gmlon_shifted(nmlon) ! edynamo lon coord shifted to 0->360 real(r8) :: gmlon_t3d_out(nmlon,nmlat) ! 2d time3d lons for bilin call nstep = get_nstep() ! write(iulog,"('Enter prep_time3d_input: nstep=',i4)") nstep if (use_time3d_gcmsim) goto 100 ! time3d not using waccm inputs ! ! Global waccm coordinates londeg,latdeg are from module commap: ! (latdeg is use-associated from commap) ! (londeg is local, londeg1 is use-associated from commap (renamed from londeg)) ! londeg(:) = londeg1(:,1) ! assume waccm lons are redundant in latitude ! write(iulog,"('prep_time3d_input: plon=',i4,' londeg=',/,(8f9.2))") plon,londeg ! write(iulog,"('prep_time3d_input: plat=',i4,' latdeg=',/,(8f9.2))") plat,latdeg ! ! Global time3d coordinates for sub bilin (see sub pregrid in time3d_gcmin.f): ! t3d_londeg is dimensioned (nlon+1) to include periodic point ! t3d_londeg_out is dimensioned (nlon,nlat) to exclude periodic point in bilin call ! do j=1,nlat t3d_londeg_out(:,j) = t3d_londeg(1:nlon) nlon_out(j) = nlon enddo do j=1,nlat ! reverse lat order from north-to-south to south-to-north. ! coord arrays must be monotonically increasing for sub bilin ! t3d_latdeg is local, so latdeg in time3d is unchanged t3d_latdeg(j) = t3d_latdeg_n2s(nlat-j+1) enddo ! write(iulog,"('prep_time3d_input: time3d nlon=',i5,' t3d_londeg_out(:,1)=',/,(8f10.3))") & ! nlon,t3d_londeg_out(:,1) ! write(iulog,"('prep_time3d_input: time3d nlat=',i5,' t3d_latdeg=',/,(8f10.3))") & ! nlat,t3d_latdeg(:) ! ! Gather waccm neutral winds to masterproc: ! call gather_field(z,i0,i1, j0,j1, z_glb ,plon,plat,plev,.true.,0) call gather_field(t,i0,i1, j0,j1, t_glb ,plon,plat,plev,.true.,0) call gather_field(u,i0,i1, j0,j1, u_glb ,plon,plat,plev,.true.,0) call gather_field(v,i0,i1, j0,j1, v_glb ,plon,plat,plev,.true.,0) call gather_field(w,i0,i1, j0,j1, w_glb ,plon,plat,plev,.true.,0) call gather_field(o2,i0,i1,j0,j1, o2_glb ,plon,plat,plev,.true.,0) call gather_field(o1,i0,i1,j0,j1, o_glb ,plon,plat,plev,.true.,0) call gather_field(n2,i0,i1,j0,j1, n2_glb ,plon,plat,plev,.true.,0) call gather_field(n1,i0,i1,j0,j1, n_glb ,plon,plat,plev,.true.,0) call gather_field(o2p,i0,i1,j0,j1,o2p_glb ,plon,plat,plev,.true.,0) call gather_field(nop,i0,i1,j0,j1,nop_glb ,plon,plat,plev,.true.,0) call gather_field(n2p,i0,i1,j0,j1,n2p_glb ,plon,plat,plev,.true.,0) call gather_field(op ,i0,i1,j0,j1, op_glb ,plon,plat,plev,.true.,0) call gather_field(lop,i0,i1,j0,j1,lop_glb ,plon,plat,plev,.true.,0) call gather_field(sigma_ped ,i0,i1,j0,j1,Ped_glb ,plon,plat,plev,.true.,0) call gather_field(sigma_hall,i0,i1,j0,j1,Hall_glb,plon,plat,plev,.true.,0) call gather_field(qepp,i0,i1,j0,j1,qep_glb,plon,plat,plev,.true.,0) call gather_field(qopp,i0,i1,j0,j1,qop_glb,plon,plat,plev,.true.,0) if (masterproc) then ! only masterproc will have 3d global fields from both models ! write(iulog,"('prep_time3d_input: plon=',i5,' plat=',i5)") plon,plat ! write(iulog,"(' waccm z_glb min,max=',2e12.4)") minval(z_glb),maxval(z_glb) ! write(iulog,"(' waccm u_glb min,max=',2e12.4)") minval(u_glb),maxval(u_glb) ! write(iulog,"(' waccm v_glb min,max=',2e12.4)") minval(v_glb),maxval(v_glb) ! write(iulog,"(' waccm w_glb min,max=',2e12.4)") minval(w_glb),maxval(w_glb) ! ! Save global waccm fields to nc file (these will be "top2bot"): ! write(iulog,"('prep_time3d_input: nstep=',i4,' nstep_savefld_t3d=',i4,' time2write=',l1)") & ! nstep,nstep_savefld_t3d,time2write(nstep,nstep_savefld_t3d) if (time2write(nstep,nstep_savefld_t3d)) then filename = 'Begin_prep_time3d_input.nc' action = 'append' if (nstep==1) action = 'create' call savefld_t3d(filename,'Op','WACCM O+',' ',& op_glb,(/plon,plat,plev/),(/'nlon','nlat','nlev'/),action,nstep) action = 'append' ! write(iulog,"('prep_time3d_input saving to ',a,' nstep=',i4,' action=',a)") & ! trim(filename),nstep,action ! call savefld_t3d(filename,'T','WACCM T','m/s',& ! t_glb,(/plon,plat,plev/),(/'nlon','nlat','nlev'/),action,nstep) ! call savefld_t3d(filename,'U','WACCM U','m/s',& ! u_glb,(/plon,plat,plev/),(/'nlon','nlat','nlev'/),action,nstep) ! call savefld_t3d(filename,'V','WACCM V','m/s',& ! v_glb,(/plon,plat,plev/),(/'nlon','nlat','nlev'/),action,nstep) ! call savefld_t3d(filename,'W','WACCM W','m/s',& ! w_glb,(/plon,plat,plev/),(/'nlon','nlat','nlev'/),action,nstep) ! call savefld_t3d(filename,'Z','WACCM Z','m',& ! z_glb,(/plon,plat,plev/),(/'nlon','nlat','nlev'/),action,nstep) call savefld_t3d(filename,'Qop','WACCM O+ Production',' ',& qop_glb,(/plon,plat,plev/),(/'nlon','nlat','nlev'/),action,nstep) call savefld_t3d(filename,'Qep','WACCM PE Production',' ',& qep_glb,(/plon,plat,plev/),(/'nlon','nlat','nlev'/),action,nstep) call savefld_t3d(filename,'lop','WACCM O+ Loss',' ',& lop_glb,(/plon,plat,plev/),(/'nlon','nlat','nlev'/),action,nstep) call savefld_t3d(filename,'PED','WACCM Pedersen Conductivity',' ',& ped_glb,(/plon,plat,plev/),(/'nlon','nlat','nlev'/),action,nstep) call savefld_t3d(filename,'HALL','WACCM Hall Conductivity',' ',& hall_glb,(/plon,plat,plev/),(/'nlon','nlat','nlev'/),action,nstep) endif ! ! Prepare waccm field for input to bilin regrid: ! Note reversal of levels order from waccm (alt decreasing from top at k=1), ! to time3d (alt increasing from bottom at k=1) ! Latitude coord is south-to-north (same as waccm), because bilin requires ! monotonically increasing coords. The field is switched to north-to-south ! for time3d after bilin call below. ! do n=1,nflds f3d_waccm_ikj = 0. ! init bilin input field select case(fnames(n)) case('z') do k=1,plev do j=1,plat f3d_waccm_ikj(:,k,j) = z_glb(:,j,plev-k+1)*1.e-3 ! m->km for time3d enddo enddo ! k=1,plev case('t') do k=1,plev do j=1,plat f3d_waccm_ikj(:,k,j) = t_glb(:,j,plev-k+1) enddo enddo ! k=1,plev case('u') do k=1,plev do j=1,plat f3d_waccm_ikj(:,k,j) = u_glb(:,j,plev-k+1) enddo enddo ! k=1,plev case('v') do k=1,plev do j=1,plat f3d_waccm_ikj(:,k,j) = v_glb(:,j,plev-k+1) enddo enddo ! k=1,plev case('w') do k=1,plev do j=1,plat f3d_waccm_ikj(:,k,j) = w_glb(:,j,plev-k+1) enddo enddo ! k=1,plev case('o2') do k=1,plev do j=1,plat f3d_waccm_ikj(:,k,j) = o2_glb(:,j,plev-k+1) enddo enddo ! k=1,plev case('o') do k=1,plev do j=1,plat f3d_waccm_ikj(:,k,j) = o_glb(:,j,plev-k+1) enddo enddo ! k=1,plev case('n2') do k=1,plev do j=1,plat f3d_waccm_ikj(:,k,j) = n2_glb(:,j,plev-k+1) enddo enddo ! k=1,plev case('n') do k=1,plev do j=1,plat f3d_waccm_ikj(:,k,j) = n_glb(:,j,plev-k+1) enddo enddo ! k=1,plev case('o2p') do k=1,plev do j=1,plat f3d_waccm_ikj(:,k,j) = o2p_glb(:,j,plev-k+1) enddo enddo ! k=1,plev case('nop') do k=1,plev do j=1,plat f3d_waccm_ikj(:,k,j) = nop_glb(:,j,plev-k+1) enddo enddo ! k=1,plev case('n2p') do k=1,plev do j=1,plat f3d_waccm_ikj(:,k,j) = n2p_glb(:,j,plev-k+1) enddo enddo ! k=1,plev case('op') do k=1,plev do j=1,plat f3d_waccm_ikj(:,k,j) = op_glb(:,j,plev-k+1) enddo enddo ! k=1,plev case('Ped') do k=1,plev do j=1,plat f3d_waccm_ikj(:,k,j) = Ped_glb(:,j,plev-k+1) enddo enddo ! k=1,plev case('Hall') do k=1,plev do j=1,plat f3d_waccm_ikj(:,k,j) = Hall_glb(:,j,plev-k+1) enddo enddo ! k=1,plev case('Qop') do k=1,plev do j=1,plat f3d_waccm_ikj(:,k,j) = qop_glb(:,j,plev-k+1) enddo enddo ! k=1,plev case('Qep') do k=1,plev do j=1,plat f3d_waccm_ikj(:,k,j) = qep_glb(:,j,plev-k+1) enddo enddo ! k=1,plev case('OpL') do k=1,plev do j=1,plat f3d_waccm_ikj(:,k,j) = lop_glb(:,j,plev-k+1) enddo enddo ! k=1,plev case default write(iulog,"('>>> prep_time3d_input 1: unknown fname=',a)") fnames(n) call endrun end select ! ! Bilinear interpolation from waccm grid to time3d grid, at plev levels: f3d_time3d_ikj = 0. ! init bilin output field on time3d grid (plev levels) do k=1,plev call bilin(f3d_waccm_ikj(:,k,:),londeg,latdeg,plon,plon,1,1,plat, & f3d_time3d_ikj(:,k,:),t3d_londeg_out,t3d_latdeg,nlon,nlon_out,1,nlat) ! ! t3d_latdeg and field is south-to-north because sub bilin requires monotonically ! increasing coords, but time3d expects fields to have lat north-to-south, so change ! to that now: f2d_time3d_ij(:,:) = f3d_time3d_ikj(:,k,:) do j=1,nlat f3d_time3d_ikj(:,k,j) = f2d_time3d_ij(:,nlat-j+1) enddo enddo ! k=1,plev ! ! Now interpolate from waccm pressure levels to time3d altitude levels. logint = 0 extrp = 0 do j=1,nlat ! ! Save waccm interpolated geopotential f2d_time3d_z ! (z must be the first field, i.e., n==1): ! (k-index already reversed above, and m->km already done above) ! ! Note where top of time3d (600 km) exceeds top of waccm (~420 km), ! assign the top waccm value to the time3d input field (see cuthtint) ! if (fnames(n)=='z') then do k=1,plev f2d_time3d_z(:,k) = f3d_time3d_ikj(:,k,j) enddo cycle ! skip further processing of Z, cycle to next field endif f2d_time3d_ik=0. ! cuthtint output ! subroutine cuthtint(fin,fht,idim1,nzp,fout,hts,nhts,logint,extrp,ier,iprnt) ! integer,intent(in) :: idim1,nzp,nhts,logint,extrp,iprnt ! real(r8),intent(in) :: fin(idim1,nzp),fht(idim1,nzp),hts(nhts) ! real(r8),intent(out) :: fout(idim1,nhts) ! real(r8) :: f3d_time3d_ikj(nlon,plev,nlat) ! 3d global time3d field for bilin output ! real(r8) :: f2d_time3d_ik(nlon,nlev) ! cuthtint output ! real(r8) :: f2d_time3d_z(nlon,plev) ! horizontally interpolated time3d geopotential call cuthtint(f3d_time3d_ikj(:,:,j),f2d_time3d_z(:,:),nlon,plev, & f2d_time3d_ik,alt,nlev,logint,extrp,ier,1) ! ! Save interpolated time3d field for input to time3d: select case(fnames(n)) ! case('t') do k=1,nlev do i=1,nlon Tn(j,i,k) = f2d_time3d_ik(i,k) enddo enddo case('u') do k=1,nlev do i=1,nlon We(j,i,k) = f2d_time3d_ik(i,k) enddo enddo case('v') ! note sign change (time3d meridional wind is positive south) do k=1,nlev do i=1,nlon Ws(j,i,k) = -f2d_time3d_ik(i,k) ! make positive south for time3d enddo enddo case('w') do k=1,nlev do i=1,nlon Wu(j,i,k) = f2d_time3d_ik(i,k) enddo enddo case('o2') do k=1,nlev do i=1,nlon cO2(j,i,k) = f2d_time3d_ik(i,k) enddo enddo case('o') do k=1,nlev do i=1,nlon cO(j,i,k) = f2d_time3d_ik(i,k) enddo enddo case('n2') do k=1,nlev do i=1,nlon cN2(j,i,k) = f2d_time3d_ik(i,k) enddo enddo case('n') do k=1,nlev do i=1,nlon cN(j,i,k) = f2d_time3d_ik(i,k) enddo enddo case('o2p') do k=1,nlev do i=1,nlon cO2p(j,i,k) = f2d_time3d_ik(i,k) enddo enddo case('nop') do k=1,nlev do i=1,nlon cNOp(j,i,k) = f2d_time3d_ik(i,k) enddo enddo case('n2p') do k=1,nlev do i=1,nlon cN2p(j,i,k) = f2d_time3d_ik(i,k) enddo enddo case('op') do k=1,nlev do i=1,nlon Op_waccm(j,i,k) = f2d_time3d_ik(i,k) enddo enddo case('Ped') do k=1,nlev do i=1,nlon Ped_waccm(j,i,k) = f2d_time3d_ik(i,k) enddo enddo case('Hall') do k=1,nlev do i=1,nlon Hall_waccm(j,i,k) = f2d_time3d_ik(i,k) enddo enddo case('Qep') do k=1,nlev do i=1,nlon Qep(j,i,k) = max(f2d_time3d_ik(i,k),1.e-10) enddo enddo case('Qop') do k=1,nlev do i=1,nlon Qop(j,i,k) = max(f2d_time3d_ik(i,k),1.e-10) enddo enddo case('OpL') do k=1,nlev do i=1,nlon OpL(j,i,k) = f2d_time3d_ik(i,k) enddo enddo case default write(iulog,"('>>> prep_time3d_input 2: unknown fname=',a)") fnames(n) call endrun end select enddo ! j=1,nlat (call cuthtint) enddo ! n=1,nflds ! ! Save interpolated fields to netcdf file: ! if (time2write(nstep,nstep_savefld_t3d)) then filename = 'End_prep_time3d_input.nc' action = 'append' if (nstep==1) action = 'create' call savefld_t3d(filename,'Op','O+ (N2S, bot2top)',& ' ',Op_waccm,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) action = 'append' call savefld_t3d(filename,'Qop','O+ Production (N2S, bot2top)',& ' ',Qop,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) call savefld_t3d(filename,'Qep','Photo-Electron PRate (N2S, bot2top)',& ' ',Qep,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) call savefld_t3d(filename,'lop','O+ Loss (N2S, bot2top)',& ' ',OpL,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) call savefld_t3d(filename,'Ped_waccm','Pedersen Conductivity (N2S, bot2top)',& ' ',Ped_waccm,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) call savefld_t3d(filename,'Hall_waccm','Hall Conductivity (N2S, bot2top)',& ' ',Hall_waccm,(/nlat,nlon,nlev/),(/'nlat','nlon','nlev'/),action,nstep) endif endif ! masterproc 100 continue ! jump to here if use_time3d_gcmsim ! ! Gather ed1,2 (edynamo output on subdomains) to the root task for input to time3d: ! ! edynamo: real(r8),allocatable,dimension(:,:),save :: ed1,ed2 ! (mlon0-1:mlon1+1,mlat0-1:mlat1+1) ! time3d: dimension(nmlat,nmlonp1) :: ed1_glb,ed2_glb ! Local: dimension(nmlonp1,nmlat) :: ed1_glb_ij,ed2_glb_ij ! call gather_field(ed1(mlon0:mlon1,mlat0:mlat1),mlon0,mlon1,mlat0,mlat1,& ed1_glb_ij,nmlonp1,nmlat,1,.false.,0) call gather_field(ed2(mlon0:mlon1,mlat0:mlat1),mlon0,mlon1,mlat0,mlat1,& ed2_glb_ij,nmlonp1,nmlat,1,.false.,0) if (masterproc) then ! ! Shift edynamo mag lon coordinate from "-180to180" to "zeroto360": ! gmlon_shifted = gmlon(1:nmlon) call lonshift_global(gmlon_shifted,nmlon,'zeroto360',.true.) ! write(iulog,"('prep_time3d: gmlon_shifted=',/,(8f9.2))") gmlon_shifted ! ! Shift edynamo global 2d ed1,2 longitudes from "-180to180" to "zeroto360": ! do j=1,nmlat call lonshift_global(ed1_glb_ij(1:nmlon,j,1),nmlon,'zeroto360',.false.) call lonshift_global(ed2_glb_ij(1:nmlon,j,1),nmlon,'zeroto360',.false.) enddo ed1_glb_ij(nmlonp1,:,1) = ed1_glb_ij(1,:,1) ed2_glb_ij(nmlonp1,:,1) = ed2_glb_ij(1,:,1) ! ! Interpolate 2d global efield from edynamo mag grid (irregular in latitude) to ! time3d mag grid (regular in latitude), if time3d is not using edyn maggrid: ! (note that both mag grids have the same dimensions) ! if (.not.use_edyn_maggrid) then do j=1,nmlat nmlon_out(j) = nmlon ! number of lons per lat for bilin gmlon_t3d_out(:,j) = gmlon_t3d(1:nmlon) ! 2d lons for bilin enddo call bilin(ed1_glb_ij(1:nmlon,:,1),gmlon_shifted,gmlat,nmlon,nmlon,1,1,nmlat,& ed1_t3d_ij(1:nmlon,:,1),gmlon_t3d_out,gmlat_t3d,nmlon,nmlon_out,1,nmlat) call bilin(ed2_glb_ij(1:nmlon,:,1),gmlon_shifted,gmlat,nmlon,nmlon,1,1,nmlat,& ed2_t3d_ij(1:nmlon,:,1),gmlon_t3d_out,gmlat_t3d,nmlon,nmlon_out,1,nmlat) ! ! Periodic point: ed1_t3d_ij(nmlonp1,:,1) = ed1_t3d_ij(1,:,1) ed2_t3d_ij(nmlonp1,:,1) = ed2_t3d_ij(1,:,1) elseif (masterproc) then ! time3d is using edynamo mag grid ed1_t3d_ij = ed1_glb_ij ed2_t3d_ij = ed2_glb_ij endif ! use_edyn_maggrid or not ! ! Transform from (nmlonp1,nmlat) to (nmlat,nmlonp1) (edynamo to time3d) ! ed1,2_glb are in edynamo module, use-associated by time3d module. ! do i=1,nmlonp1 do j=1,nmlat ed1_glb(j,i) = ed1_t3d_ij(i,j,1) ed2_glb(j,i) = ed2_t3d_ij(i,j,1) enddo enddo ed1_glb(:,nmlonp1) = ed1_t3d_ij(1,:,1) ed2_glb(:,nmlonp1) = ed2_t3d_ij(1,:,1) ! write(iulog,"('prep_time3d_input: nstep=',i4,' ed1_glb,ed2_glb=',2(1pe12.4))") & ! nstep,minval(ed1_glb),maxval(ed2_glb) endif ! masterproc end subroutine prep_time3d_input !----------------------------------------------------------------------- subroutine gather_field(fsub,ilon0,ilon1,ilat0,ilat1,fglb,nglblon,nglblat,nlev,geo,iprint) use edyn_mpi,only: mp_gather2root ! ! Gather input field fsub on subdomains to the root task in global output field fglb. ! This routine should be expanded for multiple fields. ! ! Args: integer,intent(in) :: ilon0,ilon1,ilat0,ilat1,nglblon,nglblat,nlev,iprint real(r8),intent(in) :: fsub(ilon0:ilon1,ilat0:ilat1,nlev) real(r8),intent(out) :: fglb(nglblon,nglblat,nlev) logical,intent(in) :: geo ! write(iulog,"('Enter time3d gather_field: ilon0,1=',2i4,' ilat0,1=',2i4,' nglblon,lat=',2i4,' geo=',l1,' fsub min,max=',2(1pe12.4))") & ! ilon0,ilon1,ilat0,ilat1,nglblon,nglblat,geo,minval(fsub),maxval(fsub) fglb = 0._r8 call mp_gather2root(fsub,ilon0,ilon1,ilat0,ilat1,fglb,nglblon,nglblat,nlev,geo,iprint) ! write(iulog,"('time3d gather_field after gather: fglb min,max=',2(1pe12.4))") & ! minval(fglb),maxval(fglb) end subroutine gather_field !----------------------------------------------------------------------- logical function time2write(nstep,nstep_savefld_t3d) ! ! True if global fields should be saved to nc file at nstep. ! integer,intent(in) :: nstep,nstep_savefld_t3d time2write = .false. if (nstep==1.or.mod(nstep,nstep_savefld_t3d)==0) time2write = .true. end function time2write !----------------------------------------------------------------------- end module time3d