program mgcmlbc_rotate ! ! 3/19/07 btf: ! Read an MGCM_LBC_FILE nc file, rotate the fields 12 hours, and ! write equivalent new file for use by mtgcm1.0 and later versions. ! implicit none #include "netcdf.inc" character(len=240) :: infile,outfile integer :: n,ncid_in,ncid_out,istat ! infile = ' ' outfile = ' ' ! infile = 'TGCM.nc' ! outfile= 'MGCM.nc' ! --------------------------------------------------------------------- ! JPL/CDP Studies: Fall 2010 - Spring 2012 (Rotation by 12-hours) ! --------------------------------------------------------------------- ! infile = 'TGCM_lbnd_comp_Ls270_TAU0.3-GW-VAR200-KAP2.5E-6.nc' ! outfile= 'MGCM_lbnd_comp_Ls270_TAU0.3-GW-VAR200-KAP2.5E-6.nc' ! infile = 'TGCM_lbnd_comp_Ls270_TAU0.3-noGW-noRF.nc' ! outfile= 'MGCM_lbnd_comp_Ls270_TAU0.3-noGW-noRF.nc' infile = 'TGCM_lbnd-Ls270-GW200-2.5E-6_TESyr1-CONR0.005-noswap.nc' outfile= 'MGCM_lbnd-Ls270-GW200-2.5E-6_TESyr1-CONR0.005-noswap.nc' istat = nf_open(infile,NF_NOWRITE,ncid_in) istat = nf_create(outfile,NF_CLOBBER,ncid_out) call readfile(ncid_in,ncid_out) end ! program mgcmlbc_rotate !----------------------------------------------------------------------- module globals save integer :: nlon,nlat,ntime real,allocatable,dimension(:) :: glon,glat,time end module globals !----------------------------------------------------------------------- subroutine readfile(ncid,ncid_out) use globals,only: nlon,nlat,ntime,glon,glat,time implicit none #include "netcdf.inc" integer,intent(in) :: ncid,ncid_out integer :: istat,icount3(3),istart3(3) integer :: id_time,id_lon,id_lat,id_t,id_u,id_v,id_z integer :: idv_time,idv_lon,idv_lat real,allocatable,dimension(:,:,:) :: t,u,v,z ! (nlon,nlat,ntime) istat = nf_inq_dimid(ncid,'time',id_time) istat = nf_inq_dimlen(ncid,id_time,ntime) istat = nf_inq_dimid(ncid,'lon',id_lon) istat = nf_inq_dimlen(ncid,id_lon,nlon) istat = nf_inq_dimid(ncid,'lat',id_lat) istat = nf_inq_dimlen(ncid,id_lat,nlat) istat = nf_inq_varid(ncid,'Z',id_z) istat = nf_inq_varid(ncid,'TN',id_t) istat = nf_inq_varid(ncid,'UN',id_u) istat = nf_inq_varid(ncid,'VN',id_v) allocate(time(ntime),stat=istat) istat = nf_inq_varid(ncid,'TIME',idv_time) ! get id of time var istat = nf_get_var_double(ncid,idv_time,time) ! read time var write(6,"('Read time var: time min,max=',2e12.4)") | minval(time),maxval(time) allocate(glon(nlon),stat=istat) istat = nf_inq_varid(ncid,'LON',idv_lon) ! get id of lon var istat = nf_get_var_double(ncid,idv_lon,glon) ! read lon var ! write(6,"('Read lon var: lon=',/,(6e12.4))") glon allocate(glat(nlat),stat=istat) istat = nf_inq_varid(ncid,'LAT',idv_lat) ! get id of lat var istat = nf_get_var_double(ncid,idv_lat,glat) ! read lat var ! write(6,"('Read lat var: lat=',/,(6e12.4))") glat ! write(6,"('readfile allocating t,u,v, nlon=',i4,' nlat=', ! | i4,' ntime=',i4,' product=',i8)") ! | nlon,nlat,ntime,nlon*nlat*ntime allocate(t(nlon,nlat,ntime),stat=istat) allocate(u(nlon,nlat,ntime),stat=istat) allocate(v(nlon,nlat,ntime),stat=istat) allocate(z(nlon,nlat,ntime),stat=istat) istart3(1:3) = 1 icount3(1) = nlon icount3(2) = nlat icount3(3) = ntime istat = nf_get_vara_double(ncid,id_t,istart3,icount3,t) write(6,"('Read t: min,max=',2e12.4)") minval(t),maxval(t) istat = nf_get_vara_double(ncid,id_u,istart3,icount3,u) write(6,"('Read u: min,max=',2e12.4)") minval(u),maxval(u) istat = nf_get_vara_double(ncid,id_v,istart3,icount3,v) write(6,"('Read v: min,max=',2e12.4)") minval(v),maxval(v) istat = nf_get_vara_double(ncid,id_z,istart3,icount3,z) write(6,"('Read z: min,max=',2e12.4)") minval(z),maxval(z) call rotate_fields(ncid_out,t,u,v,z,nlon,nlat,ntime) end subroutine readfile !----------------------------------------------------------------------- subroutine rotate_fields(ncid_out,t,u,v,z,nlon,nlat,ntime) implicit none ! ! Args: integer,intent(in) :: ncid_out,nlon,nlat,ntime real,intent(in),dimension(nlon,nlat,ntime) :: t,u,v,z ! ! Local: real,dimension(nlon,nlat,ntime) :: tm,um,vm,zm real,parameter :: secs_per_day = 24.*3600. real :: deg_per_sec,deg_rot,deg_per_lon integer :: nlonrot,itime,i,j,ii real,dimension(nlon) :: fz,ft,fu,fv ! deg_per_sec = 360./secs_per_day deg_rot = secs_per_day/2. * deg_per_sec deg_per_lon = 360./real(nlon) nlonrot = nint(deg_rot/deg_per_lon) ! write(6,"('rotate_fields: deg_per_sec=',f12.4,' deg_rot=',f12.4, ! | ' deg_per_lon=',f12.4,' nlonrot=',i8)") ! | deg_per_sec,deg_rot,deg_per_lon,nlonrot do itime=1,ntime do j=1,nlat fz(:) = z(:,j,itime) ft(:) = t(:,j,itime) fu(:) = u(:,j,itime) fv(:) = v(:,j,itime) do i=1,nlon ii = mod(i+nlonrot-1,nlon)+1 zm(i,j,itime) = fz(ii) tm(i,j,itime) = ft(ii) um(i,j,itime) = fu(ii) vm(i,j,itime) = fv(ii) enddo enddo enddo write(6,"('rotate_fields: zm min,max=',2e12.4)") | minval(zm),maxval(zm) write(6,"('rotate_fields: tm min,max=',2e12.4)") | minval(tm),maxval(tm) write(6,"('rotate_fields: um min,max=',2e12.4)") | minval(um),maxval(um) write(6,"('rotate_fields: vm min,max=',2e12.4)") | minval(vm),maxval(vm) call writefile(ncid_out,tm,um,vm,zm,nlon,nlat,ntime) end subroutine rotate_fields !----------------------------------------------------------------------- subroutine writefile(ncid,tm,um,vm,zm,nlon,nlat,ntime) use globals,only: glon,glat,time implicit none #include "netcdf.inc" ! ! Args: integer,intent(in) :: ncid,nlon,nlat,ntime real,intent(in),dimension(nlon,nlat,ntime) :: tm,um,vm,zm ! ! Local: integer :: i,istat,ids1(1),ids3(3),start_3d(3),count_3d(3) integer :: id_lon,id_lat,id_time integer :: idv_lon,idv_lat,idv_time integer :: idv_t,idv_u,idv_v,idv_z character(len=80) :: char80 ! ! Define time and grid dimensions: istat = nf_def_dim(ncid,"time",NF_UNLIMITED,id_time) istat = nf_def_dim(ncid,"lon",nlon,id_lon) istat = nf_def_dim(ncid,"lat",nlat,id_lat) ids1(1) = id_lon istat = nf_def_var(ncid,"lon",NF_DOUBLE,1,ids1,idv_lon) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining longitude dimension variable') write(char80,"('Geographic Longitude (-west, +east)')") istat = nf_put_att_text(ncid,idv_lon,"long_name", | len_trim(char80),trim(char80)) istat = nf_put_att_text(ncid,idv_lon,"units",12,'degrees_east') ids1(1) = id_lat istat = nf_def_var(ncid,"lat",NF_DOUBLE,1,ids1,idv_lat) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining latitude dimension variable') write(char80,"('Geographic Latitude (-south, +north)')") istat = nf_put_att_text(ncid,idv_lat,"long_name", | len_trim(char80),trim(char80)) istat = nf_put_att_text(ncid,idv_lat,"units",13,'degrees_north') ids1(1) = id_time istat = nf_def_var(ncid,"TIME",NF_DOUBLE,1,ids1,idv_time) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining time dimension variable') istat = nf_put_att_text(ncid,idv_time,"long_name",15, | "MGCM model time") istat = nf_put_att_text(ncid,idv_time,"units",4,"SOLS") write(6,"('idv_time=',i3)") idv_time ids3(1) = id_lon ids3(2) = id_lat ids3(3) = id_time istat = nf_def_var(ncid,"Z",NF_DOUBLE,3,ids3,idv_z) if (istat /= NF_NOERR) call handle_ncerr(istat,'Error defining Z') write(char80,"('GEOPOTENTIAL HEIGHT')") istat = nf_put_att_text(ncid,idv_z,"long_name",len_trim(char80), | trim(char80)) istat = nf_put_att_text(ncid,idv_z,"units",2,'KM') ! istat = nf_def_var(ncid,"TN",NF_DOUBLE,3,ids3,idv_t) if (istat /= NF_NOERR) call handle_ncerr(istat,'Error defining T') write(char80,"('NEUTRAL TEMPERATURE')") istat = nf_put_att_text(ncid,idv_t,"long_name",len_trim(char80), | trim(char80)) istat = nf_put_att_text(ncid,idv_t,"units",5,'DEG K') ! istat = nf_def_var(ncid,"UN",NF_DOUBLE,3,ids3,idv_u) if (istat /= NF_NOERR) call handle_ncerr(istat,'Error defining U') write(char80,"('NEUTRAL ZONAL WIND (-WEST, +EAST)')") istat = nf_put_att_text(ncid,idv_u,"long_name",len_trim(char80), | trim(char80)) istat = nf_put_att_text(ncid,idv_u,"units",4,'CM/S') ! istat = nf_def_var(ncid,"VN",NF_DOUBLE,3,ids3,idv_v) if (istat /= NF_NOERR) call handle_ncerr(istat,'Error defining V') write(char80,"('NEUTRAL MERIDIONAL WIND (-SOUTH, +NORTH)')") istat = nf_put_att_text(ncid,idv_v,"long_name",len_trim(char80), | trim(char80)) istat = nf_put_att_text(ncid,idv_v,"units",4,'CM/S') ! ! Exit define mode: istat = nf_enddef(ncid) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error return from nf_enddef') ! ! Write time: write(6,"('writing time: idv_time=',i4,' min,max=',2e12.4)") | idv_time,minval(time),maxval(time) do i=1,ntime istat = nf_put_var1_double(ncid,idv_time,i,time(i)) enddo ! istat = nf_put_var_double(ncid,idv_time,time) ! if (istat /= NF_NOERR) then ! call handle_ncerr(istat,'Error writing time') ! stop 'time' ! endif ! ! Write glat,glon: istat = nf_put_var_double(ncid,idv_lon,glon) if (istat /= NF_NOERR) | call handle_ncerr(istat,'Error writing glon') istat = nf_put_var_double(ncid,idv_lat,glat) if (istat /= NF_NOERR) | call handle_ncerr(istat,'Error writing glat') start_3d(1:3) = 1 count_3d(1) = nlon count_3d(2) = nlat count_3d(3) = ntime ! write(6,"('Writing tm..')") istat = nf_put_vara_double(ncid,idv_t,start_3d,count_3d,tm) if (istat /= NF_NOERR) | call handle_ncerr(istat,'Error writing tm') ! write(6,"('Writing um..')") istat = nf_put_vara_double(ncid,idv_u,start_3d,count_3d,um) if (istat /= NF_NOERR) | call handle_ncerr(istat,'Error writing um') ! write(6,"('Writing vm..')") istat = nf_put_vara_double(ncid,idv_v,start_3d,count_3d,vm) if (istat /= NF_NOERR) | call handle_ncerr(istat,'Error writing vm') ! write(6,"('Writing zm..')") istat = nf_put_vara_double(ncid,idv_z,start_3d,count_3d,zm) if (istat /= NF_NOERR) | call handle_ncerr(istat,'Error writing zm') istat = nf_close(ncid) end subroutine writefile !----------------------------------------------------------------------- subroutine handle_ncerr(istat,msg) implicit none #include "netcdf.inc" ! ! Handle a netcdf lib error: ! integer,intent(in) :: istat character(len=*),intent(in) :: msg ! write(6,"(/72('-'))") write(6,"('>>> Error from netcdf library:')") write(6,"(a)") trim(msg) write(6,"('istat=',i5)") istat write(6,"(a)") nf_strerror(istat) write(6,"(72('-')/)") return end subroutine handle_ncerr