! module gswm_module use params_module,only: nlon,nlat,nlonp4,nlonp2 use mpi_module,only: lon0,lon1,lat0,lat1 ! ! module used to read in output data from the GSWM model (Z,TN,UN and VN) ! Z [m] TN [K} UN,VN [m/s] ! converted for TIEGCM to Z [cm] TN [K} UN,VN [cm/s] ! subroutines rdgswmdi & rdgswmsdi ! so far the GSWM output for the whole year (one day per month) is read ! can be changed to only part of the year- however the interpolation ! subroutines getgswmdi & getgswmsdi will change too ! ! am 11/27/02: for using the double TIEGCM resolution grid together with the ! GSWM input, the GSWM grid is set to the double resolution grid and for ! single TIEGCM resolution only every second GSWM grid point is used ! implicit none ! #ifdef SUN #include "netcdf3.inc" #else #include "netcdf.inc" #endif ! ! Global gswm data read from gswm data file: ! (assume always one day of hourly values for each month and one year) integer,parameter :: | mxgswmyear = 1, ! maximum number of years of gswm data | mxgswmmonth= 12, ! maximum number of months of gswm data | mxgswmhour = 24, ! maximum number of hours of gswm data | mxgswmlon = nlon, ! maximum number of longitudinal grid points | mxgswmlat = nlat ! maximum number of latitudinal grid points integer :: | ngswmmonth, ! number of months of gswm data | ngswmtime ! number of hourly gswm data ! real :: | gswm_di_z(nlon,nlat,mxgswmmonth,mxgswmhour), ! diurnal tides data | gswm_di_tn(nlon,nlat,mxgswmmonth,mxgswmhour), ! diurnal tides data | gswm_di_un(nlon,nlat,mxgswmmonth,mxgswmhour), ! diurnal tides data | gswm_di_vn(nlon,nlat,mxgswmmonth,mxgswmhour), ! diurnal tides data | gswm_sdi_z(nlon,nlat,mxgswmmonth,mxgswmhour), ! semidiurnal tides data | gswm_sdi_tn(nlon,nlat,mxgswmmonth,mxgswmhour), ! semidiurnal tides data | gswm_sdi_un(nlon,nlat,mxgswmmonth,mxgswmhour), ! semidiurnal tides data | gswm_sdi_vn(nlon,nlat,mxgswmmonth,mxgswmhour), ! semidiurnal tides data | gswm_nmidi_z(nlon,nlat,mxgswmmonth,mxgswmhour), ! nonmigrating diurnal tides data | gswm_nmidi_tn(nlon,nlat,mxgswmmonth,mxgswmhour), ! nonmigrating diurnal tides data | gswm_nmidi_un(nlon,nlat,mxgswmmonth,mxgswmhour), ! nonmigrating diurnal tides data | gswm_nmidi_vn(nlon,nlat,mxgswmmonth,mxgswmhour), ! nonmigrating diurnal tides data | gswm_nmisdi_z(nlon,nlat,mxgswmmonth,mxgswmhour), ! nonmigrating semidiurnal tides data | gswm_nmisdi_tn(nlon,nlat,mxgswmmonth,mxgswmhour),! nonmigrating semidiurnal tides data | gswm_nmisdi_un(nlon,nlat,mxgswmmonth,mxgswmhour),! nonmigrating semidiurnal tides data | gswm_nmisdi_vn(nlon,nlat,mxgswmmonth,mxgswmhour),! nonmigrating semidiurnal tides data ! | zdi_gswm(nlonp4,nlat),tndi_gswm(nlonp4,nlat), ! diurnal tides boundary | undi_gswm(nlonp4,nlat),vndi_gswm(nlonp4,nlat), ! diurnal tides boundary | zsdi_gswm(nlonp4,nlat),tnsdi_gswm(nlonp4,nlat), ! semidiurnal tides boundary | unsdi_gswm(nlonp4,nlat),vnsdi_gswm(nlonp4,nlat), ! semidiurnal tides boundary | znmidi_gswm(nlonp4,nlat),tnnmidi_gswm(nlonp4,nlat), ! nonmigrating diurnal tides boundary | unnmidi_gswm(nlonp4,nlat),vnnmidi_gswm(nlonp4,nlat), ! nonmigrating diurnal tides boundary | znmisdi_gswm(nlonp4,nlat),tnnmisdi_gswm(nlonp4,nlat), ! nonmigrating semidiurnal tides boundary | unnmisdi_gswm(nlonp4,nlat),vnnmisdi_gswm(nlonp4,nlat) ! nonmigrating semidiurnal tides boundary ! contains !----------------------------------------------------------------------- subroutine rdgswmdi ! ! Obtain and read gswm_di_ncfile netcdf data file containing T,UN,VN and Z. ! These data is generated by the ouput of the GSWM. ! use input_module,only: tempdir,gswm_di_ncfile use nchist_module,only:nc_open,nc_close,handle_ncerr ! ! Local: character(len=80) :: dskfile integer :: istat,ncid,i,j,ier,i2,j2 integer :: id_nmonth,id_ntime,id_nlon,id_nlat,idv_z,idv_tn, | idv_un,idv_vn,ngswmlon,ngswmlat,lonbeg,lonend,latbeg, | lonbeg_dhz,latbeg_dhz real :: | gswm_z(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour), ! diurnal tides data org.grid | gswm_tn(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour), ! diurnal tides data org.grid | gswm_un(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour), ! diurnal tides data org.grid | gswm_vn(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour) ! diurnal tides data org.grid ! ! Acquire mss file: call mkdiskflnm(gswm_di_ncfile,dskfile) call getms(gswm_di_ncfile,dskfile,tempdir,' ') ! write(6,"(/,72('-'))") write(6,"('RDGSWMDI: read GSWM diurnal tidal data file:')") ! ! Open netcdf file: call nc_open(ncid,dskfile,'OLD','READ') if (ncid==0) then write(6,"(/,'>>> rdgswmdi: error opening netcdf gswm ', | 'diurnal file ',a)") trim(dskfile) stop 'rdgswmdi' endif ! ! Get nmonth dimension: istat = nf_inq_dimid(ncid,'nmonth',id_nmonth) istat = nf_inq_dimlen(ncid,id_nmonth,ngswmmonth) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmdi: Error getting nmonth dimension') if (ngswmmonth > mxgswmmonth) then write(6,"(/,'>>> rdgswmdi: need to increase mxgswmmonth: ', | 'mxgswmmonth=',i8,' ngswmmonth=',i8)") mxgswmmonth, | ngswmmonth stop 'rdgswmdi' endif ! ! Get ntime dimension: istat = nf_inq_dimid(ncid,'time',id_ntime) istat = nf_inq_dimlen(ncid,id_ntime,ngswmtime) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmdi: Error getting time dimension') if (ngswmtime > mxgswmhour) then write(6,"(/,'>>> rdgswmdi: input hourly data: ', | 'mxgswmhour=',i8,' ngswmtime=',i8)") mxgswmhour, | ngswmtime stop 'rdgswmdi' endif ! ! Get nlon dimension: istat = nf_inq_dimid(ncid,'lon',id_nlon) istat = nf_inq_dimlen(ncid,id_nlon,ngswmlon) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmdi: Error getting lon dimension') if (ngswmlon /= mxgswmlon) then write(6,"(/,'>>> rdgswmdi: longitudinal dim does not fit: ', | 'nlon(GSWM)=',i8,' mxgswmlon=',i8)") ngswmlon, | mxgswmlon stop 'rdgswmdi' endif ! ! Get nlat dimension: istat = nf_inq_dimid(ncid,'lat',id_nlat) istat = nf_inq_dimlen(ncid,id_nlat,ngswmlat) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmdi: Error getting lat dimension') if (ngswmlat /= mxgswmlat) then write(6,"(/,'>>> rdgswmdi: latitudinal dim does not fit: ', | 'nlat(GSWM)=',i8,' mxgswmlat=',i8)") ngswmlat, | mxgswmlat stop 'rdgswmdi' endif ! ! read whole domain ! Get Z (geopotential height perturbation [m]: istat = nf_inq_varid(ncid,'Z',idv_z) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmdi: Error getting Z var id') istat = nf_get_var_double(ncid,idv_z,gswm_z) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswm: Error getting variable Z') gswm_z = gswm_z*100. ! convert to cm ! ! Get TN (neutral temperature perturbation [deg K]: istat = nf_inq_varid(ncid,'TN',idv_tn) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmdi: Error getting TN var id') istat = nf_get_var_double(ncid,idv_tn,gswm_tn) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswm: Error getting variable TN') ! ! Get UN (neutral zonal wind perturbation [m/s]: istat = nf_inq_varid(ncid,'UN',idv_un) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmdi: Error getting UN var id') istat = nf_get_var_double(ncid,idv_un,gswm_un) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswm: Error getting variable UN') gswm_un = gswm_un*100. ! convert to cm/s ! ! Get VN (neutral meridional wind perturbation [m/s]: istat = nf_inq_varid(ncid,'VN',idv_vn) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmdi: Error getting VN var id') istat = nf_get_var_double(ncid,idv_vn,gswm_vn) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswm: Error getting variable VN') gswm_vn = gswm_vn*100. ! convert to cm/s ! ! Close the file: call nc_close(ncid) write(6,"('Completed read from GSWMDI data file ',a)") | trim(dskfile) ! ! check for the TIEGCM dimension of the grid ! if necessary copy data to right sized arrays ! copy into subdomain ! lonbeg = lon0 if (lon0==1) lonbeg = 3 lonend = lon1 if (lon1==nlonp4) lonend = lon1-2 if(nlat == mxgswmlat) then ! same dimension nlon=144,nlat=72 do i = lonbeg,lonend do j = lat0,lat1 gswm_di_z(i,j,:,:) = gswm_z(i-2,j,:,:) gswm_di_tn(i,j,:,:) = gswm_tn(i-2,j,:,:) gswm_di_un(i,j,:,:) = gswm_un(i-2,j,:,:) gswm_di_vn(i,j,:,:) = gswm_vn(i-2,j,:,:) enddo enddo endif write(6,"(72('-'),/)") end subroutine rdgswmdi !----------------------------------------------------------------------- subroutine rdgswmsdi ! ! Obtain and read gswm_sdi_ncfile netcdf data file containing T,UN,VN and Z. ! These data is generated by the ouput of the GSWM. ! use input_module,only: tempdir,gswm_sdi_ncfile use nchist_module,only:nc_open,nc_close,handle_ncerr ! ! Local: character(len=80) :: dskfile integer :: istat,ncid,i,j,ier,i2,j2 integer :: id_nmonth,id_ntime,id_nlon,id_nlat,idv_z,idv_tn, | idv_un,idv_vn,ngswmlon,ngswmlat,lonbeg,lonend,latbeg, | lonbeg_dhz,latbeg_dhz real :: | gswm_z(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour), ! semidiurnal tides data org.grid | gswm_tn(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour), ! semidiurnal tides data org.grid | gswm_un(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour), ! semidiurnal tides data org.grid | gswm_vn(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour) ! semidiurnal tides data org.grid ! ! Acquire mss file: call mkdiskflnm(gswm_sdi_ncfile,dskfile) call getms(gswm_sdi_ncfile,dskfile,tempdir,' ') ! write(6,"(/,72('-'))") write(6,"('RDGSWMSDI: read GSWM semidiurnal tidal data file:')") ! ! Open netcdf file: call nc_open(ncid,dskfile,'OLD','READ') if (ncid==0) then write(6,"(/,'>>> rdgswmdi: error opening netcdf gswm ', | 'semidiurnal file ',a)") trim(dskfile) stop 'rdgswmsdi' endif ! ! Get nmonth dimension: istat = nf_inq_dimid(ncid,'nmonth',id_nmonth) istat = nf_inq_dimlen(ncid,id_nmonth,ngswmmonth) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmsdi: Error getting nmonth dimension') if (ngswmmonth > mxgswmmonth) then write(6,"(/,'>>> rdgswmsdi: need to increase mxgswmmonth: ', | 'mxgswmmonth=',i8,' ngswmmonth=',i8)") mxgswmmonth, | ngswmmonth stop 'rdgswmsdi' endif ! ! Get ntime dimension: istat = nf_inq_dimid(ncid,'time',id_ntime) istat = nf_inq_dimlen(ncid,id_ntime,ngswmtime) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmsdi: Error getting time dimension') if (ngswmtime > mxgswmhour) then write(6,"(/,'>>> rdgswmdi: input hourly data: ', | 'mxgswmhour=',i8,' ngswmtime=',i8)") mxgswmhour, | ngswmtime stop 'rdgswmsdi' endif ! ! Get nlon dimension: istat = nf_inq_dimid(ncid,'lon',id_nlon) istat = nf_inq_dimlen(ncid,id_nlon,ngswmlon) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmsdi: Error getting lon dimension') if (ngswmlon /= mxgswmlon) then write(6,"(/,'>>> rdgswmsdi: longitudinal dim does not fit: ', | 'nlon(GSWM)=',i8,' mxgswmlon=',i8)") ngswmlon, | mxgswmlon stop 'rdgswmsdi' endif ! ! Get nlat dimension: istat = nf_inq_dimid(ncid,'lat',id_nlat) istat = nf_inq_dimlen(ncid,id_nlat,ngswmlat) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmsdi: Error getting lat dimension') if (ngswmlat /= mxgswmlat) then write(6,"(/,'>>> rdgswmsdi: latitudinal dim does not fit: ', | 'nlat(GSWM)=',i8,' mxgswmlat=',i8)") ngswmlat, | mxgswmlat stop 'rdgswmsdi' endif ! ! read whole domain ! Get Z (geopotential height perturbation [m]: istat = nf_inq_varid(ncid,'Z',idv_z) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmsdi: Error getting Z var id') istat = nf_get_var_double(ncid,idv_z,gswm_z) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmsdi: Error getting variable Z') gswm_z = gswm_z*100. ! convert to cm ! ! Get TN (neutral temperature perturbation [deg K]: istat = nf_inq_varid(ncid,'TN',idv_tn) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmdi: Error getting TN var id') istat = nf_get_var_double(ncid,idv_tn,gswm_tn) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmsdi: Error getting variable TN') ! ! Get UN (neutral zonal wind perturbation [m/s]: istat = nf_inq_varid(ncid,'UN',idv_un) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmdi: Error getting UN var id') istat = nf_get_var_double(ncid,idv_un,gswm_un) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmsdi: Error getting variable UN') gswm_un = gswm_un*100. ! convert to cm/s ! ! Get VN (neutral meridional wind perturbation [m/s]: istat = nf_inq_varid(ncid,'VN',idv_vn) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmdi: Error getting VN var id') istat = nf_get_var_double(ncid,idv_vn,gswm_vn) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmsdi: Error getting variable VN') gswm_vn = gswm_vn*100. ! convert to cm/s ! ! Close the file: call nc_close(ncid) write(6,"('Completed read from GSWMSDI data file ',a)") | trim(dskfile) ! check for the TIEGCM dimension of the grid ! if necessary copy data to right sized arrays ! copy into subdomain ! lonbeg = lon0 if (lon0==1) lonbeg = 3 lonend = lon1 if (lon1==nlonp4) lonend = lon1-2 if(nlat == mxgswmlat) then ! same dimension nlon=144,nlat=72 do i = lonbeg,lonend do j = lat0,lat1 gswm_sdi_z(i,j,:,:) = gswm_z(i-2,j,:,:) gswm_sdi_tn(i,j,:,:) = gswm_tn(i-2,j,:,:) gswm_sdi_un(i,j,:,:) = gswm_un(i-2,j,:,:) gswm_sdi_vn(i,j,:,:) = gswm_vn(i-2,j,:,:) enddo enddo endif write(6,"(72('-'),/)") end subroutine rdgswmsdi !----------------------------------------------------------------------- subroutine rdgswmnmidi ! ! Obtain and read gswm_nmidi_ncfile netcdf data file containing T,UN,VN and Z. ! These data is generated by the ouput of the GSWM. ! use input_module,only: tempdir,gswm_nmidi_ncfile use nchist_module,only:nc_open,nc_close,handle_ncerr ! ! Local: character(len=80) :: dskfile integer :: istat,ncid,i,j,ier,i2,j2 integer :: id_nmonth,id_ntime,id_nlon,id_nlat,idv_z,idv_tn, | idv_un,idv_vn,ngswmlon,ngswmlat,lonbeg,lonend,latbeg, | lonbeg_dhz,latbeg_dhz real :: ! org. GSWM grid | gswm_z(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour), ! nonmigrating diurnal tides data | gswm_tn(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour), ! nonmigrating diurnal tides data | gswm_un(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour), ! nonmigrating diurnal tides data | gswm_vn(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour) ! nonmigrating diurnal tides data ! ! Acquire mss file: call mkdiskflnm(gswm_nmidi_ncfile,dskfile) call getms(gswm_nmidi_ncfile,dskfile,tempdir,' ') ! write(6,"(/,72('-'))") write(6,"('RDGSWMNMIDI: read GSWM nonmigrating diurnal', | ' tidal data file:')") ! ! Open netcdf file: call nc_open(ncid,dskfile,'OLD','READ') if (ncid==0) then write(6,"(/,'>>> rdgswmnmidi: error opening netcdf gswm ', | 'nonmigrating diurnal file ',a)") trim(dskfile) stop 'rdgswmnmidi' endif ! ! Get nmonth dimension: istat = nf_inq_dimid(ncid,'nmonth',id_nmonth) istat = nf_inq_dimlen(ncid,id_nmonth,ngswmmonth) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmidi: Error getting nmonth dimension') if (ngswmmonth > mxgswmmonth) then write(6,"(/,'>>> rdgswmnmidi: need to increase mxgswmmonth: ', | 'mxgswmmonth=',i8,' ngswmmonth=',i8)") mxgswmmonth, | ngswmmonth stop 'rdgswmnmidi' endif ! ! Get ntime dimension: istat = nf_inq_dimid(ncid,'time',id_ntime) istat = nf_inq_dimlen(ncid,id_ntime,ngswmtime) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmidi: Error getting time dimension') if (ngswmtime > mxgswmhour) then write(6,"(/,'>>> rdgswmnmidi: input hourly data: ', | 'mxgswmhour=',i8,' ngswmtime=',i8)") mxgswmhour, | ngswmtime stop 'rdgswmnmidi' endif ! ! Get nlon dimension: istat = nf_inq_dimid(ncid,'lon',id_nlon) istat = nf_inq_dimlen(ncid,id_nlon,ngswmlon) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmidi: Error getting lon dimension') if (ngswmlon /= mxgswmlon) then write(6,"(/,'>>> rdgswmnmidi: longitudinal dim does not fit: ', | 'nlon(GSWM)=',i8,' mxgswmlon=',i8)") ngswmlon, | mxgswmlon stop 'rdgswmnmidi' endif ! ! Get nlat dimension: istat = nf_inq_dimid(ncid,'lat',id_nlat) istat = nf_inq_dimlen(ncid,id_nlat,ngswmlat) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmidi: Error getting lat dimension') if (ngswmlat /= mxgswmlat) then write(6,"(/,'>>> rdgswmnmidi: latitudinal dim does not fit: ', | 'nlat(GSWM)=',i8,' mxgswmlat=',i8)") ngswmlat, | mxgswmlat stop 'rdgswmnmidi' endif ! ! read whole domain ! Get Z (geopotential height perturbation [m]: istat = nf_inq_varid(ncid,'Z',idv_z) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmidi: Error getting Z var id') istat = nf_get_var_double(ncid,idv_z,gswm_z) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmidi: Error getting variable Z') gswm_z = gswm_z*100. ! convert to cm ! ! Get TN (neutral temperature perturbation [deg K]: istat = nf_inq_varid(ncid,'TN',idv_tn) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmidi: Error getting TN var id') istat = nf_get_var_double(ncid,idv_tn,gswm_tn) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmidi: Error getting variable TN') ! ! Get UN (neutral zonal wind perturbation [m/s]: istat = nf_inq_varid(ncid,'UN',idv_un) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmidi: Error getting UN var id') istat = nf_get_var_double(ncid,idv_un,gswm_un) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmidi: Error getting variable UN') gswm_un = gswm_un*100. ! convert to cm/s ! ! Get VN (neutral meridional wind perturbation [m/s]: istat = nf_inq_varid(ncid,'VN',idv_vn) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmidi: Error getting VN var id') istat = nf_get_var_double(ncid,idv_vn,gswm_vn) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmidi: Error getting variable VN') gswm_vn = gswm_vn*100. ! convert to cm/s ! ! Close the file: call nc_close(ncid) write(6,"('Completed read from GSWMNMIDI data file ',a)") | trim(dskfile) ! check for the TIEGCM dimension of the grid ! if necessary copy data to right sized arrays ! copy into subdomain ! lonbeg = lon0 if (lon0==1) lonbeg = 3 lonend = lon1 if (lon1==nlonp4) lonend = lon1-2 if(nlat == mxgswmlat) then ! same dimension nlon=144,nlat=72 do i = lonbeg,lonend do j = lat0,lat1 gswm_nmidi_z(i,j,:,:) = gswm_z(i-2,j,:,:) gswm_nmidi_tn(i,j,:,:) = gswm_tn(i-2,j,:,:) gswm_nmidi_un(i,j,:,:) = gswm_un(i-2,j,:,:) gswm_nmidi_vn(i,j,:,:) = gswm_vn(i-2,j,:,:) enddo enddo endif write(6,"(72('-'),/)") end subroutine rdgswmnmidi !----------------------------------------------------------------------- subroutine rdgswmnmisdi ! ! Obtain and read gswm_nmisdi_ncfile netcdf data file containing T,UN,VN and Z. ! These data is generated by the ouput of the GSWM. ! use input_module,only: tempdir,gswm_nmisdi_ncfile use nchist_module,only:nc_open,nc_close,handle_ncerr ! ! Local: character(len=80) :: dskfile integer :: istat,ncid,i,j,ier,i2,j2 integer :: id_nmonth,id_ntime,id_nlon,id_nlat,idv_z,idv_tn, | idv_un,idv_vn,ngswmlon,ngswmlat,lonbeg,lonend,latbeg, | lonbeg_dhz,latbeg_dhz real :: | gswm_z(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour), ! nonmigrating semidiurnal tides data | gswm_tn(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour),! nonmigrating semidiurnal tides data | gswm_un(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour),! nonmigrating semidiurnal tides data | gswm_vn(mxgswmlon,mxgswmlat,mxgswmmonth,mxgswmhour)! nonmigrating semidiurnal tides data ! ! Acquire mss file: call mkdiskflnm(gswm_nmisdi_ncfile,dskfile) call getms(gswm_nmisdi_ncfile,dskfile,tempdir,' ') ! write(6,"(/,72('-'))") write(6,"('RDGSWMNMISDI: read GSWM nonmigrating semidiurnal', | ' tidal data file:')") ! ! Open netcdf file: call nc_open(ncid,dskfile,'OLD','READ') if (ncid==0) then write(6,"(/,'>>> rdgswmnmisdi: error opening netcdf gswm ', | 'nonmigrating semidiurnal file ',a)") trim(dskfile) stop 'rdgswmnmisdi' endif ! ! Get nmonth dimension: istat = nf_inq_dimid(ncid,'nmonth',id_nmonth) istat = nf_inq_dimlen(ncid,id_nmonth,ngswmmonth) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmisdi: Error getting nmonth dimension') if (ngswmmonth > mxgswmmonth) then write(6,"(/,'>>> rdgswmnmisdi: need to increase mxgswmmonth: ', | 'mxgswmmonth=',i8,' ngswmmonth=',i8)") mxgswmmonth, | ngswmmonth stop 'rdgswmnmisdi' endif ! ! Get ntime dimension: istat = nf_inq_dimid(ncid,'time',id_ntime) istat = nf_inq_dimlen(ncid,id_ntime,ngswmtime) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmisdi: Error getting time dimension') if (ngswmtime > mxgswmhour) then write(6,"(/,'>>> rdgswmnmisdi: input hourly data: ', | 'mxgswmhour=',i8,' ngswmtime=',i8)") mxgswmhour, | ngswmtime stop 'rdgswmnmisdi' endif ! ! Get nlon dimension: istat = nf_inq_dimid(ncid,'lon',id_nlon) istat = nf_inq_dimlen(ncid,id_nlon,ngswmlon) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmisdi: Error getting lon dimension') if (ngswmlon /= mxgswmlon) then write(6,"(/,'>>> rdgswmnmisdi: longitudinal dim does not ', | 'fit: nlon(GSWM)=',i8,' mxgswmlon=',i8)") ngswmlon, | mxgswmlon stop 'rdgswmnmisdi' endif ! ! Get nlat dimension: istat = nf_inq_dimid(ncid,'lat',id_nlat) istat = nf_inq_dimlen(ncid,id_nlat,ngswmlat) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmisdi: Error getting lat dimension') if (ngswmlat /= mxgswmlat) then write(6,"(/,'>>> rdgswmnmisdi: latitudinal dim does not fit: ', | 'nlat(GSWM)=',i8,' mxgswmlat=',i8)") ngswmlat, | mxgswmlat stop 'rdgswmnmisdi' endif ! ! read whole domain ! Get Z (geopotential height perturbation [m]: istat = nf_inq_varid(ncid,'Z',idv_z) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmisdi: Error getting Z var id') istat = nf_get_var_double(ncid,idv_z,gswm_z) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmisdi: Error getting variable Z') gswm_z = gswm_z*100. ! convert to cm ! ! Get TN (neutral temperature perturbation [deg K]: istat = nf_inq_varid(ncid,'TN',idv_tn) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmisdi: Error getting TN var id') istat = nf_get_var_double(ncid,idv_tn,gswm_tn) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmisdi: Error getting variable TN') ! ! Get UN (neutral zonal wind perturbation [m/s]: istat = nf_inq_varid(ncid,'UN',idv_un) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmisdi: Error getting UN var id') istat = nf_get_var_double(ncid,idv_un,gswm_un) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmisdi: Error getting variable UN') gswm_un = gswm_un*100. ! convert to cm/s ! ! Get VN (neutral meridional wind perturbation [m/s]: istat = nf_inq_varid(ncid,'VN',idv_vn) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmisdi: Error getting VN var id') istat = nf_get_var_double(ncid,idv_vn,gswm_vn) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'rdgswmnmisdi: Error getting variable VN') gswm_vn = gswm_vn*100. ! convert to cm/s ! ! Close the file: call nc_close(ncid) write(6,"('Completed read from GSWMNMISDI data file ',a)") | trim(dskfile) ! check for the TIEGCM dimension of the grid ! if necessary copy data to right sized arrays ! lonbeg = lon0 if (lon0==1) lonbeg = 3 lonend = lon1 if (lon1==nlonp4) lonend = lon1-2 if(nlat == mxgswmlat) then ! same dimension nlon=144,nlat=72 do i = lonbeg,lonend do j = lat0,lat1 gswm_nmisdi_z(i,j,:,:) = gswm_z(i-2,j,:,:) gswm_nmisdi_tn(i,j,:,:) = gswm_tn(i-2,j,:,:) gswm_nmisdi_un(i,j,:,:) = gswm_un(i-2,j,:,:) gswm_nmisdi_vn(i,j,:,:) = gswm_vn(i-2,j,:,:) enddo enddo endif write(6,"(72('-'),/)") end subroutine rdgswmnmisdi !----------------------------------------------------------------------- subroutine getgswmdi(iyear,iday,iutsec,iprint) ! ! Use gswm data read from gswm_di_ncfile file to return perturbation in Z ! and T at current date and time. It is assumed that the gswm_di_ncfile ! 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 ! ! Args: integer,intent(in) :: iyear,iday,iutsec,iprint ! ! Local: integer :: i,j,ndmon_nl(13),ndmon_lp(13),nonleap, | mon_cur,mon_nxt,ihr_cur,ihr_nxt,nointp=0, | lonbeg,lonend real :: difsec,difday,sec,secint,dayint, | z_curmo(lon0:lon1,lat0:lat1), | z_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation | tn_curmo(lon0:lon1,lat0:lat1), | tn_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation | un_curmo(lon0:lon1,lat0:lat1), | un_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation | vn_curmo(lon0:lon1,lat0:lat1), | vn_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation | z_tmp(lon0:lon1,lat0:lat1), | tn_tmp(lon0:lon1,lat0:lat1), ! tmp | un_tmp(lon0:lon1,lat0:lat1), | vn_tmp(lon0:lon1,lat0:lat1) ! tmp ! ! 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 ! if (iprint > 0) then write(6,"(/,72('-'))") write(6,"('GETGSWMDI: get diurnal tidal perturbation from ', | 'database:')") write(6,"('Initial requested iyear=',i4,' iday=',i3,' iutsec=', | i10)") iyear,iday,iutsec endif ! ! 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 if(ihr_nxt == 24 ) ihr_nxt = 1 ! ! GSWM hourly data from 0 UT (i=1) to 23 UT (i=24) resp. 0 utsec ! ! Check times: ! if (modeltime(3).lt.0.or.modeltime(3).gt.59) then write(6,"('>>> getgswmdi: bad imin=',3i4)") modeltime(3) stop 'getgswmdi' endif if (modeltime(4).lt.0.or.modeltime(4).gt.59) then write(6,"('>>> getgswmdi: bad isec=',3i4)") modeltime(4) stop 'getgswmdi' endif ! Interpolate to modeltime iutsec for each month: mon_cur mon_nxt if (iutsec==0) then ! no interpolation z_curmo(:,:) =gswm_di_z (lon0:lon1,lat0:lat1,mon_cur,ihr_cur+1) ! current month tn_curmo(:,:)=gswm_di_tn(lon0:lon1,lat0:lat1,mon_cur,ihr_cur+1) un_curmo(:,:)=gswm_di_un(lon0:lon1,lat0:lat1,mon_cur,ihr_cur+1) vn_curmo(:,:)=gswm_di_vn(lon0:lon1,lat0:lat1,mon_cur,ihr_cur+1) z_nxtmo(:,:) =gswm_di_z (lon0:lon1,lat0:lat1,mon_nxt,ihr_cur+1) ! next month tn_nxtmo(:,:)=gswm_di_tn(lon0:lon1,lat0:lat1,mon_nxt,ihr_cur+1) un_nxtmo(:,:)=gswm_di_un(lon0:lon1,lat0:lat1,mon_nxt,ihr_cur+1) vn_nxtmo(:,:)=gswm_di_vn(lon0:lon1,lat0:lat1,mon_nxt,ihr_cur+1) else difsec = 60.*60. ! difference in sec between ihr_cur and ihr_nxt secint = float(modeltime(3)*60 + modeltime(4)) ! interpolation time [sec] call timeliper_array(gswm_di_z(lon0:lon1,lat0:lat1,mon_cur, ! cur. month | ihr_cur+1),gswm_di_z(lon0:lon1,lat0:lat1,mon_cur,ihr_nxt+1), | difsec,secint,z_curmo,1) call timeliper_array(gswm_di_tn(lon0:lon1,lat0:lat1,mon_cur, ! cur. month | ihr_cur+1),gswm_di_tn(lon0:lon1,lat0:lat1,mon_cur,ihr_nxt+1), | difsec,secint,tn_curmo,1) call timeliper_array(gswm_di_un(lon0:lon1,lat0:lat1,mon_cur, ! cur. month | ihr_cur+1),gswm_di_un(lon0:lon1,lat0:lat1,mon_cur,ihr_nxt+1), | difsec,secint,un_curmo,1) call timeliper_array(gswm_di_vn(lon0:lon1,lat0:lat1,mon_cur, ! cur. month | ihr_cur+1),gswm_di_vn(lon0:lon1,lat0:lat1,mon_cur,ihr_nxt+1), | difsec,secint,vn_curmo,1) call timeliper_array(gswm_di_z(lon0:lon1,lat0:lat1,mon_nxt, ! next month | ihr_cur+1),gswm_di_z(lon0:lon1,lat0:lat1,mon_nxt,ihr_nxt+1), | difsec,secint,z_nxtmo,1) call timeliper_array(gswm_di_tn(lon0:lon1,lat0:lat1,mon_nxt, ! next month | ihr_cur+1),gswm_di_tn(lon0:lon1,lat0:lat1,mon_nxt,ihr_nxt+1), | difsec,secint,tn_nxtmo,1) call timeliper_array(gswm_di_un(lon0:lon1,lat0:lat1,mon_nxt, ! next month | ihr_cur+1),gswm_di_un(lon0:lon1,lat0:lat1,mon_nxt,ihr_nxt+1), | difsec,secint,un_nxtmo,1) call timeliper_array(gswm_di_vn(lon0:lon1,lat0:lat1,mon_nxt, ! next month | ihr_cur+1),gswm_di_vn(lon0:lon1,lat0:lat1,mon_nxt,ihr_nxt+1), | difsec,secint,vn_nxtmo,1) endif ! ! Interpolate to modelday (iday) ! ! Check if interpolation 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 interpolation is necessary calculated the 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) if (nointp==0) then ! interpolation call timeliper_array(z_curmo,z_nxtmo,difday,dayint, | z_tmp,1) call timeliper_array(tn_curmo,tn_nxtmo,difday,dayint, | tn_tmp,1) call timeliper_array(un_curmo,un_nxtmo,difday,dayint, | un_tmp,1) call timeliper_array(vn_curmo,vn_nxtmo,difday,dayint, | vn_tmp,1) elseif(nointp == 1) then ! no interpolation same day as current month zdi_gswm (lon0:lon1,lat0:lat1) = z_curmo tndi_gswm(lon0:lon1,lat0:lat1) = tn_curmo undi_gswm(lon0:lon1,lat0:lat1) = un_curmo vndi_gswm(lon0:lon1,lat0:lat1) = vn_curmo goto 30 elseif(nointp == 2) then ! no interpolation same day as next month zdi_gswm (lon0:lon1,lat0:lat1) = z_nxtmo tndi_gswm(lon0:lon1,lat0:lat1) = tn_nxtmo undi_gswm(lon0:lon1,lat0:lat1) = un_nxtmo vndi_gswm(lon0:lon1,lat0:lat1) = vn_nxtmo goto 30 endif ! zdi_gswm (lon0:lon1,lat0:lat1) = z_tmp tndi_gswm(lon0:lon1,lat0:lat1) = tn_tmp undi_gswm(lon0:lon1,lat0:lat1) = un_tmp vndi_gswm(lon0:lon1,lat0:lat1) = vn_tmp 30 continue ! ! do mpi periodic points exchange for gswm with f2d(:) ! Moved here from sub getgswm* because mpi necessary when f2d data ! is allocated only for task-local subdomain block. ! lons 1,2 <- nlonp4-3,nlonp4-2 and nlonp4-1,nlonp4 <- 3,4 #ifdef MPI call mp_periodic_f2d(zdi_gswm(lon0:lon1,lat0:lat1), | lon0,lon1,lat0,lat1) call mp_periodic_f2d(tndi_gswm(lon0:lon1,lat0:lat1), | lon0,lon1,lat0,lat1) call mp_periodic_f2d(undi_gswm(lon0:lon1,lat0:lat1), | lon0,lon1,lat0,lat1) call mp_periodic_f2d(vndi_gswm(lon0:lon1,lat0:lat1), | lon0,lon1,lat0,lat1) # else set_periodic_f2d(zdi_gswm) set_periodic_f2d(tndi_gswm) set_periodic_f2d(undi_gswm) set_periodic_f2d(vndi_gswm) #endif if (iprint > 0) then write(6,"('GSWM diurnal tidal pertubation interpolated to ', | 'date and time')") endif end subroutine getgswmdi !----------------------------------------------------------------------- subroutine getgswmsdi(iyear,iday,iutsec,iprint) ! ! Use gswm data read from gswm_sdi_ncfile file to return perturbation in Z ! and T at current date and time. It is assumed that the gswm_sdi_ncfile ! 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 ! ! Args: integer,intent(in) :: iyear,iday,iutsec,iprint ! ! Local: integer :: i,j,ndmon_nl(13),ndmon_lp(13),nonleap, | mon_cur,mon_nxt,ihr_cur,ihr_nxt,nointp=0, | lonbeg,lonend real :: difsec,difday,sec,secint,dayint, | z_curmo(lon0:lon1,lat0:lat1), | z_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation | tn_curmo(lon0:lon1,lat0:lat1), | tn_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation | un_curmo(lon0:lon1,lat0:lat1), | un_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation | vn_curmo(lon0:lon1,lat0:lat1), | vn_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation | z_tmp(lon0:lon1,lat0:lat1), | tn_tmp(lon0:lon1,lat0:lat1), ! tmp | un_tmp(lon0:lon1,lat0:lat1), | vn_tmp(lon0:lon1,lat0:lat1) ! tmp ! ! 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 ! if (iprint > 0) then write(6,"(/,72('-'))") write(6,"('GETGSWMSDI: get semidiurnal tidal perturbation ', | 'from database:')") write(6,"('Initial requested iyear=',i4,' iday=',i3,' iutsec=', | i10)") iyear,iday,iutsec endif ! ! 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 if(ihr_nxt == 24 ) ihr_nxt = 1 ! ! GSWM hourly data from 0 UT (i=1) to 23 UT (i=24) resp. 0 utsec ! ! Check times: seems not to be necessary ! if (modeltime(3).lt.0.or.modeltime(3).gt.59) then write(6,"('>>> getgswmsdi: bad imin=',3i4)") modeltime(3) stop 'getgswmdi' endif if (modeltime(4).lt.0.or.modeltime(4).gt.59) then write(6,"('>>> getgswmsdi: bad isec=',3i4)") modeltime(4) stop 'getgswmdi' endif ! Interpolate to modeltime iutsec for each month: mon_cur mon_nxt if (iutsec==0) then ! no interpolation z_curmo(:,:) =gswm_sdi_z (lon0:lon1,lat0:lat1, ! current month | mon_cur,ihr_cur+1) tn_curmo(:,:)=gswm_sdi_tn(lon0:lon1,lat0:lat1, | mon_cur,ihr_cur+1) un_curmo(:,:)=gswm_sdi_un(lon0:lon1,lat0:lat1, | mon_cur,ihr_cur+1) vn_curmo(:,:)=gswm_sdi_vn(lon0:lon1,lat0:lat1, | mon_cur,ihr_cur+1) z_nxtmo(:,:) =gswm_sdi_z (lon0:lon1,lat0:lat1, ! next month | mon_nxt,ihr_cur+1) tn_nxtmo(:,:)=gswm_sdi_tn(lon0:lon1,lat0:lat1, | mon_nxt,ihr_cur+1) un_nxtmo(:,:)=gswm_sdi_un(lon0:lon1,lat0:lat1, | mon_nxt,ihr_cur+1) vn_nxtmo(:,:)=gswm_sdi_vn(lon0:lon1,lat0:lat1, | mon_nxt,ihr_cur+1) else difsec = 60.*60. ! difference in sec between ihr_cur and ihr_nxt secint = float(modeltime(3)*60 + modeltime(4)) ! interpolation time [sec] call timeliper_array(gswm_sdi_z(lon0:lon1,lat0:lat1,mon_cur, ! cur. month | ihr_cur+1),gswm_sdi_z(lon0:lon1,lat0:lat1,mon_cur,ihr_nxt+1), | difsec,secint,z_curmo,1) call timeliper_array(gswm_sdi_tn(lon0:lon1,lat0:lat1,mon_cur, ! cur. month | ihr_cur+1),gswm_sdi_tn(lon0:lon1,lat0:lat1,mon_cur,ihr_nxt+1), | difsec,secint,tn_curmo,1) call timeliper_array(gswm_sdi_un(lon0:lon1,lat0:lat1,mon_cur, ! cur. month | ihr_cur+1),gswm_sdi_un(lon0:lon1,lat0:lat1,mon_cur,ihr_nxt+1), | difsec,secint,un_curmo,1) call timeliper_array(gswm_sdi_vn(lon0:lon1,lat0:lat1,mon_cur, ! cur. month | ihr_cur+1),gswm_sdi_vn(lon0:lon1,lat0:lat1,mon_cur,ihr_nxt+1), | difsec,secint,vn_curmo,1) call timeliper_array(gswm_sdi_z(lon0:lon1,lat0:lat1,mon_nxt, ! next month | ihr_cur+1),gswm_sdi_z(lon0:lon1,lat0:lat1,mon_nxt,ihr_nxt+1), | difsec,secint,z_nxtmo,1) call timeliper_array(gswm_sdi_tn(lon0:lon1,lat0:lat1,mon_nxt, ! next month | ihr_cur+1),gswm_sdi_tn(lon0:lon1,lat0:lat1,mon_nxt,ihr_nxt+1), | difsec,secint,tn_nxtmo,1) call timeliper_array(gswm_sdi_un(lon0:lon1,lat0:lat1,mon_nxt, ! next month | ihr_cur+1),gswm_sdi_un(lon0:lon1,lat0:lat1,mon_nxt,ihr_nxt+1), | difsec,secint,un_nxtmo,1) call timeliper_array(gswm_sdi_vn(lon0:lon1,lat0:lat1,mon_nxt, ! next month | ihr_cur+1),gswm_sdi_vn(lon0:lon1,lat0:lat1,mon_nxt,ihr_nxt+1), | difsec,secint,vn_nxtmo,1) endif ! ! Interpolate to modelday (iday) ! Check if interpolation 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 interpolation is necessary calculated the 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) if (nointp==0) then ! interpolation call timeliper_array(z_curmo,z_nxtmo,difday,dayint, | z_tmp,1) call timeliper_array(tn_curmo,tn_nxtmo,difday,dayint, | tn_tmp,1) call timeliper_array(un_curmo,un_nxtmo,difday,dayint, | un_tmp,1) call timeliper_array(vn_curmo,vn_nxtmo,difday,dayint, | vn_tmp,1) elseif(nointp == 1) then ! no interpolation same day as current month zsdi_gswm (lon0:lon1,lat0:lat1) = z_curmo tnsdi_gswm(lon0:lon1,lat0:lat1) = tn_curmo unsdi_gswm(lon0:lon1,lat0:lat1) = un_curmo vnsdi_gswm(lon0:lon1,lat0:lat1) = vn_curmo goto 30 elseif(nointp == 2) then ! no interpolation same day as next month zsdi_gswm (lon0:lon1,lat0:lat1) = z_nxtmo tnsdi_gswm(lon0:lon1,lat0:lat1) = tn_nxtmo unsdi_gswm(lon0:lon1,lat0:lat1) = un_nxtmo vnsdi_gswm(lon0:lon1,lat0:lat1) = vn_nxtmo goto 30 endif ! zsdi_gswm (lon0:lon1,lat0:lat1) = z_tmp tnsdi_gswm(lon0:lon1,lat0:lat1) = tn_tmp unsdi_gswm(lon0:lon1,lat0:lat1) = un_tmp vnsdi_gswm(lon0:lon1,lat0:lat1) = vn_tmp 30 continue ! ! do mpi periodic points exchange for gswm with f2d(:) ! Moved here from sub getgswm* because mpi necessary when f2d data ! is allocated only for task-local subdomain block. ! lons 1,2 <- nlonp4-3,nlonp4-2 and nlonp4-1,nlonp4 <- 3,4 #ifdef MPI call mp_periodic_f2d(zsdi_gswm(lon0:lon1,lat0:lat1), | lon0,lon1,lat0,lat1) call mp_periodic_f2d(tnsdi_gswm(lon0:lon1,lat0:lat1), | lon0,lon1,lat0,lat1) call mp_periodic_f2d(unsdi_gswm(lon0:lon1,lat0:lat1), | lon0,lon1,lat0,lat1) call mp_periodic_f2d(vnsdi_gswm(lon0:lon1,lat0:lat1), | lon0,lon1,lat0,lat1) # else set_periodic_f2d(zsdi_gswm) set_periodic_f2d(tnsdi_gswm) set_periodic_f2d(unsdi_gswm) set_periodic_f2d(vnsdi_gswm) #endif ! if (iprint > 0) then write(6,"('GSWM semidiurnal tidal pertubation interpolated ', | 'to date and time')") endif end subroutine getgswmsdi !----------------------------------------------------------------------- subroutine getgswmnmidi(iyear,iday,iutsec,iprint) ! ! Use gswm data read from gswm_nmidi_ncfile file to return perturbation in Z ! T,UN,VN at current date and time. It is assumed that the gswm_nmidi_ncfile ! 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 ! ! Args: integer,intent(in) :: iyear,iday,iutsec,iprint ! ! Local: integer :: i,j,ndmon_nl(13),ndmon_lp(13),nonleap, | mon_cur,mon_nxt,ihr_cur,ihr_nxt,nointp=0, | lonbeg,lonend real :: difsec,difday,sec,secint,dayint, | z_curmo(lon0:lon1,lat0:lat1), | z_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation | tn_curmo(lon0:lon1,lat0:lat1), | tn_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation | un_curmo(lon0:lon1,lat0:lat1), | un_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation | vn_curmo(lon0:lon1,lat0:lat1), | vn_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation | z_tmp(lon0:lon1,lat0:lat1), | tn_tmp(lon0:lon1,lat0:lat1), ! tmp | un_tmp(lon0:lon1,lat0:lat1), | vn_tmp(lon0:lon1,lat0:lat1) ! tmp ! ! 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 ! if (iprint > 0) then write(6,"(/,72('-'))") write(6,"('GETGSWMNMIDI: get nonmigrating diurnal tidal ', | 'perturbation from database:')") write(6,"('Initial requested iyear=',i4,' iday=',i3,' iutsec=', | i10)") iyear,iday,iutsec endif ! ! 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 if(ihr_nxt == 24 ) ihr_nxt = 1 ! ! GSWM hourly data from 0 UT (i=1) to 23 UT (i=24) resp. 0 utsec ! ! Check times: ! if (modeltime(3).lt.0.or.modeltime(3).gt.59) then write(6,"('>>> getgswmnmidi: bad imin=',3i4)") modeltime(3) stop 'getgswmnmidi' endif if (modeltime(4).lt.0.or.modeltime(4).gt.59) then write(6,"('>>> getgswmnmidi: bad isec=',3i4)") modeltime(4) stop 'getgswmnmidi' endif ! Interpolate to modeltime iutsec for each month: mon_cur mon_nxt if (iutsec==0) then ! no interpolation z_curmo(:,:) = gswm_nmidi_z (lon0:lon1,lat0:lat1, ! current month | mon_cur,ihr_cur+1) tn_curmo(:,:) = gswm_nmidi_tn(lon0:lon1,lat0:lat1, | mon_cur,ihr_cur+1) un_curmo(:,:) = gswm_nmidi_un(lon0:lon1,lat0:lat1, | mon_cur,ihr_cur+1) vn_curmo(:,:) = gswm_nmidi_vn(lon0:lon1,lat0:lat1, | mon_cur,ihr_cur+1) z_nxtmo(:,:) = gswm_nmidi_z (lon0:lon1,lat0:lat1, ! next month | mon_nxt,ihr_cur+1) tn_nxtmo(:,:) = gswm_nmidi_tn(lon0:lon1,lat0:lat1, | mon_nxt,ihr_cur+1) un_nxtmo(:,:) = gswm_nmidi_un(lon0:lon1,lat0:lat1, | mon_nxt,ihr_cur+1) vn_nxtmo(:,:) = gswm_nmidi_vn(lon0:lon1,lat0:lat1, | mon_nxt,ihr_cur+1) else difsec = 60.*60. ! difference in sec between ihr_cur and ihr_nxt secint = float(modeltime(3)*60 + modeltime(4)) ! interpolation time [sec] call timeliper_array(gswm_nmidi_z(lon0:lon1,lat0:lat1, ! cur. month | mon_cur,ihr_cur+1),gswm_nmidi_z(lon0:lon1,lat0:lat1,mon_cur, | ihr_nxt+1),difsec,secint,z_curmo,1) call timeliper_array(gswm_nmidi_tn(lon0:lon1,lat0:lat1, ! cur. month | mon_cur,ihr_cur+1),gswm_nmidi_tn(lon0:lon1,lat0:lat1,mon_cur, | ihr_nxt+1),difsec,secint,tn_curmo,1) call timeliper_array(gswm_nmidi_un(lon0:lon1,lat0:lat1, ! cur. month | mon_cur,ihr_cur+1),gswm_nmidi_un(lon0:lon1,lat0:lat1,mon_cur, | ihr_nxt+1),difsec,secint,un_curmo,1) call timeliper_array(gswm_nmidi_vn(lon0:lon1,lat0:lat1, ! cur. month | mon_cur,ihr_cur+1),gswm_nmidi_vn(lon0:lon1,lat0:lat1,mon_cur, | ihr_nxt+1),difsec,secint,vn_curmo,1) call timeliper_array(gswm_nmidi_z(lon0:lon1,lat0:lat1, ! next month | mon_nxt,ihr_cur+1),gswm_nmidi_z(lon0:lon1,lat0:lat1,mon_nxt, | ihr_nxt+1),difsec,secint,z_nxtmo,1) call timeliper_array(gswm_nmidi_tn(lon0:lon1,lat0:lat1, ! next month | mon_nxt,ihr_cur+1),gswm_nmidi_tn(lon0:lon1,lat0:lat1,mon_nxt, | ihr_nxt+1),difsec,secint,tn_nxtmo,1) call timeliper_array(gswm_nmidi_un(lon0:lon1,lat0:lat1, ! next month | mon_nxt,ihr_cur+1),gswm_nmidi_un(lon0:lon1,lat0:lat1,mon_nxt, | ihr_nxt+1),difsec,secint,un_nxtmo,1) call timeliper_array(gswm_nmidi_vn(lon0:lon1,lat0:lat1, ! next month | mon_nxt,ihr_cur+1),gswm_nmidi_vn(lon0:lon1,lat0:lat1,mon_nxt, | ihr_nxt+1),difsec,secint,vn_nxtmo,1) endif ! ! Interpolate to modelday (iday) ! ! Check if interpolation 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 interpolation is necessary calculated the 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) if (nointp==0) then ! interpolation call timeliper_array(z_curmo,z_nxtmo,difday,dayint, | z_tmp,1) call timeliper_array(tn_curmo,tn_nxtmo,difday,dayint, | tn_tmp,1) call timeliper_array(un_curmo,un_nxtmo,difday,dayint, | un_tmp,1) call timeliper_array(vn_curmo,vn_nxtmo,difday,dayint, | vn_tmp,1) elseif(nointp == 1) then ! no interpolation same day as current month znmidi_gswm (lon0:lon1,lat0:lat1) = z_curmo tnnmidi_gswm(lon0:lon1,lat0:lat1) = tn_curmo unnmidi_gswm(lon0:lon1,lat0:lat1) = un_curmo vnnmidi_gswm(lon0:lon1,lat0:lat1) = vn_curmo goto 30 elseif(nointp == 2) then ! no interpolation same day as next month znmidi_gswm (lon0:lon1,lat0:lat1) = z_nxtmo tnnmidi_gswm(lon0:lon1,lat0:lat1) = tn_nxtmo unnmidi_gswm(lon0:lon1,lat0:lat1) = un_nxtmo vnnmidi_gswm(lon0:lon1,lat0:lat1) = vn_nxtmo goto 30 endif ! znmidi_gswm (lon0:lon1,lat0:lat1) = z_tmp tnnmidi_gswm(lon0:lon1,lat0:lat1) = tn_tmp unnmidi_gswm(lon0:lon1,lat0:lat1) = un_tmp vnnmidi_gswm(lon0:lon1,lat0:lat1) = vn_tmp 30 continue ! ! do mpi periodic points exchange for gswm with f2d(:) ! Moved here from sub getgswm* because mpi necessary when f2d data ! is allocated only for task-local subdomain block. ! lons 1,2 <- nlonp4-3,nlonp4-2 and nlonp4-1,nlonp4 <- 3,4 #ifdef MPI call mp_periodic_f2d(znmidi_gswm(lon0:lon1,lat0:lat1), | lon0,lon1,lat0,lat1) call mp_periodic_f2d(tnnmidi_gswm(lon0:lon1,lat0:lat1), | lon0,lon1,lat0,lat1) call mp_periodic_f2d(unnmidi_gswm(lon0:lon1,lat0:lat1), | lon0,lon1,lat0,lat1) call mp_periodic_f2d(vnnmidi_gswm(lon0:lon1,lat0:lat1), | lon0,lon1,lat0,lat1) # else set_periodic_f2d(znmidi_gswm) set_periodic_f2d(tnnmidi_gswm) set_periodic_f2d(unnmidi_gswm) set_periodic_f2d(vnnmidi_gswm) #endif ! if (iprint > 0) then write(6,"('GSWM nonmigrating diurnal tidal pertubation ', | 'interpolated to date and time')") endif end subroutine getgswmnmidi !----------------------------------------------------------------------- subroutine getgswmnmisdi(iyear,iday,iutsec,iprint) ! ! Use gswm data read from gswm_nmisdi_ncfile file to return perturbation in Z ! T,UN,VN at current date and time. It is assumed that the gswm_nmisdi_ncfile ! 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 ! ! Args: integer,intent(in) :: iyear,iday,iutsec,iprint ! ! Local: integer :: i,j,ndmon_nl(13),ndmon_lp(13),nonleap, | mon_cur,mon_nxt,ihr_cur,ihr_nxt,nointp=0, | lonbeg,lonend real :: difsec,difday,sec,secint,dayint, | z_curmo(lon0:lon1,lat0:lat1), | z_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation | tn_curmo(lon0:lon1,lat0:lat1), | tn_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation | un_curmo(lon0:lon1,lat0:lat1), | un_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation | vn_curmo(lon0:lon1,lat0:lat1), | vn_nxtmo(lon0:lon1,lat0:lat1), ! UT interpolation | z_tmp(lon0:lon1,lat0:lat1), | tn_tmp(lon0:lon1,lat0:lat1), ! tmp | un_tmp(lon0:lon1,lat0:lat1), | vn_tmp(lon0:lon1,lat0:lat1) ! tmp ! ! 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 ! if (iprint > 0) then write(6,"(/,72('-'))") write(6,"('GETGSWMNMISDI: get nonmigrating semidiurnal tidal ', | 'perturbation from database:')") write(6,"('Initial requested iyear=',i4,' iday=',i3,' iutsec=', | i10)") iyear,iday,iutsec endif ! ! 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 if(ihr_nxt == 24 ) ihr_nxt = 1 ! ! GSWM hourly data from 0 UT (i=1) to 23 UT (i=24) resp. 0 utsec ! ! Check times: ! if (modeltime(3).lt.0.or.modeltime(3).gt.59) then write(6,"('>>> getgswmnmisdi: bad imin=',3i4)") modeltime(3) stop 'getgswmnmisdi' endif if (modeltime(4).lt.0.or.modeltime(4).gt.59) then write(6,"('>>> getgswmnmisdi: bad isec=',3i4)") modeltime(4) stop 'getgswmnmisdi' endif ! Interpolate to modeltime iutsec for each month: mon_cur mon_nxt if (iutsec==0) then ! no interpolation z_curmo(:,:) = gswm_nmisdi_z (lon0:lon1,lat0:lat1, ! current month | mon_cur,ihr_cur+1) tn_curmo(:,:) = gswm_nmisdi_tn(lon0:lon1,lat0:lat1, | mon_cur,ihr_cur+1) un_curmo(:,:) = gswm_nmisdi_un(lon0:lon1,lat0:lat1, | mon_cur,ihr_cur+1) vn_curmo(:,:) = gswm_nmisdi_vn(lon0:lon1,lat0:lat1, | mon_cur,ihr_cur+1) z_nxtmo(:,:) = gswm_nmisdi_z (lon0:lon1,lat0:lat1, ! next month | mon_nxt,ihr_cur+1) tn_nxtmo(:,:) = gswm_nmisdi_tn(lon0:lon1,lat0:lat1, | mon_nxt,ihr_cur+1) un_nxtmo(:,:) = gswm_nmisdi_un(lon0:lon1,lat0:lat1, | mon_nxt,ihr_cur+1) vn_nxtmo(:,:) = gswm_nmisdi_vn(lon0:lon1,lat0:lat1, | mon_nxt,ihr_cur+1) else difsec = 60.*60. ! difference in sec between ihr_cur and ihr_nxt secint = float(modeltime(3)*60 + modeltime(4)) ! interpolation time [sec] call timeliper_array(gswm_nmisdi_z(lon0:lon1,lat0:lat1, ! cur. month | mon_cur,ihr_cur+1),gswm_nmisdi_z(lon0:lon1,lat0:lat1,mon_cur, | ihr_nxt+1),difsec,secint,z_curmo,1) call timeliper_array(gswm_nmisdi_tn(lon0:lon1,lat0:lat1, ! cur. month | mon_cur,ihr_cur+1),gswm_nmisdi_tn(lon0:lon1,lat0:lat1,mon_cur, | ihr_nxt+1),difsec,secint,tn_curmo,1) call timeliper_array(gswm_nmisdi_un(lon0:lon1,lat0:lat1, ! cur. month | mon_cur,ihr_cur+1),gswm_nmisdi_un(lon0:lon1,lat0:lat1,mon_cur, | ihr_nxt+1),difsec,secint,un_curmo,1) call timeliper_array(gswm_nmisdi_vn(lon0:lon1,lat0:lat1, ! cur. month | mon_cur,ihr_cur+1),gswm_nmisdi_vn(lon0:lon1,lat0:lat1,mon_cur, | ihr_nxt+1),difsec,secint,vn_curmo,1) call timeliper_array(gswm_nmisdi_z(lon0:lon1,lat0:lat1, ! next month | mon_nxt,ihr_cur+1),gswm_nmisdi_z(lon0:lon1,lat0:lat1,mon_nxt, | ihr_nxt+1),difsec,secint,z_nxtmo,1) call timeliper_array(gswm_nmisdi_tn(lon0:lon1,lat0:lat1, ! next month | mon_nxt,ihr_cur+1),gswm_nmisdi_tn(lon0:lon1,lat0:lat1,mon_nxt, | ihr_nxt+1),difsec,secint,tn_nxtmo,1) call timeliper_array(gswm_nmisdi_un(lon0:lon1,lat0:lat1, ! next month | mon_nxt,ihr_cur+1),gswm_nmisdi_un(lon0:lon1,lat0:lat1,mon_nxt, | ihr_nxt+1),difsec,secint,un_nxtmo,1) call timeliper_array(gswm_nmisdi_vn(lon0:lon1,lat0:lat1, ! next month | mon_nxt,ihr_cur+1),gswm_nmisdi_vn(lon0:lon1,lat0:lat1,mon_nxt, | ihr_nxt+1),difsec,secint,vn_nxtmo,1) endif ! ! Interpolate to modelday (iday) ! Check if interpolation 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 interpolation is necessary calculated the 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) if (nointp==0) then ! interpolation call timeliper_array(z_curmo,z_nxtmo,difday,dayint, | z_tmp,1) call timeliper_array(tn_curmo,tn_nxtmo,difday,dayint, | tn_tmp,1) call timeliper_array(un_curmo,un_nxtmo,difday,dayint, | un_tmp,1) call timeliper_array(vn_curmo,vn_nxtmo,difday,dayint, | vn_tmp,1) elseif(nointp == 1) then ! no interpolation same day as current month znmisdi_gswm (lon0:lon1,lat0:lat1) = z_curmo tnnmisdi_gswm(lon0:lon1,lat0:lat1) = tn_curmo unnmisdi_gswm(lon0:lon1,lat0:lat1) = un_curmo vnnmisdi_gswm(lon0:lon1,lat0:lat1) = vn_curmo goto 30 elseif(nointp == 2) then ! no interpolation same day as next month znmisdi_gswm (lon0:lon1,lat0:lat1) = z_nxtmo tnnmisdi_gswm(lon0:lon1,lat0:lat1) = tn_nxtmo unnmisdi_gswm(lon0:lon1,lat0:lat1) = un_nxtmo vnnmisdi_gswm(lon0:lon1,lat0:lat1) = vn_nxtmo goto 30 endif ! znmisdi_gswm (lon0:lon1,lat0:lat1) = z_tmp tnnmisdi_gswm(lon0:lon1,lat0:lat1) = tn_tmp unnmisdi_gswm(lon0:lon1,lat0:lat1) = un_tmp vnnmisdi_gswm(lon0:lon1,lat0:lat1) = vn_tmp 30 continue ! ! do mpi periodic points exchange for gswm with f2d(:) ! Moved here from sub getgswm* because mpi necessary when f2d data ! is allocated only for task-local subdomain block. ! lons 1,2 <- nlonp4-3,nlonp4-2 and nlonp4-1,nlonp4 <- 3,4 #ifdef MPI call mp_periodic_f2d(znmisdi_gswm(lon0:lon1,lat0:lat1), | lon0,lon1,lat0,lat1) call mp_periodic_f2d(tnnmisdi_gswm(lon0:lon1,lat0:lat1), | lon0,lon1,lat0,lat1) call mp_periodic_f2d(unnmisdi_gswm(lon0:lon1,lat0:lat1), | lon0,lon1,lat0,lat1) call mp_periodic_f2d(vnnmisdi_gswm(lon0:lon1,lat0:lat1), | lon0,lon1,lat0,lat1) # else set_periodic_f2d(znmisdi_gswm) set_periodic_f2d(tnnmisdi_gswm) set_periodic_f2d(unnmisdi_gswm) set_periodic_f2d(vnnmisdi_gswm) #endif ! if (iprint > 0) then write(6,"('GSWM nonmigrating semidiurnal tidal pertubation ', | 'interpolated to date and time')") endif end subroutine getgswmnmisdi !----------------------------------------------------------------------- subroutine timeliper_array(d1,d2,difd1d2,difint,fout, | iprnt) c c Interpolate from fields d1,d2 linearly to time difint c (d1 must be at 0 unit time), returning fout c On input: c d1,d2 = the field (d1 at 0 unit, d2 at difd1d2 units) c difd1d2 = time difference [unit] between d1 & d2 c difint = time of interpolation [unit] (counted from d1=0 time) c In output: c fout is defined at difint c ! 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 timeliper_array !----------------------------------------------------------------------- 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