! module ctmt_module ! ! This software is part of the NCAR TIE-GCM. Use is governed by the ! Open Source Academic Research License Agreement contained in the file ! tiegcmlicense.txt. ! 02/13 B Emery adapted from gswm.F reading Jens Oberheide's diur and sd .nc files ! ! Read CTMT boundary perturbations Z, T, U, V, and interpolate ! to model grid and time, where Z comes from the neutral density percent variation ! from Mack Jones: dZ = RTo(dRHO/RHOo) + RdT=8314*(181K(dRHO/RHOo)+dT) at 97km ! use params_module,only: nlon,nlat,nlonp4,nlonp2,nlevp1, | dlev,zibot,nlev use init_module,only: glon,glat use mpi_module,only: lon0,lon1,lat0,lat1 use nchist_module,only:nc_open,nc_close,handle_ncerr use addfld_module,only: addfld use cons_module,only: pi use input_module,only: | ctmt_nlev ! number of levels to nudge with CTMT ! if ctmt_nlev=1, only lbc; if >1, nudge ctmt_nlev-1 levels ! since at ctmt_nlev, the percentage of nudging is usually 0 implicit none ! integer,parameter :: nalt5=5 ! for outzav_04080_120f107_95to105km ! integer,parameter :: nalt5=11 ! for outzav_04080_120f107_95to120km integer,parameter :: nalt5=13 ! for outzav_04080_120f107_90to120km ! #ifdef SUN #include #else #include #endif ! ! CTMT boundary data at model grid, output by this module, ! to be used by dt, duv, addiag, etc: ! ! Private module data, read by sub rdctmt: ! am,ph_##_v(month,alt,glat)=(12,161,37) where month=1:1:12, alt=0:2.5:400, glat=-90:5:90 ! Easiest assumption for ctmt_nlev=1 only (ctmt tides at lbc only) ! Convert to 97km (use 97.5km at #40=nlb) at UT, for nlat glats and nlonp4 glons ! Interpolation from constant heights 95:2.5:105km for "nudging" on ctmt_nlev p sfcs ! lbc ~97km is ~94-99km, so interpolate 2.5km in ht for ctmt_nlev and add ! zonal means (zm) from files like ~emery/tidi/outzav_04080_120f107_95to105km integer,parameter,private :: nmonth= 12, nhour=24, nlat5=37 ! integer,parameter,private :: nalt2p5=161 ! Not needed for partial array ! dimension of 2 at front is nn=1 diurnal and nn=2 semidiurnal ! dimension at end nalt5 (=5 for 95:2.5:105) is number of hts to save ! NOTE: if ctmt_nlev so large z hts>~105km, must revise dim 5 and zm file ! ! nalt5 is number of hts to save (formerly nlb=1) real,dimension(2,8,4,nmonth,nlat5,nalt5),private :: amplb,phlb real,dimension(2,8),private :: ! Set in rdctmt and used in ctmtdat | fac(2,8), ! amplitude mulipliers for 1,2 di,sd and 8 cmp components | phsft(2,8) ! phase shifts in hours for 1,2 di,sd and 8 cmp real,dimension(nlat),private :: faclat1,faclat2 integer,dimension(nlat),private :: jl1 real,dimension(nlonp4),private :: glonp4 ! real,dimension(nalt5) :: ctmt_alt real,dimension(13) :: ctmt_alt real,dimension(72,nalt5) :: zm_t,zm_d,zm_u,zm_v ! zonal mean tn,un,vn,z for TIEGCM (was not calc correctly in lbc.F!) ! real,dimension(nlev,nlat) :: tzm_tie,uzm_tie,vzm_tie,zzm_tie ! allocated subdomains: ! Allocated by sub alloc_ctmt (called by allocdata): real,allocatable,dimension(:,:,:) :: ! (lon0:lon1,lat0:lat1,nalt5) | ctmt_t, ctmt_u, ctmt_v, ctmt_d ! Percentage of CTMT+zm data to nudge with from lbc to ctmt_nlev ! real,dimension(nlevp1) :: ctmt_nudge ! only defined to ctmt_nlev real,allocatable,save :: ctmt_nudge(:) ! dimensioned later to ctmt_nlev ! Arrays on constant p levels for CTMT plus zonal means (zm) (for lbc,dt,duv.F) real,allocatable,dimension(:,:,:) :: ! (lon0:lon1,lat0:lat1,ctmt_nlev) | t_zmctmt, u_zmctmt, v_zmctmt ! contains !----------------------------------------------------------------------- subroutine getctmt (itp) ! ! Module driver to read nc files, and do time interpolations. use init_module,only: istep,iyear,iday,secs ! ! Files provided by user via namelist read where have 6 di and 8 sdi tides, ! where s0 is migrating tide, nm di tides are w1,w2,e1,e2,e3, and ! nm sdi tides are w1,w2,w3,w4,e1,e2,e3,e4 where phases and amps can be revised use input_module,only: | ctmt_di_ncfile,ctmt_sdi_ncfile,zm_file ! if ctmt_nlev>0, must read a zonal mean (zm) file for the day and the F107 created ! from ~emery/msis00/msis2tie.f every 2.5km from 95:2.5:105km (5 alts) ! and from ~emery/tidi/meankm.m to create zonal mean Tn, %nden/nden, Un, Vn ! on glat dres grid of TIEGCM (72 glats -88.75:2.5:88.75) to use in 2.5 or 5deg ! Sample output: ~emery/tidi/outzav_04080_120f107_95to105km (or 95to120km) ! interface p levels for Z,DEN and midpoint p levels for TN,UN,VN ! use fields_module,only: z, ! z(nlev,nlon,nlat,2) usually allocated ! | tn,un,vn ! integer,external :: nextlu ! ! Driver for obtaining CTMT data, called once per timestep from lbc. ! ! Args: integer,intent(in) :: itp ! Local: integer :: iprint,nalt real,dimension(nlon,nlat) :: utid integer :: il,jl,nn,lu_zm,nd_ctmt real :: globtn,globden,globun,globvn,test character(len=80) :: zm_char real,dimension(72) :: glat72 integer :: istat,k real :: zp(nlevp1) ! iprint = 0 if (istep==1) iprint = 1 ! ! Get data for 6 CTMT diurnal tides and 8 CTMT semi-diurnal tides: if (istep==1) then ! Read diurnal and semidiurnal CTMT files if (len_trim(ctmt_di_ncfile)>0) | call rdctmt(ctmt_di_ncfile,'diur') if (len_trim(ctmt_sdi_ncfile)>0) | call rdctmt(ctmt_sdi_ncfile,'semi') endif ! if (istep==1) then call mkctmt(iyear,iday,int(secs),iprint) ! Calc zonal means for tn,un,vn,z for ctmt_nlev+1 p levels ! do k=1,ctmt_nlev+1 ! do jl=1,nlat ! tzm_tie(k,jl) = 0. ! uzm_tie(k,jl) = 0. ! vzm_tie(k,jl) = 0. ! zzm_tie(k,jl) = 0. ! do il=3,nlon+2 ! tzm_tie(k,jl) = tzm_tie(k,jl) + tn(k,il,jl,itp) ! uzm_tie(k,jl) = uzm_tie(k,jl) + un(k,il,jl,itp) ! vzm_tie(k,jl) = vzm_tie(k,jl) + vn(k,il,jl,itp) ! zzm_tie(k,jl) = zzm_tie(k,jl) + z(k,il,jl,itp) ! enddo ! tzm_tie(k,jl) = tzm_tie(k,jl)/nlon ! uzm_tie(k,jl) = uzm_tie(k,jl)/nlon ! vzm_tie(k,jl) = vzm_tie(k,jl)/nlon ! zzm_tie(k,jl) = zzm_tie(k,jl)/nlon ! enddo ! if (istep<3) ! | write (6,"('zm at nlev min t,u,v,z =',i2,8e12.5)") k, ! | minval(tzm_tie(k,:)),minval(uzm_tie(k,:)),minval(vzm_tie(k,:)), ! | minval(zzm_tie(k,:)),minval(tn(k,:,:,itp)), ! | minval(un(k,:,:,itp)),minval(vn(k,:,:,itp)),minval(z(k,:,:,itp)) ! enddo ! do k=1,ctmt_nlev+1 end subroutine getctmt !----------------------------------------------------------------------- subroutine rdctmt(ncfile,type) use input_module,only: mxlen_filename implicit none ! ! Args: character(len=*),intent(in) :: ncfile,type ! ! Local: integer :: ncid,istat character(len=mxlen_filename) :: dskfile integer :: nlat_rd, nmonth_rd integer :: id_nalt, id_nmonth, id_nlat, id_am, id_ph character(len=240) :: char240 real,dimension(nlat5) :: ctmt_lat ! .nc files are listed in C array order or amp(mon,alt,lat) but ncdump -bF puts ! it in fortran array order or amp(lat,alt,mon) which are dims of amptmp and phtmp! ! real,dimension(nlat5,nalt2p5,nmonth) :: amptmp,phtmp ! Only need partial slab of altitudes, say 95-120km (39-49) or 11 of 161 hts real,dimension(nlat5,13,nmonth) :: amptmp,phtmp integer :: i,im,il,jl,kl character(len=16) :: cmp character(len=4) :: var character(len=8) :: ampnm character(len=10) :: phnm ! ! nalt5 is number of hts to save (formerly nlb=1) real :: ctmt_all_alts(161),bot_alt integer :: mo,nct,nn,nc,nv,nc1,nh,nlb,nalt,mlb integer,dimension(3) :: istart,icount data cmp/'e3e2e1s0w1w2w3w4'/ ! 8 tides for sd (first 6 for diurnal) data var/'dtuv'/ ! Use 4 of 5 CTMT variables (skip 'w') ! NOTE: density is NOT in percent here, but is in ratio d_pert/d_bkg or [-] ! Set amplitude factors for CTMT diurnal and semi-diurnal tides (default=1) ! and set phase shifts (default=0 h, + are shifted earlier) do nn=1,2 do nc=1,8 fac(nn,nc) = 1. ! 8/14: Test with only 1 component ! fac(nn,nc) = 0. phsft(nn,nc) = 0. enddo enddo ! 8/14: Test with only 1 component ! fac(1,5) = 1. ! 3/12: migrating DW1*2. shift 5h earlier, migrating SW2 shift 5h earlier ! fac(1,5) = 2. ! phsft(1,5) = 5. ! phsft(2,6) = 5. ! 3/12: migr SW2*2 w 3nudge_rev 1,.6,.15 since DW1*2 gave NEG day Viz! This -+-+(PRE) ! fac(2,6) = 2. ! d -32% to +29%, T=143-227K, U=-73to+63, V=-103to+91, z=94.4-98.3km ! Pick partial array of 11 heights from 95-120km (#39-49) of 161 hts total 0:2.5:400km ! amptmp(37lat,161ht,12mo) nlb=40 (97.5km) for full array, nlb=2 for partial array istart(1) = 1 istart(2) = 39 if (nalt5==13) istart(2) = 37 istart(3) = 1 icount(1) = nlat5 icount(2) = nalt5 icount(3) = nmonth ! Set nlb for 97.5km (nlb not used now, and save 95:2.5:105km) nlb = 2 ! TEMP use full array for testing! !! nlb = 40 nn = 0 if (trim(type)=='diur') then nn = 1 nct = 6 endif if (trim(type)=='semi') then nn = 2 nct = 8 endif if (nn==0) then write(6,"(/,'>>> rdctmt: unknown type=',a)") type call shutdown('rdctmt') endif dskfile = ' ' call getfile(ncfile,dskfile) write(6,"(/,72('-'))") write(6,"('Reading CTMT data file ',a)") trim(ncfile) ! ! Open netcdf file: ! am,ph_##_v(month,alt,glat)=(12,161,37) where month=1:1:12, alt=0:2.5:400, glat=-90:5:90 call nc_open(ncid,dskfile,'OLD','READ') if (ncid==0) then write(6,"(/,'>>> rdctmt: error opening netcdf CTMT ', | 'file ',a,' dskfile ',a)") trim(ncfile),trim(dskfile) call shutdown('rdctmt') endif ! Get partial altitude grid ! istat = nf_inq_dimid(ncid,'altitude',id_nalt) ! CANNOT BE FOUND?? ! istat = nf_inq_dimid(ncid,'alt',id_nalt) ! istat = nf_get_var_double(ncid,id_nalt,ctmt_all_alts) ! if (istat /= NF_NOERR) call handle_ncerr(istat, ! | 'rdctmt: Error getting variable ctmt_all_alts') ! do i=istart(2),istart(2)+icount(2)-1 ! ctmt_alt(i-istart(2)+1) = ctmt_all_alts(i) ! enddo ! 6/13: Give up and will write 13 altitudes from 90-105km or 11 from 95-105 bot_alt = 95. if (nalt5==13) bot_alt = 90. do i=1,nalt5 ctmt_alt(i) = bot_alt + (i-1)*2.5 enddo write (6,"(1x,'from istart for count: ctmt altitude =',2i3, | 15f6.1)") istart(2),icount(2),(ctmt_alt(i),i=1,icount(2)) ! ! Check nmonth dimension: istat = nf_inq_dimid(ncid,'month',id_nmonth) istat = nf_inq_dimlen(ncid,id_nmonth,nmonth_rd) if (istat /= NF_NOERR) then write(char240,"('rdctmt: Error getting nmonth dimension from ', | 'file ',a)") trim(ncfile) call handle_ncerr(istat,char240) endif if (nmonth_rd /= nmonth) then write(6,"(/,'>>> rdctmt: nmonth_rd=',i4,' not equal to nmonth=', | i4)") nmonth_rd,nmonth write(6,"('ctmt data file: ',a)") trim(ncfile) call shutdown('rdctmt') endif ! ! Get nlat5 dimension (37 values of lat=glat_ctmt=-90:5:90): istat = nf_inq_dimid(ncid,'lat',id_nlat) istat = nf_inq_dimlen(ncid,id_nlat,nlat_rd) if (istat /= NF_NOERR) then write(char240,"('rdctmt: Error getting nlat dimension from ', | 'file ',a)") trim(ncfile) call handle_ncerr(istat,char240) endif if (nlat_rd /= nlat5) then write(6,"(/,'>>> rdctmt: nlat_rd=',i4,' not equal to nlat5=', | i4)") nlat_rd,nlat5 write(6,"('CTMT data file: ',a)") trim(ncfile) call shutdown('rdctmt') endif ! Get CTMT glat ctmt_lat - should be -90:5:90 istat = nf_get_var_double(ncid,id_nlat,ctmt_lat) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdctmt: Error getting variable ctmt_lat') ! Find interpolation factors for lat do jl=1,nlat do kl=1,nlat5-1 if (glat(jl) .ge. ctmt_lat(kl) .and. | glat(jl) .le. ctmt_lat(kl+1)) then faclat1(jl) = (ctmt_lat(kl+1)-glat(jl)) / | (ctmt_lat(kl+1)-ctmt_lat(kl)) faclat2(jl) = 1.-faclat1(jl) jl1(jl) = kl endif enddo enddo ! TEMP ! write (6,"(1x,'jl1 and faclat1 between ctmt_lat and glat')") ! write (6,"(1x,10(i3,f7.2))") (jl1(jl),faclat1(jl),jl=1,nlat) ! write (6,"(1x,10f7.2)") (ctmt_lat(jl),jl=1,nlat5) ! ! Read amplitudes and phases (month,alt,glat)=(12,161,37) ! where month=1:1:12, alt=0:2.5:400, glat=-90:5:90; Pick alt=#40 at 97.5km ! Migrating tide w1 for diurnal and w2 for semi-diurnal do nv=1,4 do nc=1,nct nc1 = (nc-1)*2 + 1 ampnm = 'amp_'//cmp(nc1:nc1+1)//'_'//var(nv:nv) phnm = 'phase_'//cmp(nc1:nc1+1)//'_'//var(nv:nv) istat = nf_inq_varid(ncid,ampnm,id_am) if (istat /= NF_NOERR) then write (6,"(1x,'nn nc nv =',3i2)") nn,nc,nv call handle_ncerr(istat,'rdctmt: Error getting amp id') endif ! istat = nf_get_var_double(ncid,id_am,amptmp) ! get full array ! Pick partial array of 11 heights from 95-120km (#39-49) of 161 hts total 0:2.5:400km istat = nf_get_vara_double(ncid,id_am,istart,icount,amptmp) if (istat /= NF_NOERR) then write (6,"(1x,'nn nc nv =',3i2)") nn,nc,nv call handle_ncerr(istat,'rdctmt: Error getting amp') endif ! amplb(nn,nc,nv,:,:)=amptmp(:,nlb,:) ! Now have wrong array order istat = nf_inq_varid(ncid,phnm,id_ph) if (istat /= NF_NOERR) then write (6,"(1x,'nn nc nv =',3i2)") nn,nc,nv call handle_ncerr(istat,'rdctmt: Error getting phase id') endif ! istat = nf_get_var_double(ncid,id_ph,phtmp) ! get full array istat = nf_get_vara_double(ncid,id_ph,istart,icount,phtmp) if (istat /= NF_NOERR) then write (6,"(1x,'nn nc nv =',3i2)") nn,nc,nv call handle_ncerr(istat,'rdctmt: Error getting phase') endif ! phlb(nn,nc,nv,:,:)=phtmp(:,nlb,:) ! Now have wrong array order do jl=1,nlat5 do mo=1,12 ! Save at 97.5km ! amplb(nn,nc,nv,mo,jl,1)=amptmp(jl,nlb,mo) ! phlb(nn,nc,nv,mo,jl,1)=phtmp(jl,nlb,mo) ! TEMP save 105,112.5,120km ! amplb(nn,nc,nv,mo,jl,2)=amptmp(jl,nlb+3,mo) ! phlb(nn,nc,nv,mo,jl,2)=phtmp(jl,nlb+3,mo) ! amplb(nn,nc,nv,mo,jl,3)=amptmp(jl,nlb+6,mo) ! phlb(nn,nc,nv,mo,jl,3)=phtmp(jl,nlb+6,mo) ! amplb(nn,nc,nv,mo,jl,4)=amptmp(jl,nlb+9,mo) ! phlb(nn,nc,nv,mo,jl,4)=phtmp(jl,nlb+9,mo) ! Now save nalt5 hts do mlb=1,nalt5 ! TEMP: Print out amp/ph at mo=3,jl=19(eq) for nalt5 hts from istart to icount hts ! if (mo==3 .and. jl==19) ! | write (6,"(1x,'Mar 0glat: nn nc nv amp ph ht=',4i3,2f8.2)") ! | nn,nc,nv,mlb,amptmp(jl,mlb,mo),phtmp(jl,mlb,mo) amplb(nn,nc,nv,mo,jl,mlb)=amptmp(jl,mlb,mo) phlb(nn,nc,nv,mo,jl,mlb)=phtmp(jl,mlb,mo) enddo enddo enddo ! TEMP ! write (6,"(1x,4i2,1x,a8,1x,a10,4e12.4,' nlb nn nc nv ampnm ', ! | 'phnm (11b1)a,p')") nlb,nn,nc,nv,ampnm,phnm,amptmp(1,1,1), ! | amptmp(1,nlb,1),phtmp(1,1,1),phtmp(1,nlb,1) ! if (nn==2 .and. nv==4 .and. nc==4) then ! write (6,"(1x,4i2,1x,a8,1x,a10' nlb nn nc nv ampnm phnm')") ! | nlb,nn,nc,nv,ampnm,phnm ! write (6,"(1x,'Sep 95-105km (#1-5) 60S-60N (#7-31)')") ! do jl=7,31 ! write (6,"(f6.1,5f7.2)") ctmt_lat(jl),(amptmp(jl,nh,9),nh=1,5) ! enddo ! write (6,"(1x,i1,'=nn Sep s0 v at 97.5km w lat: amp ph')")nn ! write (6,"(1x,10f7.2)") (amptmp(jl,nlb,9),jl=1,nlat5) ! write (6,"(1x,10f7.2)") (phtmp(jl,nlb,9),jl=1,nlat5) ! endif ! if (nn==1 .and. nv==1 .and. nc==2) then ! write (6,"(1x,i1,'=nn DE2 den #157 390km w mo w lat: amp ph')")nn ! write (6,"(12f7.2)")((amptmp(jl,157,mo)*100.,mo=1,12),jl=1,nlat5) ! write (6,"(12f7.2)")((phtmp(jl,157,mo),mo=1,12),jl=1,nlat5) ! endif ! if (nn==1 .and. nv==1 .and. nc==1) then ! write (6,"(1x,3i2,1x,a8,1x,a10' nn nc nv ampnm phnm')") ! | nn,nc,nv,ampnm,phnm ! write (6,"(1x,'amp den DE3 mo=1-12, nh=40=nlb, jl=19 equ')") ! write (6,"(12f7.4)")(amptmp(19,nlb,mo)*100.,mo=1,12) ! write (6,"(1x,i1,'=nn DE3 den #157 390km w mo w lat: amp ph')")nn ! write (6,"(12f7.2)")((amptmp(jl,157,mo)*100.,mo=1,12),jl=1,nlat5) ! write (6,"(12f7.2)")((phtmp(jl,157,mo),mo=1,12),jl=1,nlat5) ! endif ! if (nn==2 .and. nv==1 .and. nc==4) then ! write (6,"(1x,i1,'=nn S0 den #157 390km w mo w lat: amp ph')")nn ! write (6,"(12f7.2)")((amptmp(jl,157,mo)*100.,mo=1,12),jl=1,nlat5) ! write (6,"(12f7.2)")((phtmp(jl,157,mo),mo=1,12),jl=1,nlat5) ! endif ! if (nn==2 .and. nv==1 .and. nc==2) then ! write (6,"(1x,i1,'=nn SE2 den #157 390km w mo w lat: amp ph')")nn ! write (6,"(12f7.2)")((amptmp(jl,157,mo)*100.,mo=1,12),jl=1,nlat5) ! write (6,"(12f7.2)")((phtmp(jl,157,mo),mo=1,12),jl=1,nlat5) ! endif enddo ! do nc=1,nct enddo ! do nv=1,4 ! Set glonp4 or glon plus wrap-around points do il=1,nlon glonp4(il+2) = glon(il) enddo glonp4(1) = glonp4(nlon+1) glonp4(2) = glonp4(nlon+2) glonp4(nlon+3) = glonp4(3) glonp4(nlon+4) = glonp4(4) ! ! Close the file: call nc_close(ncid) write(6,"('Completed read from CTMT data file ',a)") trim(ncfile) write(6,"(/,72('-'))") end subroutine rdctmt !----------------------------------------------------------------------- subroutine mkctmt(iyear,iday,isecs,iprint) ! ! Use data read from ctmt file to return perturbations in d,t,u,v ! at current model date and time. It is assumed that the ctmt file ! provides data for a whole year with one day of data for each month ! with hourly values (0 UT to 23 UT). ! The data is linearly interpolated first to the model UT for the ! corresponding months (current and next month) and then linearly ! interpolated to the modelday ! use hist_module,only: modeltime #ifdef MPI use mpi_module,only: mp_periodic_f2d #endif implicit none ! ! Args: integer,intent(in) :: iyear,iday,isecs,iprint ! ! Local: real :: difsec,difday,sec,secint,dayint integer :: i,j,ndmon_nl(13),ndmon_lp(13),nonleap, | mon_cur,mon_nxt,ihr_cur,ihr_nxt,nointp=0, | lonbeg,lonend,i0,i1,j0,j1,malt,nalt real,dimension(lon0:lon1,lat0:lat1) :: t,u,v,d, | d_curmo, t_curmo, u_curmo, v_curmo, | d_nxtmo, t_nxtmo, u_nxtmo, v_nxtmo, | d_curmo_nxthr, t_curmo_nxthr, u_curmo_nxthr, v_curmo_nxthr, | d_nxtmo_nxthr, t_nxtmo_nxthr, u_nxtmo_nxthr, v_nxtmo_nxthr real :: fctmt(lon0:lon1,lat0:lat1,4) ! for mp calls ! ! Data: CTMT data given at the middle of each month ! -> day of year for this data assuming non-leap year ! J F M A M J J A S O N D J data ndmon_nl ! non_leap | /15,46,74,105,135,166,196,227,258,288,319,349,380/ ! ! External: real,external :: finterp ! ! For addfld calls: i0 = lon0 ; i1 = lon1 j0 = lat0 ; j1 = lat1 ! ! Get month of model run do i = 1,13 if(iday.le.ndmon_nl(i)) goto 10 enddo 10 mon_nxt = i ! next month if(mon_nxt == 13 ) mon_nxt = 1 mon_cur = i-1 ! current month if(mon_cur == 0 ) mon_cur = 12 ! ! Get hour of model run (model hours from 0 to 23 ) ihr_cur = modeltime(2) ! current hour ! ihr_nxt = modeltime(2)+1 ! next hour ! ihr_cur = modeltime(2)+2 ! current hour - shift 2 hr earlier ALL components! if(ihr_cur >=24) ihr_cur=ihr_cur-24 ihr_nxt = ihr_cur+1 ! next hour if(ihr_nxt == 24 ) ihr_nxt = 1 ! ! Calculate nalt5 altitudes (95:2.5:105km) do nalt=1,nalt5 ! ! Subdomains at current month, current hour: d_curmo = ctmtdat(mon_cur,ihr_cur+1,'d',nalt) t_curmo = ctmtdat(mon_cur,ihr_cur+1,'t',nalt) ! TEMP write (6,"(1x,'nalt t_curmo(lon0,lat0) =',i2,f7.1)") | nalt,t_curmo(lon0,lat0) u_curmo = ctmtdat(mon_cur,ihr_cur+1,'u',nalt) v_curmo = ctmtdat(mon_cur,ihr_cur+1,'v',nalt) ! call addfld('d_curmo','d_curmo',' ',d_curmo, ! | 'lon',i0,i1,'lat',j0,j1,0) ! call addfld('t_curmo','t_curmo',' ',t_curmo, ! | 'lon',i0,i1,'lat',j0,j1,0) ! call addfld('u_curmo','u_curmo',' ',u_curmo, ! | 'lon',i0,i1,'lat',j0,j1,0) ! call addfld('v_curmo','v_curmo',' ',v_curmo, ! | 'lon',i0,i1,'lat',j0,j1,0) ! ! Subdomains at next month, current hour: d_nxtmo = ctmtdat(mon_nxt,ihr_cur+1,'d',nalt) t_nxtmo = ctmtdat(mon_nxt,ihr_cur+1,'t',nalt) ! TEMP write (6,"(1x,'t_nxtmo(lon0,lat0) =',f7.1)") | nalt,t_nxtmo(lon0,lat0) u_nxtmo = ctmtdat(mon_nxt,ihr_cur+1,'u',nalt) v_nxtmo = ctmtdat(mon_nxt,ihr_cur+1,'v',nalt) ! call addfld('d_nxtmo','d_nxtmo',' ',d_nxtmo, ! | 'lon',i0,i1,'lat',j0,j1,0) ! call addfld('t_nxtmo','z_nxtmo',' ',t_nxtmo, ! | 'lon',i0,i1,'lat',j0,j1,0) ! call addfld('u_nxtmo','z_nxtmo',' ',u_nxtmo, ! | 'lon',i0,i1,'lat',j0,j1,0) ! call addfld('v_nxtmo','z_nxtmo',' ',v_nxtmo, ! | 'lon',i0,i1,'lat',j0,j1,0) ! ! Interpolate to month: if (isecs /= 0) then difsec = 60.*60. ! difference in sec between ihr_cur and ihr_nxt secint = float(modeltime(3)*60 + modeltime(4)) ! interpolation time [sec] ! ! Current month, next hour: d = d_curmo ; t = t_curmo ; u = u_curmo ; v = v_curmo d_curmo_nxthr = ctmtdat(mon_cur,ihr_nxt+1,'d',nalt) t_curmo_nxthr = ctmtdat(mon_cur,ihr_nxt+1,'t',nalt) ! TEMP write (6,"(1x,'nalt t_curmo_nxthr(lon0,lat0) =',i2,f7.1)") | nalt,t_curmo_nxthr(lon0,lat0) u_curmo_nxthr = ctmtdat(mon_cur,ihr_nxt+1,'u',nalt) v_curmo_nxthr = ctmtdat(mon_cur,ihr_nxt+1,'v',nalt) call timeinterp(d, d_curmo_nxthr, difsec, secint, d_curmo, 1) call timeinterp(t, t_curmo_nxthr, difsec, secint, t_curmo, 1) call timeinterp(u, u_curmo_nxthr, difsec, secint, u_curmo, 1) call timeinterp(v, v_curmo_nxthr, difsec, secint, v_curmo, 1) ! ! Interpolate to next month: d = d_nxtmo ; t = t_nxtmo ; u = u_nxtmo ; v = v_nxtmo ! ! Next month, next hour: d_nxtmo_nxthr = ctmtdat(mon_nxt,ihr_nxt+1,'d',nalt) t_nxtmo_nxthr = ctmtdat(mon_nxt,ihr_nxt+1,'t',nalt) ! TEMP write (6,"(1x,'t_nxtmo_nxthr(lon0,lat0) =',f7.1)") | nalt,t_nxtmo_nxthr(lon0,lat0) u_nxtmo_nxthr = ctmtdat(mon_nxt,ihr_nxt+1,'u',nalt) v_nxtmo_nxthr = ctmtdat(mon_nxt,ihr_nxt+1,'v',nalt) call timeinterp(d, d_nxtmo_nxthr, difsec, secint, d_nxtmo,1) call timeinterp(t, t_nxtmo_nxthr, difsec, secint, t_nxtmo,1) call timeinterp(u, u_nxtmo_nxthr, difsec, secint, u_nxtmo,1) call timeinterp(v, v_nxtmo_nxthr, difsec, secint, v_nxtmo,1) endif ! interpolate to month or not ! ! Check if interpolation to ut is necessary nointp = 0 if(iday.eq.ndmon_nl(mon_cur)) then ! same day as cur. month nointp = 1 goto 20 endif if(iday.eq.ndmon_nl(mon_nxt)) then ! same day as next month nointp = 2 goto 20 endif ! ! If ut interpolation is necessary, calculate time differences if(mon_cur /= 12) then ! not December difday = ndmon_nl(mon_nxt)-ndmon_nl(mon_cur) ! difference in days dayint = iday - ndmon_nl(mon_cur) ! difference to interpolation day else ! December wrap around difday = ndmon_nl(mon_cur+1)-ndmon_nl(mon_cur) ! difference in days if(iday.lt.ndmon_nl(mon_nxt)) then ! difference to interpolation day dayint = 365. - ndmon_nl(mon_cur)+ iday else dayint = iday - ndmon_nl(mon_cur) endif endif 20 continue ! if no interpolation is necessary (nointp /= 0) ! ! Interpolate to ut if necessary to local d,t,u,v: ! TEMP for 1 MPI processor, do whole at once to find out why 1,17=0 always!? ! write (6,"(1x,'nalt t_curmo,nxtmo_nxthr(1,17) =', ! | i2,4f7.1)") nalt,t_curmo(1,17),t_curmo_nxthr(1,17), ! | t_nxtmo(1,17),t_nxtmo_nxthr(1,17) select case (nointp) case (0) ! interpolate call timeinterp(d_curmo,d_nxtmo,difday,dayint,d,1) call timeinterp(t_curmo,t_nxtmo,difday,dayint,t,1) call timeinterp(u_curmo,u_nxtmo,difday,dayint,u,1) call timeinterp(v_curmo,v_nxtmo,difday,dayint,v,1) case (1) ! no interp (same day as current month) d = d_curmo t = t_curmo u = u_curmo v = v_curmo case (2) ! no interp (same day as next month) d = d_nxtmo t = t_nxtmo u = u_nxtmo v = v_nxtmo case default write(6,"(/,'>>> mkctmt: unknown nointp=',i4)") nointp call shutdown('mkctmt') end select ! ! Transfer to 3D module data ctmt_d(lon0:lon1,lat0:lat1,nalt) = d ctmt_t(lon0:lon1,lat0:lat1,nalt) = t ctmt_u(lon0:lon1,lat0:lat1,nalt) = u ctmt_v(lon0:lon1,lat0:lat1,nalt) = v ! ! Do mpi periodic points exchange for gswm with f2d(:) ! lons 1,2 <- nlonp4-3,nlonp4-2 and nlonp4-1,nlonp4 <- 3,4 ! #ifdef MPI fctmt(:,:,1) = ctmt_d(lon0:lon1,lat0:lat1,nalt) fctmt(:,:,2) = ctmt_t(lon0:lon1,lat0:lat1,nalt) fctmt(:,:,3) = ctmt_u(lon0:lon1,lat0:lat1,nalt) fctmt(:,:,4) = ctmt_v(lon0:lon1,lat0:lat1,nalt) call mp_periodic_f2d(fctmt,lon0,lon1,lat0,lat1,4) ctmt_d(lon0:lon1,lat0:lat1,nalt) = fctmt(:,:,1) ctmt_t(lon0:lon1,lat0:lat1,nalt) = fctmt(:,:,2) ctmt_u(lon0:lon1,lat0:lat1,nalt) = fctmt(:,:,3) ctmt_v(lon0:lon1,lat0:lat1,nalt) = fctmt(:,:,4) #else ! lons 1,2 <- nlonp4-3,nlonp4-2 and nlonp4-1,nlonp4 <- 3,4 ctmt_d(1:2,:,nalt) = ctmt_d(nlonp4-3:nlonp4-2,:,nalt) ctmt_d(nlonp4-1:nlonp4,:,nalt) = ctmt_d(3:4,:,nalt) ctmt_t(1:2,:,nalt) = ctmt_t(nlonp4-3:nlonp4-2,:,nalt) ctmt_t(nlonp4-1:nlonp4,:,nalt) = ctmt_t(3:4,:,nalt) ctmt_u(1:2,:,nalt) = ctmt_u(nlonp4-3:nlonp4-2,:,nalt) ctmt_u(nlonp4-1:nlonp4,:,nalt) = ctmt_u(3:4,:,nalt) ctmt_v(1:2,:,nalt) = ctmt_v(nlonp4-3:nlonp4-2,:,nalt) ctmt_v(nlonp4-1:nlonp4,:,nalt) = ctmt_v(3:4,:,nalt) ! call set_periodic_f2d(ctmt_d2d) ! call set_periodic_f2d(ctmt_t2d) ! call set_periodic_f2d(ctmt_u2d) ! call set_periodic_f2d(ctmt_v2d) #endif ! TEMP write-out (but not zero here??) ! write (6,"(1x,'mkctmt m,nalt ctmt_t(1,17,nalt)=',2i2,3f7.1)") ! | malt,nalt,ctmt_t(1,17,nalt) enddo ! malt=1,nalt5 ! TEMP write-out and save altitudes do nalt=1,nalt5 ! write (6,"(1x,'mkctmt nalt ctmt_t(1,17,nalt)=',i2,3f7.1)") ! | nalt,ctmt_t(1,17,nalt) call addfld('TIDE_D','TIDE D','percent', | ctmt_d(i0:i1,j0:j1,nalt),'lon',i0,i1,'lat',j0,j1,nalt) call addfld('TIDE_T','TIDE TN','K', | ctmt_t(i0:i1,j0:j1,nalt),'lon',i0,i1,'lat',j0,j1,nalt) call addfld('TIDE_U','TIDE UN','cm/s', | ctmt_u(i0:i1,j0:j1,nalt),'lon',i0,i1,'lat',j0,j1,nalt) call addfld('TIDE_V','TIDE VN','cm/s', | ctmt_v(i0:i1,j0:j1,nalt),'lon',i0,i1,'lat',j0,j1,nalt) enddo end subroutine mkctmt !----------------------------------------------------------------------- function ctmtdat(mon,ihr,fname,nalt) ! ! Return subdomain array from read arrays according to requested field, ! month, and hour. ! The 4 global fields t_all, etc, are private module data, read by ! sub rdctmt at beginning of model run. ! ! Args: character(len=*),intent(in) :: fname integer,intent(in) :: mon,ihr,nalt ! ! Function output dimension: real :: ctmtdat(lon0:lon1,lat0:lat1) ! Local: real :: ut,amplat,phlat,facunits(4) integer :: nct,nn,nc,nv,ss,jl,il data facunits/1.,1.,100.,100./ ! Convert u,v from m/s to cm/s ! ! Field fname must be t,u,v, or d: if (fname/='t'.and.fname/='u'.and.fname/='v'.and.fname/='d') then write(6,"(/,'>>> ctmtdat: unknown fname=',a)") fname write(6,"('Must be t, u, v, or d')") call shutdown('ctmtdat') endif if (fname=='d') nv = 1 if (fname=='t') nv = 2 if (fname=='u') nv = 3 if (fname=='v') nv = 4 ! Calculate the total migrating tide and non-migrating tide at UT=ihr-1 ! migrating tide is W1 for diurnal, W2 for sd, where ss=nn (don't separate out) ! Set tides to zero since 6 or 8 tides must be added together ut = ihr - 1 do jl=lat0,lat1 do il=lon0,lon1 ctmtdat(il,jl) = 0. enddo enddo do jl=lat0,lat1 do nn=1,2 nct = 6 if (nn==2) nct = 8 do nc=1,nct ss = nc - 4 ! Find amp and phase interpolated in lat to use multiplied by fac and by ! and by facunits (m/s to cm/s) and shifted by phsft amplat = (faclat1(jl)*amplb(nn,nc,nv,mon,jl1(jl),nalt) | + faclat2(jl)*amplb(nn,nc,nv,mon,jl1(jl)+1,nalt) ) * | fac(nn,nc) * facunits(nv) ! NOT DONE YET - CHECK TO SEE IF PHASE CROSSES 0/24UT TIME!! phlat = faclat1(jl)*phlb(nn,nc,nv,mon,jl1(jl),nalt) + | faclat2(jl)*phlb(nn,nc,nv,mon,jl1(jl)+1,nalt) - phsft(nn,nc) ! tide=amp*cos[n*omega*LT +(s-n)glon -phase] or amp*cos[n*omega*UT +s*glon -phase do il=lon0,lon1 ctmtdat(il,jl) = ctmtdat(il,jl) + amplat * | cos(nn*pi*(ut-phlat)/12. + ss*glonp4(il)*pi/180.) ! Assume wrap-around ctmtdat enddo ! il=lon0:lon1 enddo ! nc=1,nct ! TEMP - write out for equ nlat/2 or 36 for dres if (nn==1 .and. nv==3 .and. jl==36) write (6,"(1x, | 'nn=1(D) Un amplat phlat=',2f9.3)") amplat,phlat if (nn==2 .and. nv==3 .and. jl==36) write (6,"(1x, | 'nn=2(SD) Un amplat phlat=',2f9.3)") amplat,phlat enddo ! nn=1,2 enddo ! jl=1,nlat end function ctmtdat !----------------------------------------------------------------------- subroutine timeinterp(d1,d2,difd1d2,difint,fout, | iprnt) ! ! Interpolate from fields d1,d2 linearly to time difint ! (d1 must be at 0 unit time), returning fout ! On input: ! d1,d2 = the field (d1 at 0 unit, d2 at difd1d2 units) ! difd1d2 = time difference [unit] between d1 & d2 ! difint = time of interpolation [unit] (counted from d1=0 time) ! In output: ! fout is defined at difint ! ! Args: integer,intent(in) :: iprnt real,intent(in) :: difd1d2,difint,d1(lon0:lon1,lat0:lat1), | d2(lon0:lon1,lat0:lat1) real,intent(out) :: fout(lon0:lon1,lat0:lat1) ! ! Local: integer :: i,j real :: frac0,difd1d2_inv ! fout = 0. ! initialize output ! ! Interpolation: ! from data_cur/d1 to data_nxt/d2 the difference is in time difd1d2 ! difd1d2_inv = 1./difd1d2 frac0 = difint*difd1d2_inv ! x(int)/[x(d2)-x(d1)] ! linear interpolation: special case x(d1) = 0 ! fout = (d2-d1)*frac0 + d1 - (d2-d1)*difd1d2_inv*x(d1) ! (d2-d1)*difd1d2_inv*x(d1) = 0 since x(d1) = 0 do i = lon0,lon1 do j = lat0,lat1 fout(i,j) = (d2(i,j)-d1(i,j))*frac0 + d1(i,j) enddo enddo end subroutine timeinterp !----------------------------------------------------------------------- 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,:) end subroutine set_periodic_f2d !----------------------------------------------------------------------- subroutine alloc_ctmt(lon0,lon1,lat0,lat1) ! Args: integer,intent(in) :: lon0,lon1,lat0,lat1 ! ! Local: integer :: istat ! ! CTMT data in altitude to be use-associated to lbc.F: allocate(ctmt_d (lon0:lon1,lat0:lat1,nalt5),stat=istat) allocate(ctmt_t (lon0:lon1,lat0:lat1,nalt5),stat=istat) allocate(ctmt_u (lon0:lon1,lat0:lat1,nalt5),stat=istat) allocate(ctmt_v (lon0:lon1,lat0:lat1,nalt5),stat=istat) write(6,"('Allocated ctmt t,u,v,z (lon0:lon1,lat0:lat1,nalt5)')") end subroutine alloc_ctmt !----------------------------------------------------------------------- end module ctmt_module !----------------------------------------------------------------------- subroutine cal_ctmt_z(istep,itp,tbgrd,zprt,tprt,rhoprt | ,izm,zmtbgrd,zmzprt,zmtprt,zmrhoprt,k1) ! cal_ctmt_z based on Maute's cal_hme_z and called from lbc.F ! with 2 extra arguements istep,itp in tuzvz_lbc called from advance.F ! ilev = interface points for Z,ZG,DEN,NE,UI_ExB,VI_ExB,WI_ExB ! lev(=ilev) midpoints for UN,VN,WN,TN,TI,TE,O1,O2,NO use mpi_module,only: lon0,lon1,lat0,lat1 use init_module,only: glat,zpint ! added for output to fort.41 for nproc=1 use params_module,only: nlonp4,nlat,nlon,dlat use cons_module,only: | grav, ! 870. [cm/s2] accel due to gravity | boltz, ! 1.38E-16 [ergs/K] boltzman's constant | rmassinv, | avo, ! 6.023e23 [#/mol] Avogadro constant | pi use fields_module,only: o2,o1,barm implicit none ! ! calculates the zpert ! zprt = -R/g0*tprt - R/g0*tbgrd*(pprt/p0 + tprt/tbgrd) ! with R = kb/m [cm2/s2] kb Boltzman constant; m mass ! integer,intent(in):: istep, ! step | itp ! index of previous step | ,izm ! flag=1 if doing zonal mean | ,k1 ! value of bottom interface zpint(k1) level where 1=extrap real,intent(in) :: | tbgrd(nlonp4,nlat), ! t_lbc | tprt(lon0:lon1,lat0:lat1), ! tn pert CTMT [K] | rhoprt(lon0:lon1,lat0:lat1) ! rho_pertr/rho_bgrd CTMT [-] real,intent(in) :: | zmtbgrd, ! t_lbc | zmtprt(72), ! tn pert zonal mean [K] | zmrhoprt(72) ! rho_pertr/rho_bgrd zonal mean [%] real,intent(out):: zprt(lon0:lon1,lat0:lat1) real,intent(out):: zmzprt(72) integer :: i,j,k,k0 real :: barm_lb(lon0:lon1,lat0:lat1),barm_zmlbc(lat0:lat1), | barm_tmp(2),barm_lbc(lon0:lon1,lat0:lat1), | barm_zmlb(lat0:lat1) real :: dtr, re, areatot, glob_barm, dglat, p_lev, da(72) ! if(istep == 1) then ! ! barm = mean molecular weight (k+1/2): [g/mol] ! barm_lb [g] ! do j=lat0,lat1 do i=lon0,lon1 ! 10/14: k=1,2 means barm ONLY correct for lbc! k0 = 0 do k=k1,k1+1 k0 = k0 + 1 barm_tmp(k0) = 1./ | (o2(k,i,j,itp)*rmassinv(1)+o1(k,i,j,itp)*rmassinv(2)+ | (1.-o2(k,i,j,itp)-o1(k,i,j,itp))*rmassinv(3)) enddo if (k1==1) then ! Extrapolate below from mid-points barm_lbc(i,j) = 1.5*barm_tmp(1)-0.5*barm_tmp(2) else ! Interpolate half-way barm_lbc(i,j) = 0.5*(barm_tmp(1)+barm_tmp(2)) endif barm_lb(i,j) = barm_lbc(i,j)/avo enddo enddo ! else ! not first timestep ! ! barm1 = barm(k=0) (linear extrapolation) ! ! do j=lat0,lat1 ! do i=lon0,lon1 ! barm_lb(i,j) = 1.5*barm(1,i,j,itp)-0.5*barm(2,i,j,itp) ! barm_lb(i,j) = barm_lb(i,j)/avo ! enddo ! enddo ! endif ! end first timestep ! if (izm==1) then ! 10/14: Calc global mean barm at several P levels dtr = pi/180. re = 6371.e5 areatot = 0. glob_barm = 0. do j=lat0,lat1 barm_zmlbc(j) = 0. barm_zmlb(j) = 0. da(j) = re*cos(glat(j)*dtr)*(dlat*dtr)*(re*dlat*dtr) areatot = areatot + da(j)*nlon do i=1,nlon barm_zmlbc(j) = barm_zmlbc(j) + barm_lbc(i,j) barm_zmlb(j) = barm_zmlb(j) + barm_lb(i,j) enddo ! NOTE: ZM ONLY correct if nproc=1!!! barm_zmlbc(j) = barm_zmlbc(j)/nlon barm_zmlb(j) = barm_zmlb(j)/nlon glob_barm = glob_barm + barm_zmlbc(j)*da(j)*nlon enddo glob_barm = glob_barm/areatot p_lev = 5e-7*exp(-zpint(k1)) write (41,"(i4,e12.4,3f9.3,' k1,p_lev,glob_barm,zpint(k1),', | 'zmtbgrd')") k1,p_lev,glob_barm,zpint(k1),zmtbgrd endif ! kb/m/g*Tn [ergs/K /g / (cm/s2) *K] = [cm] do j=lat0,lat1 do i=lon0,lon1 if (izm==0) then ! 12/14 bae: This is wrong, with -2Tpert - revise! ! zprt(i,j) = - boltz/grav/barm_lb(i,j)*tprt(i,j)- ! | boltz/grav/barm_lb(i,j)*tbgrd(i,j)* ! | (rhoprt(i,j)+tprt(i,j)/tbgrd(i,j)) zprt(i,j) = - boltz/grav/barm_lb(i,j)*tprt(i,j)- | boltz/grav/barm_lb(i,j)*tbgrd(i,j)*rhoprt(i,j) endif ! write(6,"('cal_ctmt_z: istep=',i3,' j=',i3,' i=',i3, ! | ' barm_lb=',e12.4,' tprt=',e12.4,' tbgrd=',e12.4, ! | ' rhoprt=',e12.4,' zprt=',e12.4)") istep,j,i, ! | barm_lb(i,j),tprt(i,j),tbgrd(i,j),rhoprt(i,j),zprt(i,j) ! write(6,*) 'CTMT1',i,j,boltz/grav/barm_lb(i,j)*tprt(i,j) ! write(6,*) 'CTMT2',i,j, ! | boltz/grav/barm_lb(i,j)*tbgrd(i,j)* ! | (rhoprt(i,j)) ! write(6,*) 'CTMT3',i,j, ! | boltz/grav/barm_lb(i,j)*tbgrd(i,j)* ! | (tprt(i,j)/tbgrd(i,j)) enddo if (izm==1) then ! 10/14 bae: This is wrong, with -2Tpert/Tbckg - always neg - revise! ! zmzprt(j) = - boltz/grav/barm_zmlb(j)*zmtprt(j)- ! | boltz/grav/barm_zmlb(j)*zmtbgrd* ! | (zmrhoprt(j)+zmtprt(j)/zmtbgrd) zmzprt(j) = - boltz/grav/barm_zmlb(j)*zmtprt(j)- | boltz/grav/barm_zmlb(j)*zmtbgrd*zmrhoprt(j) write (41,"(1x,i3,f8.3,2e14.5,3f8.3,3e14.5)") | j,barm_zmlbc(j),barm_zmlb(j),boltz/grav/barm_zmlb(j), | zmtprt(j),zmrhoprt(j),zmtbgrd*zmrhoprt(j),zmzprt(j), | -boltz/grav/barm_zmlb(j)*zmtprt(j), | -boltz/grav/barm_zmlb(j)*zmtbgrd*zmrhoprt(j) endif enddo ! write(6,"('cal_cmtm_z: istep=',i3,' barm_lb min,max=',2e12.4, ! | ' zprt min,max=',2e12.4)") istep,minval(barm_lb), ! | maxval(barm_lb),minval(zprt),maxval(zprt) end subroutine cal_ctmt_z