! 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 ! ! Variable naming convention for CTMT tidal components: ! _mi_di_ Migrating diurnal ! _mi_sdi_ Migrating semi-diurnal ! _nm_di_ Non-migrating diurnal ! _nm_sdi_ Non-migrating semi-diurnal ! use params_module,only: nlon,nlat,nlonp4,nlonp2 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 implicit none ! #ifdef SUN #include #else #include #endif ! ! CTMT boundary data at model grid, output by this module, ! to be used by dt, duv, addiag, etc: ! ! allocated subdomains: ! real,dimension(:,:) :: ! will be (lon0:lon1,lat0,lat1) ! real,dimension(nlonp4,nlat) :: | ctmt_mi_di_z, ctmt_mi_sdi_z, ctmt_nm_di_z, ctmt_nm_sdi_z, | ctmt_mi_di_t, ctmt_mi_sdi_t, ctmt_nm_di_t, ctmt_nm_sdi_t, | ctmt_mi_di_u, ctmt_mi_sdi_u, ctmt_nm_di_u, ctmt_nm_sdi_u, | ctmt_mi_di_v, ctmt_mi_sdi_v, ctmt_nm_di_v, ctmt_nm_sdi_v ! ! 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 ! Convert to 97km (use 97.5km at #40=nlb) at UT, for nlat glats and nlonp4 glons integer,parameter,private :: nmonth= 12, nhour=24, nlat5=37, | nalt2p5=161, nlb=40 real,dimension(2,nlonp4,nlat,nmonth,nhour),private :: | z_mi, z_nm, t_mi, t_nm, u_mi, u_nm, v_mi, v_nm, d_mi, d_mn ! contains !----------------------------------------------------------------------- subroutine getctmt(istep,iyear,iday,secs) ! ! Module driver to read nc files, and do time interpolations. ! ! 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 ! ! Integer flags set according to user-requested files: use init_module,only: | ictmt ! 0/1 flag to get CTMT data for 6 diurnal and 8 semidiurnal tides ! ! Driver for obtaining CTMT data, called once per timestep from advance. ! ! Args: integer,intent(in) :: istep,iyear,iday real,intent(in) :: secs ! ! Local: integer :: iprint ! iprint = 0 if (istep==1) iprint = 1 ! ! Get data for 6 CTMT diurnal tides and 8 CTMT semi-diurnal tides: if (ictmt > 0) then if (istep==1) then call rdctmt(ctmt_di_ncfile,'diur') call rdctmt(ctmt_sd_ncfile,'semi') end call mkctmt(iyear,iday,int(secs),iprint) endif 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_nmonth, id_nlat, id_am, id_ph character(len=240) :: char240 real,dimension(nlat) :: faclat1,faclat2 integer,dimension(nlat) :: jl1 real,dimension(nlat5) :: ctmt_lat real,dimension(nmonth,nalt2p5,nlat5) :: amptmp,phtmp real,dimension(nlon) :: aa,bb,ampnm,phnm integer :: ishift,i,im,il,jl character(len=16) :: cmp character(len=4) :: var character(len=8) :: ampnm character(len=10) :: phnm real :: fac(2,8) ! amplitude mulipliers for 1,2 di,sd and 8 cmp components real :: phsft(2,8) ! phase shifts in hours for 1,2 di,sd and 8 cmp real :: ut,pi,amplat(4),phlat(4),facunits(4) integer :: mo,nut,nct,nn,nc,nv,ss data cmp/'e3e2e1s0w1w2w3w4'/ ! 8 tides for sd (first 6 for diurnal) data var/'dtuv'/ ! Use 4 of 5 CTMT variables (skip 'w') data facunits/1.,1.,100.,100./ ! Convert u,v from m/s to cm/s pi = 4.*atan(1.) ! Set amplitude factors for CTMT diurnal and semi-diurnal tides (default=1) ! and set phase shifts (default=0 h) do nn=1,2 do nc=1,8 fac(nn,nc) = 1. phsft(nn,nc) = 0. enddo enddo 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 ! ! 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(ncfiledi) 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(ncfiledi) call shutdown('rdctmt') endif ! Get CTMT glat ctmt_lat - should be -90:5:90 istat = nf_get_var(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)) 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 ! ! 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(ncid,id_am,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,:) 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(ncid,id_am,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 enddo enddo ! Calculate the migrating tide nmig and non-migrating tide at 24 UTs ! migrating tide is W1 for diurnal, W2 for sd, where ss=nn ! Set non-migrating tides to zero since 5 or 7 tides must be added together do mo=1,12 do jl=1,nlat do il=3,nlon+2 do nut=1,24 d_nm(nn,il,jl,mo,nut) = 0. t_nm(nn,il,jl,mo,nut) = 0. u_nm(nn,il,jl,mo,nut) = 0. v_nm(nn,il,jl,mo,nut) = 0. enddo enddo enddo enddo do mo=1,12 do jl=1,nlat do nc=1,nct ss = nc - 4 ! Find amp and phase interpolated in lat to use multiplied by fac and shifted by phsft do nv=1,4 amplat(nv) = faclat1(jl)*amplb(nn,nc,nv,mo,jl1(jl)) + | faclat2(jl)*amplb(nn,nmig,nv,mo,jl1(jl)+1) ) * fac(nn,nmig) | * facunits(nv) phlat(nv) = faclat1(jl)*phlb(nn,nc,nv,mo,jl1(jl)) + | faclat2(jl)*phlb(nn,nmig,nv,mo,jl1(jl)+1) - phsft(nn,nmig) enddo ! nv=1,4 ! tide=amp*cos[n*omega*LT +(s-n)glon -phase] or amp*cos[n*omega*UT +s*glon -phase] do nut=1,24 ut = nut - 1 if (nn==ss) then ! Set migrating tides do il=3,nlon+2 d_mi(nn,il,jl,mo,nut) = amplat(1) * | cos(nn*pi*(ut-phlat(1))/24.) + ss*glon(il-2)*pi/180.) t_mi(nn,il,jl,mo,nut) = amplat(2) * | cos(nn*pi*(ut-phlat(2))/24.) + ss*glon(il-2)*pi/180.) u_mi(nn,il,jl,mo,nut) = amplat(3) * | cos(nn*pi*(ut-phlat(3))/24.) + ss*glon(il-2)*pi/180.) v_mi(nn,il,jl,mo,nut) = amplat(4) * | cos(nn*pi*(ut-phlat(4))/24.) + ss*glon(il-2)*pi/180.) enddo ! il=3,nlon+2 else ! Set non-migrating tides do il=3,nlon+2 d_nm(nn,il,jl,mo,nut) = d_nm(nn,il,jl,mo,nut) + amplat(1) | * cos(nn*pi*(ut-phlat(1))/24.) + ss*glon(il-2)*pi/180.) t_nm(nn,il,jl,mo,nut) = t_nm(nn,il,jl,mo,nut) + amplat(2) | * cos(nn*pi*(ut-phlat(2))/24.) + ss*glon(il-2)*pi/180.) u_nm(nn,il,jl,mo,nut) = u_nm(nn,il,jl,mo,nut) + amplat(3) | * cos(nn*pi*(ut-phlat(3))/24.) + ss*glon(il-2)*pi/180.) v_nm(nn,il,jl,mo,nut) = v_nm(nn,il,jl,mo,nut) + amplat(4) | * cos(nn*pi*(ut-phlat(4))/24.) + ss*glon(il-2)*pi/180.) enddo ! il=3,nlon=2 endif ! if,then,else enddo ! nut=1,24 enddo ! nc=1,nct ! Add wrap-around pts do nut=1,24 d_mi(nn,1,jl,mo,nut) = d_mi(nn,nlon+1,jl,mo,nut) d_mi(nn,2,jl,mo,nut) = d_mi(nn,nlon+2,jl,mo,nut) d_mi(nn,nlon+1,jl,mo,nut) = d_mi(nn,3,jl,mo,nut) d_mi(nn,nlon+2,jl,mo,nut) = d_mi(nn,4,jl,mo,nut) t_mi(nn,1,jl,mo,nut) = t_mi(nn,nlon+1,jl,mo,nut) t_mi(nn,2,jl,mo,nut) = t_mi(nn,nlon+2,jl,mo,nut) t_mi(nn,nlon+1,jl,mo,nut) = t_mi(nn,3,jl,mo,nut) t_mi(nn,nlon+2,jl,mo,nut) = t_mi(nn,4,jl,mo,nut) u_mi(nn,1,jl,mo,nut) = u_mi(nn,nlon+1,jl,mo,nut) u_mi(nn,2,jl,mo,nut) = u_mi(nn,nlon+2,jl,mo,nut) u_mi(nn,nlon+1,jl,mo,nut) = u_mi(nn,3,jl,mo,nut) u_mi(nn,nlon+2,jl,mo,nut) = u_mi(nn,4,jl,mo,nut) v_mi(nn,1,jl,mo,nut) = v_mi(nn,nlon+1,jl,mo,nut) v_mi(nn,2,jl,mo,nut) = v_mi(nn,nlon+2,jl,mo,nut) v_mi(nn,nlon+1,jl,mo,nut) = v_mi(nn,3,jl,mo,nut) v_mi(nn,nlon+2,jl,mo,nut) = v_mi(nn,4,jl,mo,nut) d_nm(nn,1,jl,mo,nut) = d_nm(nn,nlon+1,jl,mo,nut) d_nm(nn,2,jl,mo,nut) = d_nm(nn,nlon+2,jl,mo,nut) d_nm(nn,nlon+1,jl,mo,nut) = d_nm(nn,3,jl,mo,nut) d_nm(nn,nlon+2,jl,mo,nut) = d_nm(nn,4,jl,mo,nut) t_nm(nn,1,jl,mo,nut) = t_nm(nn,nlon+1,jl,mo,nut) t_nm(nn,2,jl,mo,nut) = t_nm(nn,nlon+2,jl,mo,nut) t_nm(nn,nlon+1,jl,mo,nut) = t_nm(nn,3,jl,mo,nut) t_nm(nn,nlon+2,jl,mo,nut) = t_nm(nn,4,jl,mo,nut) u_nm(nn,1,jl,mo,nut) = u_nm(nn,nlon+1,jl,mo,nut) u_nm(nn,2,jl,mo,nut) = u_nm(nn,nlon+2,jl,mo,nut) u_nm(nn,nlon+1,jl,mo,nut) = u_nm(nn,3,jl,mo,nut) u_nm(nn,nlon+2,jl,mo,nut) = u_nm(nn,4,jl,mo,nut) v_nm(nn,1,jl,mo,nut) = v_nm(nn,nlon+1,jl,mo,nut) v_nm(nn,2,jl,mo,nut) = v_nm(nn,nlon+2,jl,mo,nut) v_nm(nn,nlon+1,jl,mo,nut) = v_nm(nn,3,jl,mo,nut) v_nm(nn,nlon+2,jl,mo,nut) = v_nm(nn,4,jl,mo,nut) enddo ! nut=1,24 enddo ! jl=1,nlat enddo ! mo=1,12 ! Get Z geopotential height perturbation [m] and convert to cm (and U,V in m/s to cm/s): ! ! z = z*100. ! convert from m to cm ! ! 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 mkgswm(iyear,iday,isecs,iprint,type) ! ! Use data read from gswm file to return perturbation in Z ! and T at current model date and time. It is assumed that the gswm 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 character(len=*),intent(in) :: type ! ! 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 real,dimension(lon0:lon1,lat0:lat1) :: t,u,v,z, | z_curmo, t_curmo, u_curmo, v_curmo, | z_nxtmo, t_nxtmo, u_nxtmo, v_nxtmo, | z_curmo_nxthr, t_curmo_nxthr, u_curmo_nxthr, v_curmo_nxthr, | z_nxtmo_nxthr, t_nxtmo_nxthr, u_nxtmo_nxthr, v_nxtmo_nxthr real :: fgswm(lon0:lon1,lat0:lat1,4) ! for mp calls ! ! Data: GSWM 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 if(ihr_nxt == 24 ) ihr_nxt = 1 ! ! Subdomains at current month, current hour: z_curmo = gswmdat(type,mon_cur,ihr_cur+1,'z') t_curmo = gswmdat(type,mon_cur,ihr_cur+1,'t') u_curmo = gswmdat(type,mon_cur,ihr_cur+1,'u') v_curmo = gswmdat(type,mon_cur,ihr_cur+1,'v') ! call addfld('z_curmo','z_curmo',' ',z_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: z_nxtmo = gswmdat(type,mon_nxt,ihr_cur+1,'z') t_nxtmo = gswmdat(type,mon_nxt,ihr_cur+1,'t') u_nxtmo = gswmdat(type,mon_nxt,ihr_cur+1,'u') v_nxtmo = gswmdat(type,mon_nxt,ihr_cur+1,'v') ! call addfld('z_nxtmo','z_nxtmo',' ',z_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: z = z_curmo ; t = t_curmo ; u = u_curmo ; v = v_curmo z_curmo_nxthr = gswmdat(type,mon_cur,ihr_nxt+1,'z') t_curmo_nxthr = gswmdat(type,mon_cur,ihr_nxt+1,'t') u_curmo_nxthr = gswmdat(type,mon_cur,ihr_nxt+1,'u') v_curmo_nxthr = gswmdat(type,mon_cur,ihr_nxt+1,'v') call timeinterp(z, z_curmo_nxthr, difsec, secint, z_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: z = z_nxtmo ; t = t_nxtmo ; u = u_nxtmo ; v = v_nxtmo ! ! Next month, next hour: z_nxtmo_nxthr = gswmdat(type,mon_nxt,ihr_nxt+1,'z') t_nxtmo_nxthr = gswmdat(type,mon_nxt,ihr_nxt+1,'t') u_nxtmo_nxthr = gswmdat(type,mon_nxt,ihr_nxt+1,'u') v_nxtmo_nxthr = gswmdat(type,mon_nxt,ihr_nxt+1,'v') call timeinterp(z, z_nxtmo_nxthr, difsec, secint, z_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 z,t,u,v: select case (nointp) case (0) ! interpolate call timeinterp(z_curmo,z_nxtmo,difday,dayint,z,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) z = z_curmo t = t_curmo u = u_curmo v = v_curmo case (2) ! no interp (same day as next month) z = z_nxtmo t = t_nxtmo u = u_nxtmo v = v_nxtmo case default write(6,"(/,'>>> mkgswm: unknown nointp=',i4)") nointp call shutdown('mkgswm') end select ! ! Transfer to module data according to type: select case (type) case ('mi_di') gswm_mi_di_z(lon0:lon1,lat0:lat1) = z gswm_mi_di_t(lon0:lon1,lat0:lat1) = t gswm_mi_di_u(lon0:lon1,lat0:lat1) = u gswm_mi_di_v(lon0:lon1,lat0:lat1) = v case ('mi_sdi') gswm_mi_sdi_z(lon0:lon1,lat0:lat1) = z gswm_mi_sdi_t(lon0:lon1,lat0:lat1) = t gswm_mi_sdi_u(lon0:lon1,lat0:lat1) = u gswm_mi_sdi_v(lon0:lon1,lat0:lat1) = v case ('nm_di') gswm_nm_di_z(lon0:lon1,lat0:lat1) = z gswm_nm_di_t(lon0:lon1,lat0:lat1) = t gswm_nm_di_u(lon0:lon1,lat0:lat1) = u gswm_nm_di_v(lon0:lon1,lat0:lat1) = v case ('nm_sdi') gswm_nm_sdi_z(lon0:lon1,lat0:lat1) = z gswm_nm_sdi_t(lon0:lon1,lat0:lat1) = t gswm_nm_sdi_u(lon0:lon1,lat0:lat1) = u gswm_nm_sdi_v(lon0:lon1,lat0:lat1) = v case default write(6,"(/,'>>> mkgswm: unknown type=',a)") type call shutdown('mkgswm') end select ! ! Do mpi periodic points exchange for gswm with f2d(:) ! lons 1,2 <- nlonp4-3,nlonp4-2 and nlonp4-1,nlonp4 <- 3,4 ! #ifdef MPI select case (type) case ('mi_di') fgswm(:,:,1) = gswm_mi_di_z(lon0:lon1,lat0:lat1) fgswm(:,:,2) = gswm_mi_di_t(lon0:lon1,lat0:lat1) fgswm(:,:,3) = gswm_mi_di_u(lon0:lon1,lat0:lat1) fgswm(:,:,4) = gswm_mi_di_v(lon0:lon1,lat0:lat1) call mp_periodic_f2d(fgswm,lon0,lon1,lat0,lat1,4) gswm_mi_di_z(lon0:lon1,lat0:lat1) = fgswm(:,:,1) gswm_mi_di_t(lon0:lon1,lat0:lat1) = fgswm(:,:,2) gswm_mi_di_u(lon0:lon1,lat0:lat1) = fgswm(:,:,3) gswm_mi_di_v(lon0:lon1,lat0:lat1) = fgswm(:,:,4) case ('mi_sdi') fgswm(:,:,1) = gswm_mi_sdi_z(lon0:lon1,lat0:lat1) fgswm(:,:,2) = gswm_mi_sdi_t(lon0:lon1,lat0:lat1) fgswm(:,:,3) = gswm_mi_sdi_u(lon0:lon1,lat0:lat1) fgswm(:,:,4) = gswm_mi_sdi_v(lon0:lon1,lat0:lat1) call mp_periodic_f2d(fgswm,lon0,lon1,lat0,lat1,4) gswm_mi_sdi_z(lon0:lon1,lat0:lat1) = fgswm(:,:,1) gswm_mi_sdi_t(lon0:lon1,lat0:lat1) = fgswm(:,:,2) gswm_mi_sdi_u(lon0:lon1,lat0:lat1) = fgswm(:,:,3) gswm_mi_sdi_v(lon0:lon1,lat0:lat1) = fgswm(:,:,4) case ('nm_di') fgswm(:,:,1) = gswm_nm_di_z(lon0:lon1,lat0:lat1) fgswm(:,:,2) = gswm_nm_di_t(lon0:lon1,lat0:lat1) fgswm(:,:,3) = gswm_nm_di_u(lon0:lon1,lat0:lat1) fgswm(:,:,4) = gswm_nm_di_v(lon0:lon1,lat0:lat1) call mp_periodic_f2d(fgswm,lon0,lon1,lat0,lat1,4) gswm_nm_di_z(lon0:lon1,lat0:lat1) = fgswm(:,:,1) gswm_nm_di_t(lon0:lon1,lat0:lat1) = fgswm(:,:,2) gswm_nm_di_u(lon0:lon1,lat0:lat1) = fgswm(:,:,3) gswm_nm_di_v(lon0:lon1,lat0:lat1) = fgswm(:,:,4) case ('nm_sdi') fgswm(:,:,1) = gswm_nm_sdi_z(lon0:lon1,lat0:lat1) fgswm(:,:,2) = gswm_nm_sdi_t(lon0:lon1,lat0:lat1) fgswm(:,:,3) = gswm_nm_sdi_u(lon0:lon1,lat0:lat1) fgswm(:,:,4) = gswm_nm_sdi_v(lon0:lon1,lat0:lat1) call mp_periodic_f2d(fgswm,lon0,lon1,lat0,lat1,4) gswm_nm_sdi_z(lon0:lon1,lat0:lat1) = fgswm(:,:,1) gswm_nm_sdi_t(lon0:lon1,lat0:lat1) = fgswm(:,:,2) gswm_nm_sdi_u(lon0:lon1,lat0:lat1) = fgswm(:,:,3) gswm_nm_sdi_v(lon0:lon1,lat0:lat1) = fgswm(:,:,4) case default end select #else select case (type) case ('mi_di') call set_periodic_f2d(gswm_mi_di_z) call set_periodic_f2d(gswm_mi_di_t) call set_periodic_f2d(gswm_mi_di_u) call set_periodic_f2d(gswm_mi_di_v) case ('mi_sdi') call set_periodic_f2d(gswm_mi_sdi_z) call set_periodic_f2d(gswm_mi_sdi_t) call set_periodic_f2d(gswm_mi_sdi_u) call set_periodic_f2d(gswm_mi_sdi_v) case ('nm_di') call set_periodic_f2d(gswm_nm_di_z) call set_periodic_f2d(gswm_nm_di_t) call set_periodic_f2d(gswm_nm_di_u) call set_periodic_f2d(gswm_nm_di_v) case ('nm_sdi') call set_periodic_f2d(gswm_nm_sdi_z) call set_periodic_f2d(gswm_nm_sdi_t) call set_periodic_f2d(gswm_nm_sdi_u) call set_periodic_f2d(gswm_nm_sdi_v) case default end select #endif select case (type) case ('mi_di') ! call addfld('mi_di_z','GSWM migrating diurnal Z','cm', ! | gswm_mi_di_z(i0:i1,j0:j1),'lon',i0,i1,'lat',j0,j1,0) ! call addfld('mi_di_t','GSWM migrating diurnal TN','K', ! | gswm_mi_di_t(i0:i1,j0:j1),'lon',i0,i1,'lat',j0,j1,0) ! call addfld('mi_di_u','GSWM migrating diurnal UN','cm/s', ! | gswm_mi_di_u(i0:i1,j0:j1),'lon',i0,i1,'lat',j0,j1,0) ! call addfld('mi_di_v','GSWM migrating diurnal VN','cm/s', ! | gswm_mi_di_v(i0:i1,j0:j1),'lon',i0,i1,'lat',j0,j1,0) case ('mi_sdi') ! call addfld('mi_sdi_z','GSWM migrating semi-diurnal Z','cm', ! | gswm_mi_sdi_z(i0:i1,j0:j1),'lon',i0,i1,'lat',j0,j1,0) ! call addfld('mi_sdi_t','GSWM migrating semi-diurnal TN','K', ! | gswm_mi_sdi_t(i0:i1,j0:j1),'lon',i0,i1,'lat',j0,j1,0) ! call addfld('mi_sdi_u','GSWM migrating semi-diurnal UN','cm/s', ! | gswm_mi_sdi_u(i0:i1,j0:j1),'lon',i0,i1,'lat',j0,j1,0) ! call addfld('mi_sdi_v','GSWM migrating semi-diurnal VN','cm/s', ! | gswm_mi_sdi_v(i0:i1,j0:j1),'lon',i0,i1,'lat',j0,j1,0) case ('nm_di') ! call addfld('nm_di_z','GSWM non-migrating diurnal Z','cm', ! | gswm_nm_di_z(i0:i1,j0:j1),'lon',i0,i1,'lat',j0,j1,0) ! call addfld('nm_di_t','GSWM non-migrating diurnal TN','K', ! | gswm_nm_di_t(i0:i1,j0:j1),'lon',i0,i1,'lat',j0,j1,0) ! call addfld('nm_di_u','GSWM non-migrating diurnal UN','cm/s', ! | gswm_nm_di_u(i0:i1,j0:j1),'lon',i0,i1,'lat',j0,j1,0) ! call addfld('nm_di_v','GSWM non-migrating diurnal VN','cm/s', ! | gswm_nm_di_v(i0:i1,j0:j1),'lon',i0,i1,'lat',j0,j1,0) case ('nm_sdi') ! call addfld('gswm_nm_sdi_z','GSWM non-migrating semi-diurnal Z', ! | 'cm',gswm_nm_sdi_z(i0:i1,j0:j1),'lon',i0,i1,'lat',j0,j1,0) ! call addfld('gswm_nm_sdi_t','GSWM non-migrating semi-diurnal TN', ! | 'K',gswm_nm_sdi_t(i0:i1,j0:j1),'lon',i0,i1,'lat',j0,j1,0) ! call addfld('gswm_nm_sdi_u','GSWM non-migrating semi-diurnal UN', ! | 'cm/s',gswm_nm_sdi_u(i0:i1,j0:j1),'lon',i0,i1,'lat',j0,j1,0) ! call addfld('gswm_nm_sdi_v','GSWM non-migrating semi-diurnal VN', ! | 'cm/s',gswm_nm_sdi_v(i0:i1,j0:j1),'lon',i0,i1,'lat',j0,j1,0) case default end select end subroutine mkgswm !----------------------------------------------------------------------- function gswmdat(type,mon,ihr,fname) ! ! Return subdomain array from read arrays according to requested field, ! month, and hour. ! The 16 global fields t_mi_di, etc, are private module data, read by ! sub rdgswm at beginning of model run. ! ! Args: character(len=*),intent(in) :: type,fname integer,intent(in) :: mon,ihr ! ! Function output dimension: real :: gswmdat(lon0:lon1,lat0:lat1) ! ! Field type must be t,u,v, or z: if (fname/='t'.and.fname/='u'.and.fname/='v'.and.fname/='z') then write(6,"(/,'>>> gswmdat: unknown fname=',a)") fname write(6,"('Must be t, u, v, or z')") call shutdown('gswmdat') endif select case(trim(type)) case('mi_di') ! migrating diurnal if (fname=='t') gswmdat(:,:)= t_mi_di(lon0:lon1,lat0:lat1,mon,ihr) if (fname=='u') gswmdat(:,:)= u_mi_di(lon0:lon1,lat0:lat1,mon,ihr) if (fname=='v') gswmdat(:,:)= v_mi_di(lon0:lon1,lat0:lat1,mon,ihr) if (fname=='z') gswmdat(:,:)= z_mi_di(lon0:lon1,lat0:lat1,mon,ihr) case('mi_sdi') ! migrating semi-diurnal if (fname=='t') gswmdat(:,:)=t_mi_sdi(lon0:lon1,lat0:lat1,mon,ihr) if (fname=='u') gswmdat(:,:)=u_mi_sdi(lon0:lon1,lat0:lat1,mon,ihr) if (fname=='v') gswmdat(:,:)=v_mi_sdi(lon0:lon1,lat0:lat1,mon,ihr) if (fname=='z') gswmdat(:,:)=z_mi_sdi(lon0:lon1,lat0:lat1,mon,ihr) case('nm_di') ! non-migrating diurnal if (fname=='t') gswmdat(:,:)=t_nm_di(lon0:lon1,lat0:lat1,mon,ihr) if (fname=='u') gswmdat(:,:)=u_nm_di(lon0:lon1,lat0:lat1,mon,ihr) if (fname=='v') gswmdat(:,:)=v_nm_di(lon0:lon1,lat0:lat1,mon,ihr) if (fname=='z') gswmdat(:,:)=z_nm_di(lon0:lon1,lat0:lat1,mon,ihr) case('nm_sdi') ! non-migrating semi-diurnal if (fname=='t') gswmdat(:,:)=t_nm_sdi(lon0:lon1,lat0:lat1,mon,ihr) if (fname=='u') gswmdat(:,:)=u_nm_sdi(lon0:lon1,lat0:lat1,mon,ihr) if (fname=='v') gswmdat(:,:)=v_nm_sdi(lon0:lon1,lat0:lat1,mon,ihr) if (fname=='z') gswmdat(:,:)=z_nm_sdi(lon0:lon1,lat0:lat1,mon,ihr) ! ! Unknown type: case default write(6,"(/,'>>> gswmdat: unknown type=',a)") type call shutdown('gswmdat') end select end function gswmdat !----------------------------------------------------------------------- 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 !----------------------------------------------------------------------- end module gswm_module