! module nchist_module ! ! Module to create, read and write netcdf history files: ! use params_module,only: nlon,nlev,nmlev,nmlevp1,nlat,nlevp1,spval, | nlonp4,nlonp1,nmlonp1,nmlat,nmlon,nmagphrlon,nmagphrlat,mxfsech, | nlonp2,zmbot,zmtop,zibot,zitop use fields_module,only: nf4d,nf4d_hist,f4d,fakeflds,shortname_len, | itc,itp,fields_4d,fields_3d,fields_2d,poten,foutput, | tlbc,ulbc,vlbc,tlbc_glb,ulbc_glb,vlbc_glb, | tlbc_nm,ulbc_nm,vlbc_nm,tlbc_nm_glb,ulbc_nm_glb,vlbc_nm_glb, | nf3d,f3d,fzg use input_module,only: mxlen_filename, | gswm_mi_di_ncfile, gswm_mi_sdi_ncfile, | gswm_nm_di_ncfile, gswm_nm_sdi_ncfile, | wacxg5_nfiles, wacxg5_nlev use init_module,only: glon,glat,zpmid,zpint,gmlon,gmlat,istep, | start_mtime,zpimag,zpmag implicit none ! ! Included file svn_version.inc contains the svn revision number used in the build. ! (character(len=8) svn_version) ! #include "svn_version.inc" ! #ifdef SUN #include #else #include #endif ! integer,private :: | id_time, ! dimension id for time (unlimited) | id_mtimedim, ! dimension id for model time (3) | id_lon,id_lat, ! horizontal dimension ids | id_lev,id_ilev, ! vertical dimension ids (midpoints,interfaces) | id_mlon,id_mlat, ! horizontal dimension ids for magnetic grid | id_mlev,id_imlev, ! vertical dimension ids for magnetic grid | id_magphrlon,id_magphrlat, ! dimension ids for 2d magnetosphere grid | id_latlon, ! dimension id for coords (2) | id_dtide, ! dimension id for dtide (2) | id_sdtide, ! dimension id for sdtide (10) | id_fake, ! fake dimension id (1) | id_datelen, ! length of date and time | id_filelen, ! max length of a file name | ids1(1),ids2(2),ids3(3),ids4(4), ! vectors of dim id's | start_2d(2),count_2d(2), ! start loc and count for 2d vars | start_3d(3),count_3d(3), ! start loc and count for 3d vars | start_4d(4),count_4d(4), ! start loc and count for 4d vars ! ! Variable id's (idv_xxxx). Note lon, lat, lev, and time are ! coordinate variables (lon(lon), lat(lat), lev(lev), time(time)). ! Variables starting at mtime are "primary" variables. ! idv_f4d(nf4d) are id's for prognostic 4d fields. ! | idv_time,idv_mtime,idv_ut,idv_calendar_adv, ! variable id's | idv_lon,idv_lat,idv_lev,idv_ilev, | idv_lbc,idv_tlbc,idv_ulbc,idv_vlbc, | idv_tlbc_nm,idv_ulbc_nm,idv_vlbc_nm, | idv_mlon,idv_mlat,idv_mlev,idv_imlev, | idv_magphrlon,idv_magphrlat,idv_grav, | idv_step,idv_iter,idv_year,idv_day,idv_p0,idv_p0_model, | idv_hpower,idv_ctpoten,idv_byimf,idv_colfac,idv_ncep_ncfile, | idv_ncep_reanalysis,idv_ecmwf_ncfile, | idv_gswm_mi_di_ncfile,idv_gswm_mi_sdi_ncfile, | idv_gswm_nm_di_ncfile,idv_gswm_nm_sdi_ncfile, | idv_f107d,idv_f107a,idv_gpi_ncfile, | idv_mag,idv_dtide,idv_sdtide,idv_nflds,idv_f4d(nf4d), | idv_alfa30,idv_e30,idv_ed2,idv_alfad2, | idv_writedate,idv_wacxg5_ncfile integer :: idv_fsech(mxfsech) ! ! If check_nan is set, check for presence of NaN's in data when writing ! history. Sub check_nans (util.F) will be called from subs wrf4d, wrf3d, ! wrf2d, and rdf4d. This sub works only on IBM-AIX. ! #ifdef AIX logical :: check_nan = .true. #else logical :: check_nan = .false. #endif contains !----------------------------------------------------------------------- subroutine nc_open(ncid,flnm,status,rdwr) ! ! Open or create a netcdf dataset flnm, returning netcdf file ! id handle in ncid. ! If status=='NEW' or 'REPLACE', a new netcdf file is created, ! otherwise (status=='OLD'), open existing dataset). ! If an existing file is opened and rdwr=='READ', the dataset is ! opened read-only (NF_NOWRITE), otherwise (rdwr=='WRITE' or ! 'APPEND') the dataset is opened read/write (NF_WRITE). ! If the open fails, print error msg and return ncid == 0. ! ! Args: integer,intent(out) :: ncid ! netcdf file id character(len=*),intent(in) :: flnm,status,rdwr ! ! Local: integer :: istat character(len=120) :: char120 ! ncid = 0 if (status=='NEW'.or.status=='REPLACE') then ! ! Create new dataset (clobber old one if it pre-exists): ! istat = nf_create(flnm,NF_CLOBBER,ncid) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_create for netcdf ', | 'file ',a,' status=',a)") trim(flnm),status call handle_ncerr(istat,char120) ncid = 0 else write(6,"('Created netcdf file ',a)") trim(flnm) endif ! ! Open pre-existing dataset: ! NF_WRITE opens with read-write, where write means append, change, etc ! NF_NOWRITE is read-only. ! elseif (status=='OLD') then if (rdwr=='WRITE'.or.rdwr=='APPEND') then istat = nf_open(flnm,NF_WRITE,ncid) elseif (rdwr=='READ') then istat = nf_open(flnm,NF_NOWRITE,ncid) else write(6,"('>>> nc_open: unrecognized rdwr=',a)") trim(rdwr) write(6,"(' Will open NF_WRITE (read/write).')") istat = nf_open(flnm,NF_WRITE,ncid) endif if (istat /= NF_NOERR) then write(char120,"('Error return from nf_open for netcdf ', | 'file ',a,' status=',a)") trim(flnm),status call handle_ncerr(istat,char120) ncid = 0 else ! write(6,"('Opened existing netcdf file ',a,' ncid=',i8)") ! | trim(flnm),ncid endif ! ! Unknown status request: else write(6,"(/,'>>> nc_open: unrecognized status = ',a)") | status ncid = 0 endif end subroutine nc_open !------------------------------------------------------------------- subroutine nc_close(ncid) ! ! Close a netcdf dataset file: ! ! Arg: integer,intent(in) :: ncid ! ! Local: integer :: istat ! istat = nf_close(ncid) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error return from nf_close') end subroutine nc_close !----------------------------------------------------------------------- subroutine nc_rdhist(ncid,diskfile,mtime,itime,ier) ! ! Read history at model time mtime from netcdf file diskfile, ! currently open under ncid. Define output history in global ! history structure h (module hist_module in history.f). ! Prognostic fields are transferred to global fg-array (wrf4d). ! If there is a problem (e.g., mtime not found), return ier > 0 ! use hist_module,only: sh,h,hist_initype,hist_print,nsource, | nsecsource use input_module,only: start_day,start_year,start,output use params_module,only: ispval,mxhvols use mk_hvols,only: mkhvols ! ! Args: integer,intent(in) :: ncid,mtime(3) character(len=*) :: diskfile integer,intent(out) :: itime,ier ! ! Local: integer,parameter :: mxdims=20 ! assume <= mxdims dimensions integer :: istat,id_time,id_mtimedim,mtime_read(3),ntimes, | n,i,ix,istart2(2),icount2(2),j,irdf4d(nf4d_hist),nstarts integer :: ndims,nvars,ngatts,idunlim,itype,natts,len, | iddims(mxdims),iscalar,nhists_total,nfiles_prim,nout character(len=80) :: varname,varnametmp real :: scalar,var1(1),var22(2,2),var2(2),var10(10),varlbc real,dimension(nlevp1) :: zpmid_rd,zpint_rd character(len=1024) :: hvols(mxhvols) ! ! External: logical,external :: arrayeq ! util.F integer,external :: numfiles ! write(6,"('Reading source history from diskfile ',a,':')") | trim(diskfile) ier = 0 ! ! Get id for mtime variable: istat = nf_inq_varid(ncid,"mtime",idv_mtime) if (istat /= NF_NOERR) then call handle_ncerr(istat,'Error getting id of variable mtime') ier = 1 return endif ! ! Get number of histories (times) on the file: istat = nf_inq_unlimdim(ncid,id_time) ! id of unlimited record var if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error getting unlimited dimension id') istat = nf_inq_dimlen(ncid,id_time,ntimes) ! length of time var if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error getting length of time record variable') ! ! If source file was not provided, determine if enough output files were ! provided (inp_primhist in input.F could not do this): if (nsource == 0) then n = size(start)-count(start==ispval) nstarts = n/3 ! number of start times given nout = mkhvols(output,hvols,mxhvols) nfiles_prim=numfiles('prim',nstarts,nsource,ntimes,nhists_total) if (nout < nfiles_prim) then write(6,"(/,'>>> INPUT: Will need ',i3,' OUTPUT files, but', | ' read only ',i3)") nfiles_prim,nout call shutdown('OUTPUT') endif endif ! ! Search for requested model time. When found, retain time index itime: istart2(1) = 1 icount2(1) = 3 icount2(2) = 1 itime = 0 do i=1,ntimes istart2(2) = i istat = nf_get_vara_int(ncid,idv_mtime,istart2,icount2, | mtime_read) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error getting mtime values') write(6,"('nc_rdhist: seeking ',3i4,' found ',3i4,' n=',i3)") | mtime,mtime_read,i if (all(mtime_read==mtime)) then itime = i exit endif enddo if (itime==0) then write(6,"(/,'>>> nc_rdhist: mtime ',3i4,' NOT found on ', | 'diskfile ',a)") mtime,trim(diskfile) ier = 1 return endif ! ! Define history structure from netcdf vars: ! ! Initialize history structure: call hist_initype(sh,istep) ! initialize source history structure ! ! Define vars already found above: sh%ihist = itime sh%modeltime(1:3) = mtime sh%modeltime(4) = 0 ! ! Read global attributes: istat = nf_get_att_text(ncid,NF_GLOBAL,"label",sh%label) istat = nf_get_att_text(ncid,NF_GLOBAL,"create_date", | sh%create_date) istat = nf_get_att_text(ncid,NF_GLOBAL,"logname",sh%logname) istat = nf_get_att_text(ncid,NF_GLOBAL,"host",sh%host) istat = nf_get_att_text(ncid,NF_GLOBAL,"system",sh%system) istat = nf_get_att_text(ncid,NF_GLOBAL,"model_name",sh%model_name) istat = nf_get_att_text(ncid,NF_GLOBAL,"model_version", | sh%model_version) istat = nf_get_att_int(ncid,NF_GLOBAL,"source_mtime", | sh%source_mtime) istat = nf_get_att_text(ncid,NF_GLOBAL,"initial_file", | sh%initial_file) if (istat /= NF_NOERR) ! was not on "old" histories | sh%initial_file = ' ' ! ! History_type is either "primary" or "secondary": istat=nf_get_att_text(ncid,NF_GLOBAL,"history_type",sh%hist_type) ! ! Check for tuv_lbc_intop (lbc of t,u,v stored in top slot) ! tuv_lbc_intop == 0 --> new history (tuv lbc *not* stored in top slot) ! tuv_lbc_intop == 1 --> old history (tuv lbc *is* stored in top slot) ! If tuv_lbc_intop global attribute is not on the file, assume ! "old" history, i.e., tuv_lbc_intop==1: ! istat=nf_get_att_int(ncid,NF_GLOBAL,"tuv_lbc_intop", | sh%tuv_lbc_intop) if (istat /= NF_NOERR) then ! old hist write(6,"('nc_rdhist: did not find tuv_lbc_intop in source', | ' history ',a,' --> assuming tuv_lbc_intop==1')") | trim(diskfile) sh%tuv_lbc_intop = 1 endif ! ! History_type is either "primary" or "secondary": istat = nf_get_att_text(ncid,NF_GLOBAL,"history_type", | sh%hist_type) ! ! Get number of dims, vars, atts, and id of unlimited dimension: istat = nf_inq(ncid,ndims,nvars,ngatts,idunlim) ! ! Get dimensions and check against current model parameters: ! (i.e., current model run at resolution A should not read a ! source history from a model that was run with resolution B) ! do i=1,ndims istat = nf_inq_dim(ncid,i,varname,len) if (trim(varname)=='lon') then if (len /= nlon) then write(6,"(/,'>>> nc_rdhist: length of lon read from', | ' history does not match parameter nlon.')") write(6,"(' len=',i4,' nlon=',i4)") len,nlon call shutdown('nlon') endif sh%nlon = len endif if (trim(varname)=='lat') then if (len /= nlat) then write(6,"(/,'>>> nc_rdhist: length of lat read from', | ' history does not match parameter nlat.')") write(6,"(' len=',i4,' nlat=',i4)") len,nlat call shutdown('nlat') endif sh%nlat = len endif if (trim(varname)=='lev') then ! ! If number of levels read from source history (len) is less than the ! number of levels allocated in the model (nlev), this is fatal: ! if (len < nlevp1) then write(6,"(/,'>>> nc_rdhist: length of lev dimension', | ' of history read is less than model parameter nlevp1.')") write(6,"(' len=',i4,' nlevp1=',i4)") len,nlevp1 call shutdown('nlevp1') ! ! If number of levels read from source history (len) is greater than ! the number of levels allocated in the model (nlevp1), then we will ! read the first nlevp1 levels (e.g., if nlevp1==88 and len==96, then ! read levels 1-88 from the history i.e., the bottom 88 levels). ! elseif (len > nlevp1) then write(6,"(/,'NOTE nc_rdhist: number of levels in source', | ' history (len=',i3,') exceeds the ',/,2x, | 'number of levels allocated in the model (nlevp1=', | i3,')')") len,nlevp1 write(6,"(' Will read the first ',i3,' levels from the', | ' source into the history buffer.')") nlevp1 endif sh%nlev = len endif enddo ! ! Get variables (should not need to get coord vars lat,lon,lev) do i=1,nvars istat = nf_inq_var(ncid,i,varname,itype,ndims,iddims,natts) select case(trim(varname)) case('lon') istat = nf_get_var_double(ncid,i,glon) case('lat') istat = nf_get_var_double(ncid,i,glat) ! ! New hist: lev(lev) is midpoints case('lev') if (sh%tuv_lbc_intop <= 0) then ! new hist (lev is midpoints) istat = nf_get_var_double(ncid,i,zpmid_rd) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading zpmid (new history)') ! write(6,"('nc_rdhist sh%nlev=',i3,': read zpmid from new', ! | ' hist: zpmid_rd=',/,(10f7.3))") sh%nlev,zpmid_rd if (.not.arrayeq(zpmid,zpmid_rd,nlevp1)) then write(6,"(/,'>>> nc_rdhist: zpmid_rd (as read from new', | ' history) is not equal to model zpmid')") write(6,"('nlevp1=',i3,' zpmid_rd=',/,(10f7.3))") | nlevp1,zpmid_rd write(6,"('nlevp1=',i3,' zpmid =',/,(10f7.3))") | nlevp1,zpmid call shutdown('zpmid') endif ! ! Old hist lev(lev) is interfaces: else ! old hist (lev is interfaces) istat = nf_get_var_double(ncid,i,zpint_rd) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading zpint (old history)') ! write(6,"('nc_rdhist sh%nlev=',i3,': read zpint from old', ! | ' hist: zpint_rd=',/,(10f7.3))") sh%nlev,zpint_rd if (.not.arrayeq(zpint,zpint_rd,nlevp1)) then write(6,"(/,'>>> nc_rdhist: zpint_rd (as read from old', | ' history) is not equal to model zpint')") write(6,"('nlevp1=',i3,' zpint_rd=',/,(10f7.3))") | nlevp1,zpint_rd write(6,"('nlevp1=',i3,' zpint =',/,(10f7.3))") | nlevp1,zpint call shutdown('zpint') endif endif ! ! ilev(ilev) interface levels are on new histories, or histories written ! by mksrc.ncl (i.e., new lbc from old histories): case('ilev') if (sh%tuv_lbc_intop <= 0) then ! new hist (lev is midpoints) istat = nf_get_var_double(ncid,i,zpint_rd) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading zpint (new history)') if (.not.arrayeq(zpint,zpint_rd,nlevp1)) then write(6,"(/,'>>> nc_rdhist: zpint_rd (as read from new', | ' history) is not equal to model zpint')") write(6,"('nlevp1=',i3,' zpint_rd=',/,(10f7.3))") | nlevp1,zpint_rd write(6,"('nlevp1=',i3,' zpint =',/,(10f7.3))") | nlevp1,zpint call shutdown('zpint') endif else ! old hist ilev interfaces (probably written by mksrc.ncl) istat = nf_get_var_double(ncid,i,zpint_rd) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading zpint (old history)') if (.not.arrayeq(zpint,zpint_rd,nlevp1)) then write(6,"(/,'>>> nc_rdhist: zpint_rd (as read from new', | ' history) is not equal to model zpint')") write(6,"('nlevp1=',i3,' zpint_rd=',/,(10f7.3))") | nlevp1,zpint_rd write(6,"('nlevp1=',i3,' zpint =',/,(10f7.3))") | nlevp1,zpint call shutdown('zpint') endif endif ! ! Magnetic coordinates: case('mlon') istat = nf_get_var_double(ncid,i,gmlon) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading gmlon') case('mlat') istat = nf_get_var_double(ncid,i,gmlat) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading gmlat') case('mlev') ! mag midpoints istat = nf_get_var_double(ncid,i,zpmag) ! zpmag is in init module if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading zpmag') case('imlev') ! mag interfaces istat = nf_get_var_double(ncid,i,zpimag) ! zpimag is in init module if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading zpimag') ! ! Time: case('time') istat = nf_get_var1_double(ncid,i,itime,scalar) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading time') sh%time = scalar ! ! If this is a continuation run, read initial times ! Otherwise (this is an initial run) the initial time was set to the ! model start time by hist_initype. ! if (nsource <= 0) then istat=nf_get_att_int(ncid,i,"initial_year", | sh%initial_year) istat=nf_get_att_int(ncid,i,"initial_mtime", | sh%initial_mtime) endif case('mtime') ! already got it from hist search above case('year') istat = nf_get_var1_int(ncid,i,itime,iscalar) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading year') sh%year = iscalar case('day') istat = nf_get_var1_int(ncid,i,itime,iscalar) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading day') sh%day = iscalar case('iter') istat = nf_get_var1_int(ncid,i,itime,iscalar) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading iter') sh%iter = iscalar case('ut') istat = nf_get_var1_double(ncid,i,itime,scalar) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading ut') sh%ut = scalar case('dtide') istart2(1) = 1 istart2(2) = sh%ihist icount2(1) = 2 icount2(2) = 1 istat = nf_get_vara_double(ncid,i,istart2,icount2,var2) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading dtide') sh%dtide(:) = var2(:) case('sdtide') istart2(1) = 1 istart2(2) = sh%ihist icount2(1) = 10 icount2(2) = 1 istat = nf_get_vara_double(ncid,i,istart2,icount2,var10) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading sdtide') sh%sdtide(:) = var10(:) case('f107d') istat = nf_get_var1_double(ncid,i,itime,scalar) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading f107d') sh%f107d = scalar case('f107a') istat = nf_get_var1_double(ncid,i,itime,scalar) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading f107a') sh%f107a = scalar case('hpower') istat = nf_get_var1_double(ncid,i,itime,scalar) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading hpower') sh%hpower = scalar case('ctpoten') istat = nf_get_var1_double(ncid,i,itime,scalar) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading ctpoten') sh%ctpoten = scalar case('byimf') istat = nf_get_var1_double(ncid,i,itime,scalar) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading byimf') sh%byimf = scalar case('colfac') istat = nf_get_var1_double(ncid,i,itime,scalar) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading colfac') sh%colfac = scalar case('alfa30') istat = nf_get_var1_double(ncid,i,itime,scalar) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading alfa30') sh%alfa30 = scalar case('e30') istat = nf_get_var1_double(ncid,i,itime,scalar) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading e30') sh%e30 = scalar case('ed2') istat = nf_get_var1_double(ncid,i,itime,scalar) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading ed2') sh%ed2 = scalar case('alfad2') istat = nf_get_var1_double(ncid,i,itime,scalar) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading alfad2') sh%alfad2 = scalar case('p0') istat = nf_get_var_double(ncid,i,scalar) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading p0') sh%p0 = scalar case('p0_model') istat = nf_get_var_double(ncid,i,var1) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading p0_model') sh%p0_model = var1(1) case('timestep') istat = nf_get_var_int(ncid,i,iscalar) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading timestep') sh%step = iscalar ! case('nflds') ! istat = nf_get_var_int(ncid,i,iscalar) ! if (istat /= NF_NOERR) ! | call handle_ncerr(istat,'reading nflds') ! sh%nflds = iscalar ! ! GPI data file: gpi_ncfile(time,filelen), where filelen=mxlen_filename case('gpi_ncfile') ! new histories start_2d(1) = 1 start_2d(2) = sh%ihist count_2d(1) = mxlen_filename count_2d(2) = 1 istat = nf_get_vara_text(ncid,i,start_2d,count_2d, | sh%gpi_ncfile) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading gpi_ncfile') if (len_trim(sh%gpi_ncfile)==0) sh%gpi_ncfile='[none]' ! ! NCEP data file: ncep_ncfile(time,filelen), where filelen=mxlen_filename ! (ncep files used only by timegcm -- will always be '[none]' for tiegcm) case('ncep_ncfile') ! new histories start_2d(1) = 1 start_2d(2) = sh%ihist count_2d(1) = mxlen_filename count_2d(2) = 1 istat = nf_get_vara_text(ncid,i,start_2d,count_2d, | sh%ncep_ncfile) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading ncep_ncfile') if (len_trim(sh%ncep_ncfile)==0) sh%ncep_ncfile='[none]' ! ! Read path to ncep reanalysis data file (T,Z lbc): ! (as of 1/07, this is timegcm 2.5x2.5 deg grid only) case('ncep_reanalysis') start_2d(1) = 1 start_2d(2) = sh%ihist count_2d(1) = mxlen_filename count_2d(2) = 1 istat = nf_get_vara_text(ncid,i,start_2d,count_2d, | sh%ncep_reanalysis) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading ncep_reanalysis') if (len_trim(sh%ncep_reanalysis)==0) | sh%ncep_reanalysis='[none]' ! ! Read path to ECMWF ERA40 data file (T,Z lbc): ! (as of 1/07, this is timegcm 2.5x2.5 deg grid only) case('ecmwf_ncfile') start_2d(1) = 1 start_2d(2) = sh%ihist count_2d(1) = mxlen_filename count_2d(2) = 1 istat = nf_get_vara_text(ncid,i,start_2d,count_2d, | sh%ecmwf_ncfile) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading ecmwf_ncfile') if (len_trim(sh%ecmwf_ncfile)==0) | sh%ecmwf_ncfile='[none]' ! ! wacxg5 nudging files: ! case('wacxg5_ncfile') start_2d(1) = 1 start_2d(2) = sh%ihist count_2d(1) = mxlen_filename count_2d(2) = 1 istat = nf_get_vara_text(ncid,i,start_2d,count_2d, | sh%wacxg5_ncfile) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading wacxg5_ncfile') if (len_trim(sh%wacxg5_ncfile)==0) | sh%wacxg5_ncfile='[none]' ! ! GSWM data files (from namelist input): ! character(len=mxlen_filename) :: ! gswm_mi_di_ncfile, ! GSWM migrating diurnal data file ! gswm_mi_sdi_ncfile,! GSWM migrating semi-diurnal data file ! gswm_nm_di_ncfile, ! GSWM non-migrating diurnal data file ! gswm_nm_sdi_ncfile,! GSWM non-migrating semi-diurnal data file ! ! gswm_mi_di_ncfile: case('gswm_mi_di_ncfile') istat = nf_get_vara_text(ncid,i,start_2d,count_2d, | sh%gswm_mi_di_ncfile) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading gswm_mi_di_ncfile') ! ! gswm_mi_sdi_ncfile: case('gswm_mi_sdi_ncfile') istat = nf_get_vara_text(ncid,i,start_2d,count_2d, | sh%gswm_mi_sdi_ncfile) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading gswm_mi_sdi_ncfile') ! ! gswm_nm_di_ncfile: case('gswm_nm_di_ncfile') istat = nf_get_vara_text(ncid,i,start_2d,count_2d, | sh%gswm_nm_di_ncfile) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading gswm_nm_di_ncfile') ! ! gswm_nm_sdi_ncfile: case('gswm_nm_sdi_ncfile') istat = nf_get_vara_text(ncid,i,start_2d,count_2d, | sh%gswm_nm_sdi_ncfile) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading gswm_nm_sdi_ncfile') case('maglon') case('maglat') ! ! Interface level of lower boundary (scalar coord for TLBC,ULBC,VLBC): case('LBC') istat = nf_get_var1_double(ncid,i,itime,varlbc) if (istat /= NF_NOERR) | call handle_ncerr(istat,'reading LBC') if (varlbc /= zibot) then write(6,"(/,'>>> WARNING nc_rdhist read LBC: varlbc=', | f8.2,' but should equal zibot=',f8.2)") varlbc,zibot else sh%lbc = varlbc write(6,"('Read LBC from source history: i=',i3, | ' sh%lbc=',f8.2)") i,sh%lbc endif ! ! Read lbc's for t,u,v. See fields.F. These will exist on source history ! only if sh%tuv_lbc_intop==0, i.e., "new histories"). Read into global ! arrays, then store subdomains. ! case('TLBC') call rdlbc(ncid,i,tlbc,varname,itime) case('ULBC') call rdlbc(ncid,i,ulbc,varname,itime) case('VLBC') call rdlbc(ncid,i,vlbc,varname,itime) case('TLBC_NM') call rdlbc(ncid,i,tlbc_nm,varname,itime) case('ULBC_NM') call rdlbc(ncid,i,ulbc_nm,varname,itime) case('VLBC_NM') call rdlbc(ncid,i,vlbc_nm,varname,itime) ! ! Check for variable name in 4-d fields: case default if (any(f4d(1:nf4d_hist)%short_name==varname)) then do j=1,nf4d_hist if (trim(varname)==f4d(j)%short_name) | call rdf4d(ncid,varname,itime,j) enddo ! ! 10/27/05 btf: This conditional covers "W" from old histories. ! (var name "W" has been changed to "OMEGA" on the new histories) ! elseif (varname=='W') then do j=1,nf4d_hist if (f4d(j)%short_name=='OMEGA') | call rdf4d(ncid,varname,itime,j) enddo ! ! Variable is unknown to the model: else write(6,"('Note nc_rdhist: unused variable: ',a)") | trim(varname) endif end select enddo ! i=1,nvars sh%zmtop = zmtop sh%zmbot = zmbot sh%zitop = zitop sh%zibot = zibot ! ! Report to stdout: call hist_print(sh,'READ',diskfile) ! end subroutine nc_rdhist !----------------------------------------------------------------------- subroutine rdlbc(ncid,idv,flbc,vname,itime) ! ! Read global TLBC, ULBC, or VLBC from current history, and transfer to ! subdomains tlbc, ulbc, or vlbc (fields.F). ! use mpi_module,only: lon0,lon1,lat0,lat1 use fields_module,only: lond0,lond1,latd0,latd1 ! Routine to add fields to secondary histories: ! use addfld_module,only: addfld implicit none ! ! Args: integer,intent(in) :: ncid,idv,itime character(len=*),intent(in) :: vname real,intent(out) :: flbc(lond0:lond1,latd0:latd1) ! ! Local: integer :: istat,i,j,lonbeg,lonend real,dimension(nlon,nlat) :: glbc ! lbc at global domain character(len=80) :: char80 ! start_3d(1:2) = 1 start_3d(3) = itime count_3d(1) = nlon count_3d(2) = nlat count_3d(3) = 1 ! ! Read global domain into local glbc: istat = nf_get_vara_double(ncid,idv,start_3d,count_3d,glbc) if (istat /= NF_NOERR) then write(char80,"('Error reading ',a,' ncid=',i3)") | trim(vname),ncid call handle_ncerr(istat,trim(char80)) endif ! ! Store subdomain in flbc (actual arg is tlbc, ulbc, or vlbc): lonbeg = lon0 if (lon0==1) lonbeg = 3 lonend = lon1 if (lon1==nlonp4) lonend = lon1-2 do j=lat0,lat1 do i=lonbeg,lonend flbc(i,j) = glbc(i-2,j) enddo enddo #ifndef MPI ! Periodic points for non-MPI runs: flbc(1:2,:) = flbc(nlonp4-3:nlonp4-2,:) flbc(nlonp4-1:nlonp4,:) = flbc(3:4,:) #endif write(6,"('Read ',a,' subdomain from source history: min,max=', | 2f8.2)") trim(vname), minval(flbc(lon0:lon1,lat0:lat1)), | maxval(flbc(lon0:lon1,lat0:lat1)) ! if (trim(vname)=='TLBC') then ! do i=lon0,lon1 ! write(6,"('rdlbc: i=',i3,' lat0,1=',2i3,' tlbc(i,lat0:lat1)=', ! | /,(6e12.4))") i,lat0,lat1,tlbc(i,lat0:lat1) ! enddo ! j=lat0,lat1 ! endif ! if (trim(vname)=='ULBC') then ! do i=lon0,lon1 ! write(6,"('rdlbc: i=',i3,' lat0,1=',2i3,' ulbc(i,lat0:lat1)=', ! | /,(6e12.4))") i,lat0,lat1,ulbc(i,lat0:lat1) ! enddo ! j=lat0,lat1 ! endif ! if (trim(vname)=='VLBC') then ! do i=lon0,lon1 ! write(6,"('rdlbc: i=',i3,' lat0,1=',2i3,' vlbc(i,lat0:lat1)=', ! | /,(6e12.4))") i,lat0,lat1,vlbc(i,lat0:lat1) ! enddo ! j=lat0,lat1 ! endif end subroutine rdlbc !----------------------------------------------------------------------- subroutine rdf4d(ncid,name,itime,ix) ! ! Read prognostic field "name" from netcdf file on ncid into into local ! f3drd at full domain. Then subdomains are transferred from f3rd into ! 4-d arrays for use in the model. It would be interesting to try reading ! only the subdomain directly from the file rather than the full domain. ! use input_module,only: dynamo use hist_module,only: sh use init_module,only: glat use mpi_module,only: lon0,lon1,lat0,lat1 use fields_module,only: dynpot ! ! Args: character(len=*),intent(in) :: name integer,intent(in) :: ncid,itime,ix ! ! Local: integer :: i,ii,j,k,istat,idvar,itype,ndims,iddims(4),natts, | iscalar,idimsizes(4),nx,nxk,lonbeg,lonend,nnans,nanfatal character(len=16) :: rdname,tmpname character(len=120) :: char120 real :: fakevar(1,1,1),fmin,fmax,fmin1,fmax1 real :: f3drd(nlon,nlat,nlevp1) ! (i,j,k) as on history real :: zmj2,zmj1,fac ! ! External: integer,external :: strloc ! ! Get info about the field: istat = nf_inq_varid(ncid,name,idvar) istat = nf_inq_var(ncid,idvar,rdname,itype,ndims,iddims,natts) if (ndims /= 4) then write(6,"(/,'>>> WARNING rdf4d: bad ndims=',i3, | ' (every prognostic should have 4 dimensions)')") ndims endif ! ! Get info about the dimensions: do i=1,ndims istat = nf_inq_dim(ncid,iddims(i),tmpname,idimsizes(i)) enddo ! ! Check dimension sizes against current model grid resolution: if (.not.fakeflds) then if (idimsizes(1) /= nlon .or. idimsizes(2) /= nlat .or. | idimsizes(3) /= nlevp1) then write(6,"(/,'>>> WARNING rdf4d: bad dimension sizes', | ' for field ',a)") trim(rdname) write(6,"(' idimsizes(1:3)=',3i4,' but should be ', | 'nlon,nlat,nlevp1=',3i4)") idimsizes(1:3),nlon,nlat, | nlevp1 call shutdown('source history coordinate dimensions') endif endif ! if (fakeflds) then start_4d(1:3) = 1 start_4d(4) = itime count_4d(1) = 1 count_4d(2) = 1 count_4d(3) = 1 count_4d(4) = 1 istat = nf_get_vara_double(ncid,idvar,start_4d,count_4d,fakevar) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_get_vara_double', + ' for fake field var ',a,' itime=',i2)") name,itime call handle_ncerr(istat,char120) endif return endif ! ! Zero out electric potential and return if DYNAMO flag is not set: if (trim(name)=='POTEN'.and.dynamo <= 0) then f4d(ix)%data = 0. return endif ! ! Read field into data component of field structure: start_4d(1:3) = 1 start_4d(4) = itime count_4d(1) = nlon count_4d(2) = nlat count_4d(3) = nlevp1 count_4d(4) = 1 istat = nf_get_vara_double(ncid,idvar,start_4d,count_4d,f3drd) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_get_vara_double', + ' for field var ',a,' itime=',i2)") name,itime call handle_ncerr(istat,char120) endif ! ! Check for NaNs in the source history: ! subroutine check_nans(f,id1,id2,id3,name,n_total,ispval,spval, ! | iprint,ifatal) ! real :: f3drd(nlon,nlat,nlevp1) ! (i,j,k) as on history ! ! nanfatal = 0 ! if 0, then reading nans from source history is non-fatal nanfatal = 1 ! if 1, then reading nans from source history is fatal if (check_nan) | call check_nans(f3drd,nlon,nlat,nlevp1,name(1:8),nnans,0,0.,1, | nanfatal) ! ! Transform from f3drd(i,j,k) to f4d(k,i,j): lonbeg = lon0 if (lon0==1) lonbeg = 3 lonend = lon1 if (lon1==nlonp4) lonend = lon1-2 do j=lat0,lat1 do i=lonbeg,lonend do k=1,nlevp1 ! ! Local : real :: f3drd(nlon,nlat,nlevp1) ! (i,j,k) as on history ! Global: real,pointer :: f4d(ix)%data(levd0:levd1,lond0:lond1,latd0:latd1) ! f4d(ix)%data(k,i,j,itp) = f3drd(i-2,j,k) enddo ! k=1,nlevp1 enddo ! i=lonbeg,lonend enddo ! j=lat0,lat1 ! ! Define t,u,v lbc subdomains from top slot if tuv_lbc_intop==1: ! If tuv_lbc_intop==0, these arrays were read by nc_rdhist. ! These arrays are in fields.F, allocated by sub init_lbc. ! if (sh%tuv_lbc_intop > 0) then ! (is "old" history) select case(trim(name)) case('TN') do j=lat0,lat1 tlbc(:,j) = f4d(ix)%data(nlevp1,:,j,itp) enddo ! j=lat0,lat1 write(6,"('rdf4d: tlbc from top slot min,max=',2f10.2)") | minval(tlbc(lon0:lon1,lat0:lat1)), | maxval(tlbc(lon0:lon1,lat0:lat1)) f4d(ix)%data(nlevp1,:,:,itp) = spval ! top slot is now missing data case('UN') do j=lat0,lat1 ulbc(:,j) = f4d(ix)%data(nlevp1,:,j,itp) enddo write(6,"('rdf4d: ulbc from top slot min,max=',2f10.2)") | minval(ulbc(lon0:lon1,lat0:lat1)), | maxval(ulbc(lon0:lon1,lat0:lat1)) f4d(ix)%data(nlevp1,:,:,itp) = spval ! top slot is now missing data case('VN') do j=lat0,lat1 vlbc(:,j) = f4d(ix)%data(nlevp1,:,j,itp) enddo write(6,"('rdf4d: vlbc from top slot min,max=',2f10.2)") | minval(vlbc(lon0:lon1,lat0:lat1)), | maxval(vlbc(lon0:lon1,lat0:lat1)) f4d(ix)%data(nlevp1,:,:,itp) = spval ! top slot is now missing data case('TN_NM') do j=lat0,lat1 tlbc_nm(:,j) = f4d(ix)%data(nlevp1,:,j,itp) enddo ! j=lat0,lat1 write(6,"('rdf4d: tlbc_nm from top slot min,max=',2f10.2)") | minval(tlbc_nm(lon0:lon1,lat0:lat1)), | maxval(tlbc_nm(lon0:lon1,lat0:lat1)) f4d(ix)%data(nlevp1,:,:,itp) = spval ! top slot is now missing data case('UN_NM') do j=lat0,lat1 ulbc_nm(:,j) = f4d(ix)%data(nlevp1,:,j,itp) enddo write(6,"('rdf4d: ulbc_nm from top slot min,max=',2f10.2)") | minval(ulbc_nm(lon0:lon1,lat0:lat1)), | maxval(ulbc_nm(lon0:lon1,lat0:lat1)) f4d(ix)%data(nlevp1,:,:,itp) = spval ! top slot is now missing data case('VN_NM') do j=lat0,lat1 vlbc_nm(:,j) = f4d(ix)%data(nlevp1,:,j,itp) enddo write(6,"('rdf4d: vlbc_nm from top slot min,max=',2f10.2)") | minval(vlbc_nm(lon0:lon1,lat0:lat1)), | maxval(vlbc_nm(lon0:lon1,lat0:lat1)) f4d(ix)%data(nlevp1,:,:,itp) = spval ! top slot is now missing data case default ! do nothing if not one of these fields end select ! case(trim(name)) endif ! if (sh%tuv_lbc_intop > 0) ! ! Report to stdout: ! call fminmax(f4d(ix)%data(:,lonbeg:lonend,lat0:lat1,itp), | nlevp1*(lonend-lonbeg+1)*(lat1-lat0+1),fmin,fmax) ! write(6,"('lonbeg,lonend=',2i3)") lonbeg,lonend write(6,"('Read field ',a,' 3d min,max=',2e12.4)") | f4d(ix)%short_name,fmin,fmax ! ! For MPI runs, subdomains of poten are read from the source history above, ! then mp_dynpot is called from readsource to define dynpot at the full ! domain for each task. For non-MPI runs, define dynpot here from source. ! #ifdef MPI #else ! Serial non-MPI: ! Define dynpot(nlonp1,0:nlatp1,nlevp1) from electric potential ! if (trim(name)=='POTEN') then write(6,"('rdf4d defining dynpot from f4d..')") do k=1,nlevp1 do j=lat0,lat1 ! do i=3,nlonp1+1 ! dynpot out-of-bounds do i=3,nlonp1 dynpot(i,j,k) = f3drd(i-2,j,k) enddo ! k=1,nlevp1 enddo ! i=lonbeg,lonend enddo ! j=lat0,lat1 endif #endif ! end subroutine rdf4d !----------------------------------------------------------------------- subroutine nc_define(ncid) ! ! Define dimensions, variables and attributes in a new history ! netcdf dataset (file has been opened, but no histories written yet). ! On input: ! ncid = file id for open netcdf dataset (in define mode). ! h = defined history structure (in hist_module) ! On output: ! Variables and dimensions for the netcdf file are defined. ! Values are given to the dimension variables, and the file ! is taken out of define mode (ready to receive histories) ! use hist_module,only: sh,h,nsource,nsecsource use init_module,only: glon,glat,gmlon,gmlat use cons_module,only: pi,ylatmagphr,ylonmagphr use input_module,only: potential_model,ncep_ncfile, | ncep_reanalysis,ecmwf_ncfile use params_module,only: tgcm_name,tgcm_version ! ! Args: integer,intent(in) :: ncid ! ! Local: integer :: i,istat,idum,ivar1(1),imo,ida,startmtime(4),nbytes character(len=120) :: char120 character(len=80) :: char80 character(len=24) :: create_date,create_time real :: var1(1),rdays real,external :: mtime_to_datestr ! ! Define dimensions (field vars will be (time,lev,lat,lon) in ! CDL notation. Mag fields will be (time,mlev,mlat,mlon): ! istat = nf_def_dim(ncid,"time",NF_UNLIMITED,id_time) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining time dimension') ! ! Geographic grid dimensions: istat = nf_def_dim(ncid,"lon",nlon,id_lon) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining longitude dimension') istat = nf_def_dim(ncid,"lat",nlat,id_lat) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining latitude dimension') ! ! Midpoint levels dimension "lev"=nlevp1: istat = nf_def_dim(ncid,"lev",nlevp1,id_lev) ! midpoint levels dimension if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining midpoint levels dimension') ! ! Interface levels dimension "ilev"=nlevp1: istat = nf_def_dim(ncid,"ilev",nlevp1,id_ilev) ! interface levels dimension if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining interface levels dimension') ! ! Magnetic grid dimensions: istat = nf_def_dim(ncid,"mlon",nmlonp1,id_mlon) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining magnetic longitude dimension') istat = nf_def_dim(ncid,"mlat",nmlat,id_mlat) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining magnetic latitude dimension') ! ! Mag midpoint levels dimension: istat = nf_def_dim(ncid,"mlev",nmlevp1,id_mlev) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining magnetic midpoints levels dimension') ! ! Mag interface levels dimension: istat = nf_def_dim(ncid,"imlev",nmlevp1,id_imlev) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining magnetic interface levels dimension') ! ! Msc dimensions mtimedim, latlon, dtidedim, sdtidedim: istat = nf_def_dim(ncid,"mtimedim",3,id_mtimedim) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining mtime dimension') istat = nf_def_dim(ncid,"latlon",2,id_latlon) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining latlon dimension') istat = nf_def_dim(ncid,"dtidedim",2,id_dtide) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining dtide dimension') istat = nf_def_dim(ncid,"sdtidedim",10,id_sdtide) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining sdtide dimension') if (fakeflds) then ! fake dimension for testing istat = nf_def_dim(ncid,"fakedim",1,id_fake) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining fakedim dimension') endif ! ! String length of date and time: istat = nf_def_dim(ncid,"datelen",len(h%write_date),id_datelen) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining datelen dimension') ! ! String length of file names (mxlen_filename is parameter in input.F): istat = nf_def_dim(ncid,"filelen",mxlen_filename,id_filelen) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining filelen dimension') ! ! Define dimension (coordinate) variables and attributes: ! ! Time (coordinate variable time(time)). This is days since ! the initial run's start time. The units string is: yyyy-m-d, ! where yyyy is the year, m is month, and d is day of the source ! start time. ! ids1(1) = id_time ! for 1d time-unlimited vars 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') startmtime(1:3) = sh%initial_mtime(1:3) ; startmtime(4) = 0 rdays = mtime_to_datestr(sh%initial_year,startmtime,imo,ida, | char80) istat = nf_put_att_text(ncid,idv_time,"long_name",4,"time") if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining long_name of time dimension variable') istat = nf_put_att_text(ncid,idv_time,"units",len_trim(char80), | trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining units of time dimension variable') istat = nf_put_att_int(ncid,idv_time,"initial_year",NF_INT,1, | h%initial_year) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining initial_year attribute of time coord var') istat = nf_put_att_int(ncid,idv_time,"initial_mtime",NF_INT,3, | h%initial_mtime) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining initial_mtime attribute of time coord var') ! ! Geographic longitude (deg) (coordinate variable lon(lon)): 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)) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining long_name of longitude dimension variable') istat = nf_put_att_text(ncid,idv_lon,"units",12,'degrees_east') if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining units of longitude dimension variable') ! ! Geographic latitude (deg) (coordinate variable lat(lat)): 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)) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining long_name of latitude dimension variable') istat = nf_put_att_text(ncid,idv_lat,"units",13,'degrees_north') if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining units of latitude dimension variable') ! ! Midpoint levels coordinate variable lev(lev): ids1(1) = id_lev istat = nf_def_var(ncid,"lev",NF_DOUBLE,1,ids1,idv_lev) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining midpoint levels coordinate variable lev(lev)') ! long name of lev coord var: write(char80,"('midpoint levels')") istat = nf_put_att_text(ncid,idv_lev,"long_name", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, |'Error defining long_name of midpoint levels coordinate variable') ! short name = ln(p0/p) istat = nf_put_att_text(ncid,idv_lev,"short_name", | 8,"ln(p0/p)") if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining short_name of midpoint levels coord var') ! lev coord var is unitless: istat = nf_put_att_text(ncid,idv_lev,"units",0," ") if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining units of midpoint levels coord var') ! positive='up' istat = nf_put_att_text(ncid,idv_lev,"positive",2,'up') if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining positive attribute of midpoint levels coord var') ! standard name of midpoints lev coord var: write(char80,"('atmosphere_ln_pressure_coordinate')") istat = nf_put_att_text(ncid,idv_lev,"standard_name", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining standard_name of midpoint levels coord var') ! formula terms for midpoints lev coord var: write(char80,"('p0:p0 lev:lev')") istat = nf_put_att_text(ncid,idv_lev,"formula_terms", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining formula terms of midpoint levels coord var') ! formula to obtain pressure from midpoint lev: write(char80,"('p(k) = p0 * exp(-lev(k))')") istat = nf_put_att_text(ncid,idv_lev,"formula", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining formula for midpoint levels coord var') ! ! Interface levels coordinate array: ids1(1) = id_ilev istat = nf_def_var(ncid,"ilev",NF_DOUBLE,1,ids1,idv_ilev) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining interface levels coordinate variable ilev(ilev)') ! long name of ilev coord var: write(char80,"('interface levels')") istat = nf_put_att_text(ncid,idv_ilev,"long_name", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, |'Error defining long_name of interface levels coordinate var') ! short name = ln(p0/p) istat = nf_put_att_text(ncid,idv_ilev,"short_name", | 8,"ln(p0/p)") if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining short_name of interface levels coord var') ! ilev coord var is unitless: istat = nf_put_att_text(ncid,idv_ilev,"units",0," ") if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining units of interface levels coord var') ! positive='up' istat = nf_put_att_text(ncid,idv_ilev,"positive",2,'up') if (istat /= NF_NOERR) call handle_ncerr(istat, |'Error defining positive attribute of interface levels coord var') ! standard name of interface ilev coord var: write(char80,"('atmosphere_ln_pressure_coordinate')") istat = nf_put_att_text(ncid,idv_ilev,"standard_name", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining standard_name of interface levels coord var') ! formula terms for interface ilev coord var: write(char80,"('p0:p0 ilev:ilev')") istat = nf_put_att_text(ncid,idv_ilev,"formula_terms", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining formula terms of interface levels coord var') ! formula to obtain pressure from interface ilev: write(char80,"('p(k) = p0 * exp(-ilev(k))')") istat = nf_put_att_text(ncid,idv_ilev,"formula", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining formula for interface levels coord var') ! ! Magnetic longitude (deg) (coordinate variable mlon(mlon)): ids1(1) = id_mlon istat = nf_def_var(ncid,"mlon",NF_DOUBLE,1,ids1,idv_mlon) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining magnetic longitude dimension variable') write(char80,"('magnetic longitude (-west, +east)')") istat = nf_put_att_text(ncid,idv_mlon,"long_name", + len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining long_name of mag longitude dimension variable') istat = nf_put_att_text(ncid,idv_mlon,"units",12,'degrees_east') if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining units of magnetic longitude dimension variable') ! ! Magnetic latitude (deg) (coordinate variable mlat(mlat)): ids1(1) = id_mlat istat = nf_def_var(ncid,"mlat",NF_DOUBLE,1,ids1,idv_mlat) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining magnetic latitude dimension variable') write(char80,"('magnetic latitude (-south, +north)')") istat = nf_put_att_text(ncid,idv_mlat,"long_name", + len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining long_name of mag latitude dimension variable') istat = nf_put_att_text(ncid,idv_mlat,"units",13,'degrees_north') if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining units of magnetic latitude dimension variable') ! ! Magnetic midpoint coordinate variable mlev(mlev)): ids1(1) = id_mlev istat = nf_def_var(ncid,"mlev",NF_DOUBLE,1,ids1,idv_mlev) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining magnetic levels coord variable mlev') ! long name of mlev: write(char80,"('magnetic midpoint levels')") istat = nf_put_att_text(ncid,idv_mlev,"long_name", + len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining long_name of mlev coordinate variable') ! short name = ln(p0/p) istat = nf_put_att_text(ncid,idv_mlev,"short_name", | 8,"ln(p0/p)") if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining short_name of mlev levels coord var') ! units of mlev (unitless): istat = nf_put_att_text(ncid,idv_mlev,"units",0," ") if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining units of mlev coordinate variable') ! positive attribute of mlev: istat = nf_put_att_text(ncid,idv_mlev,"positive",2,'up') if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining positive attribute for mlev coord var') ! standard_name of mlev: write(char80,"('atmosphere_ln_pressure_coordinate')") istat = nf_put_att_text(ncid,idv_mlev,"standard_name", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining standard_name attribute for mlev coord var') ! formula terms for mlev: write(char80,"('p0:p0 mlev:mlev')") istat = nf_put_att_text(ncid,idv_mlev,"formula_terms", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining formula terms of mlev coord var') ! formula to obtain pressure from mlev: write(char80,"('p(k) = p0 * exp(-mlev(k))')") istat = nf_put_att_text(ncid,idv_mlev,"formula", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining formula for mlev coord var') ! ! Magnetic interface coordinate variable imlev(imlev)): ids1(1) = id_imlev istat = nf_def_var(ncid,"imlev",NF_DOUBLE,1,ids1,idv_imlev) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining magnetic levels coord variable imlev') ! long name of imlev: write(char80,"('magnetic interface levels')") istat = nf_put_att_text(ncid,idv_imlev,"long_name", + len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining long_name of imlev coordinate variable') ! short name = ln(p0/p) istat = nf_put_att_text(ncid,idv_imlev,"short_name", | 8,"ln(p0/p)") if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining short_name of imlev levels coord var') ! units of imlev (unitless): istat = nf_put_att_text(ncid,idv_imlev,"units",0," ") if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining units of imlev coordinate variable') ! positive attribute of imlev: istat = nf_put_att_text(ncid,idv_imlev,"positive",2,'up') if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining positive attribute for imlev coord var') ! standard_name of imlev: write(char80,"('atmosphere_ln_pressure_coordinate')") istat = nf_put_att_text(ncid,idv_imlev,"standard_name", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining standard_name attribute for imlev coord var') ! formula terms for imlev: write(char80,"('p0:p0 imlev:imlev')") istat = nf_put_att_text(ncid,idv_imlev,"formula_terms", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining formula terms of imlev coord var') ! formula to obtain pressure from imlev: write(char80,"('p(k) = p0 * exp(-imlev(k))')") istat = nf_put_att_text(ncid,idv_imlev,"formula", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining formula for imlev coord var') ! ! ids1(1) is id of time dimension for the following defines: ids1(1) = id_time ! for 1d time-unlimited vars ! ! Source date and time are attributes of time var ! (they are also global file attributes): istat = nf_put_att_int(ncid,idv_time,"initial_year",NF_INT,1, | h%initial_year) istat = nf_put_att_int(ncid,idv_time,"initial_day",NF_INT,1, | h%initial_day) istat = nf_put_att_int(ncid,idv_time,"initial_mtime",NF_INT,3, | h%initial_mtime) ! ! Define "primary" (non-coordinate) time dependent variables: ! (these variables get a value for each history) ! ! Model time variable mtime(3,time), where 3 is the mtime dimension ! (id_mtimedim) (day,hour,minute), and time is the unlimited dimension ! (id_time): ! ids2(1) = id_mtimedim ids2(2) = id_time istat = nf_def_var(ncid,"mtime",NF_INT,2,ids2,idv_mtime) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining mtime variable') write(char80,"('model times (day, hour, minute)')") istat = nf_put_att_text(ncid,idv_mtime,"long_name", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining long_name attribute of mtime variable') istat = nf_put_att_text(ncid,idv_mtime,"units",17, | 'day, hour, minute') if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining units of mtime variable') ! ! Calendar year and day (this will be advanced if model is advancing ! in calendar time, otherwise will be constant from input) istat = nf_def_var(ncid,"year",NF_INT,1,ids1,idv_year) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining calendar year variable') istat = nf_put_att_text(ncid,idv_year,"long_name", | 13,'calendar year') if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining long_name of calendar year variable') ! istat = nf_def_var(ncid,"day",NF_INT,1,ids1,idv_day) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining calendar day variable') istat = nf_put_att_text(ncid,idv_day,"long_name", | 12,'calendar day') if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining long_name of calendar day variable') ! ! Calendar advance: istat = nf_def_var(ncid,"calendar_advance",NF_INT,1,ids1, | idv_calendar_adv) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining integer calendar advance variable') write(char80,"('calendar advance flag', | ' (1 if advancing calendar time)')") istat = nf_put_att_text(ncid,idv_calendar_adv,"long_name", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining long_name of calendar advance variable') ! ! Date each history was written: ids2(1) = id_datelen ids2(2) = id_time istat= nf_def_var(ncid,"write_date",NF_CHAR,2,ids2,idv_writedate) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining history write_date variable') write(char80,"('Date and time each history was written')") istat = nf_put_att_text(ncid,idv_writedate,"long_name", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining long_name of write_date') ! ! Number of time steps from model time 0,0,0 to current model time ! (iter(time)): istat = nf_def_var(ncid,"iter",NF_INT,1,ids1,idv_iter) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining iter variable') istat = nf_put_att_text(ncid,idv_iter,"long_name", | 42,'number of time steps from model time 0,0,0') if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining long_name of iter variable') ! ! Decimal ut(time) (ids1==id_time is unlimited time dimension): call defvar_time_dbl(ncid,"ut",ids1,idv_ut, | "universal time (from model time hour and minute)", | "hours") ! ! Time step (seconds) (integer scalar) (time-independent): istat = nf_def_var(ncid,"timestep",NF_INT,0,idum,idv_step) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining timestep variable') istat = nf_put_att_text(ncid,idv_step,"long_name",8,'timestep') if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining long_name of timestep variable') istat = nf_put_att_text(ncid,idv_step,"units",7,'seconds') if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining units of timestep variable') ! ! f107 daily and f107 average solar fluxes: call defvar_time_dbl(ncid,"f107d",ids1,idv_f107d, | "Daily 10.7 cm solar flux","1.e-22 W/m2/Hz") call defvar_time_dbl(ncid,"f107a",ids1,idv_f107a, | "81-day average 10.7 cm solar flux","1.e-22 W/m2/Hz") ! ! Hemispheric power: call defvar_time_dbl(ncid,"hpower",ids1,idv_hpower, | "Hemispheric power","GW") ! ! Cross-tail potential: call defvar_time_dbl(ncid,"ctpoten",ids1,idv_ctpoten, | "Cross-tail potential","kV") ! ! Byimf: call defvar_time_dbl(ncid,"byimf",ids1,idv_byimf, | "BY component of IMF"," ") ! ! Collision factor: call defvar_time_dbl(ncid,"colfac",ids1,idv_colfac, | "ion/neutral collision factor"," ") ! ! 1/31/01: Making tides time-dependent (dtide(time,2) and sdtide(time,10): ! Diurnal tide amp/phase: ids2(1) = id_dtide ids2(2) = id_time istat = nf_def_var(ncid,"dtide",NF_DOUBLE,2,ids2,idv_dtide) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining variable dtide') write(char80,"('amplitude and phase of diurnal tide')") istat = nf_put_att_text(ncid,idv_dtide,"long_name", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining long_name attribute of variable dtide') ! ! Semi-diurnal tide amp/phase: ids2(1) = id_sdtide ids2(2) = id_time istat = nf_def_var(ncid,"sdtide",NF_DOUBLE,2,ids2,idv_sdtide) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining variable sdtide') write(char80,"('amplitude and phase of semi-diurnal tide')") istat = nf_put_att_text(ncid,idv_sdtide,"long_name", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining long_name attribute of variable sdtide') ! ! Define ncep data file (used by timegcm only): call define_ncfile(ncid,"ncep_ncfile",idv_ncep_ncfile, | "long_name","NCEP data file", | "data","10 mb lower boundary for TN and Z") ! ! Define ncep reanalysis data file (used by timegcm only): ! (Only at 2.5 deg grid, as of Jan, 07) call define_ncfile(ncid,"ncep_reanalysis",idv_ncep_reanalysis, | "long_name","NCEP reanalysis data file", | "data","10 mb lower boundary for TN and Z") ! ! Define ecmwf nc data file (used by timegcm only): ! (Only at 2.5 deg grid, as of Jan, 07) call define_ncfile(ncid,"ecmwf_ncfile",idv_ecmwf_ncfile, | "long_name","ECMWF ERA40 data file", | "data","10 mb lower boundary for TN and Z") ! ! Define wacxg5 nudging file: call define_ncfile(ncid,"wacxg5_ncfile",idv_wacxg5_ncfile, | "long_name","WACCM-X nudging file","data", | "wacxg5 file used for nudging upward from lower boundary") ! ! Define gpi data file: call define_ncfile(ncid,"gpi_ncfile",idv_gpi_ncfile, | "long_name","GeoPhysical Indices data file", | "data","f107, f107a, ctpoten, power") ! ! Define gswm data files: call define_ncfile(ncid,"gswm_mi_di_ncfile", | idv_gswm_mi_di_ncfile, | "long_name","GSWM migrating diurnal tides data file", | "data","lower boundary perturbations for TN, UN, VN, and Z") call define_ncfile(ncid,"gswm_mi_sdi_ncfile", | idv_gswm_mi_sdi_ncfile, | "long_name","GSWM migrating semi-diurnal tides data file", | "data","lower boundary perturbations for TN, UN, VN, and Z") call define_ncfile(ncid,"gswm_nm_di_ncfile", | idv_gswm_nm_di_ncfile, | "long_name","GSWM non-migrating diurnal tides data file", | "data","lower boundary perturbations for TN, UN, VN, and Z") call define_ncfile(ncid,"gswm_nm_sdi_ncfile", | idv_gswm_nm_sdi_ncfile, | "long_name","GSWM non-migrating semi-diurnal tides data file", | "data","lower boundary perturbations for TN, UN, VN, and Z") ! ! 2/14/01: Add auroral parameters (from input): ! ! alfa30: call defvar_time_dbl(ncid,"alfa30",ids1,idv_alfa30, | "Characteristic energy of high-energy auroral electrons", | "KeV") ! ! e30: call defvar_time_dbl(ncid,"e30",ids1,idv_e30, | "Energy flux of high-energy auroral electrons", | "ergs/cm2/s") ! ! alfad2: call defvar_time_dbl(ncid,"alfad2",ids1,idv_alfad2, | "Characteristic energy of solar protons in the polar cap", | "KeV") ! ! ed2: call defvar_time_dbl(ncid,"ed2",ids1,idv_ed2, | "Energy flux of solar protons in the polar cap", | "ergs/cm2/s") ! ! Time-independent "primary" variables: ! ! Magnetic pole coordinates (real(2,2)): ids2(1) = id_latlon ids2(2) = id_latlon istat = nf_def_var(ncid,"mag",NF_DOUBLE,2,ids2,idv_mag) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining variable mag') write(char80,"('lat,lon coordinates of S,N magnetic poles')") istat = nf_put_att_text(ncid,idv_mag,"long_name", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining long_name attribute of variable mag') istat = nf_put_att_text(ncid,idv_mag,"units",7, | 'degrees') if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining units of variable mag') ! ! Reference pressure used to convert ln(p0/p) to pressure (millibars (hPa)): istat = nf_def_var(ncid,"p0",NF_DOUBLE,0,idum,idv_p0) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining variable p0') istat = nf_put_att_text(ncid,idv_p0,"long_name", | 18,"Reference pressure") ! for conversion of vertical coordinates ! ! 4/26/06 btf: p0 units can be either hPa or millibars. Using the latter ! now to be consistent w/ microbar units of p0_model istat = nf_put_att_text(ncid,idv_p0,"units",9,"millibars") ! istat = nf_put_att_text(ncid,idv_p0,"units",3,"hPA") ! ! Reference pressure used by the model (dynes/cm2): ! p0_model units can be either dynes/cm2 or microbars. Using the latter ! to be consistent w/ 1981 paper by Dickinson, Ridley, and Roble ! istat = nf_def_var(ncid,"p0_model",NF_DOUBLE,0,idum,idv_p0_model) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining variable p0') char80 = ' ' write(char80,"('Reference pressure (as used by the model)')") istat = nf_put_att_text(ncid,idv_p0_model,"long_name", | len_trim(char80),trim(char80)) istat = nf_put_att_text(ncid,idv_p0_model,"units",9,"microbars") ! istat = nf_put_att_text(ncid,idv_p0_model,"units",9,"dynes/cm2") ! ! Gravity (constant used in the model, independent of height): istat = nf_def_var(ncid,"grav",NF_DOUBLE,0,idum,idv_grav) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error defining variable grav') istat = nf_put_att_text(ncid,idv_grav,"long_name", | 26,"gravitational acceleration") istat = nf_put_att_text(ncid,idv_grav,"units", | 4,"cm/s") write(char80,"('constant used in the model, independent of ', | 'height')") istat = nf_put_att_text(ncid,idv_grav,"info", | len_trim(char80),trim(char80)) ! ! All prognostic fields are on primary history: if (h%hist_type(1:3)=='pri') then nbytes = 8 do i=1,nf4d call defvar_f4d(ncid,f4d(i),idv_f4d(i),nbytes) enddo ! ! Define secondary history fields (geographic, magnetic and magnetospheric): elseif (h%hist_type(1:3)=='sec') then call def_fsech(ncid) else write(6,"('>>> nc_define: unknown h%hist_type=',a)") h%hist_type endif ! ! Define lower boundaries of t,u,v (both secondary and primary): call deflbc(ncid,"TLBC",idv_tlbc,h%hist_type) call deflbc(ncid,"ULBC",idv_ulbc,h%hist_type) call deflbc(ncid,"VLBC",idv_vlbc,h%hist_type) call deflbc(ncid,"TLBC_NM",idv_tlbc_nm,h%hist_type) call deflbc(ncid,"ULBC_NM",idv_ulbc_nm,h%hist_type) call deflbc(ncid,"VLBC_NM",idv_vlbc_nm,h%hist_type) ! ! Define LBC: bottom interface level (scalar coord of TLBC,ULBC,VLBC): istat = nf_def_var(ncid,"LBC",NF_DOUBLE,0,idum,idv_lbc) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining var LBC') write(char80, | "('Interface level of t,u,v lower boundary condition')") istat = nf_put_att_text(ncid,idv_lbc,"long_name", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining long_name attribute of variable LBC') ! write(6,"('defined LBC: idv_lbc=',i3)") idv_lbc ! ! Define global file attributes: ! ! User-defined job comment (namelist parameter "label"): if (len_trim(h%label) > 0) then istat = nf_put_att_text(ncid,NF_GLOBAL,"label", | len_trim(h%label),trim(h%label)) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_text ', | 'for label global attribute')") call handle_ncerr(istat,char120) endif endif ! ! Create date and time: call datetime(create_date,create_time) h%create_date = trim(create_date)//' '//trim(create_time) istat = nf_put_att_text(ncid,NF_GLOBAL,"create_date", + len_trim(h%create_date),trim(h%create_date)) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_text ', + 'for create_date global attribute')") call handle_ncerr(istat,char120) endif ! ! User login name: istat = nf_put_att_text(ncid,NF_GLOBAL,"logname", + len_trim(h%logname),trim(h%logname)) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_text ', + 'for logname global attribute')") call handle_ncerr(istat,char120) endif ! ! Machine host name: istat = nf_put_att_text(ncid,NF_GLOBAL,"host", + len_trim(h%host),trim(h%host)) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_text ', + 'for host global attribute')") call handle_ncerr(istat,char120) endif ! ! Machine operating system (from pre-proc macro) istat = nf_put_att_text(ncid,NF_GLOBAL,"system", + len_trim(h%system),trim(h%system)) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_text ', + 'for system global attribute')") call handle_ncerr(istat,char120) endif ! ! Model name string: istat = nf_put_att_text(ncid,NF_GLOBAL,"model_name", + len_trim(tgcm_name),trim(tgcm_name)) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_text ', + 'for model name global attribute')") call handle_ncerr(istat,char120) endif ! ! Model version string: istat = nf_put_att_text(ncid,NF_GLOBAL,"model_version", + len_trim(tgcm_version),trim(tgcm_version)) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_text ', + 'for model version global attribute')") call handle_ncerr(istat,char120) endif ! ! SVN version (see included file revision.inc, created by the build): istat = nf_put_att_text(ncid,NF_GLOBAL,"svn_version", + len_trim(svn_version),trim(svn_version)) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_text ', + 'for SVN version global attribute')") call handle_ncerr(istat,char120) endif ! ! Output path (from namelist): istat = nf_put_att_text(ncid,NF_GLOBAL,"output_file", | len_trim(h%output_file),trim(h%output_file)) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_text ', | 'for output_file global attribute')") call handle_ncerr(istat,char120) endif ! ! ! Primary or Secondary histories (h%hist_type): istat = nf_put_att_text(ncid,NF_GLOBAL,"history_type", | len_trim(h%hist_type),trim(h%hist_type)) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_text ', | 'for history_type global attribute')") call handle_ncerr(istat,char120) endif ! ! Run type ('initial' or 'continuation') istat = nf_put_att_text(ncid,NF_GLOBAL,"run_type", | len_trim(h%run_type),trim(h%run_type)) if (istat /= NF_NOERR) then write(char120,"('nc_define: Error return from nf_put_att_text ', | 'for run_type global attribute')") call handle_ncerr(istat,char120) endif ! ! Source history file (for this run, either initial or continuation): if (nsource > 0) then ! initial source file istat = nf_put_att_text(ncid,NF_GLOBAL,"source_file", | len_trim(h%source_file)+10,trim(h%source_file)//' (initial)') else ! continuation source file istat = nf_put_att_text(ncid,NF_GLOBAL,"source_file", | len_trim(h%source_file)+15, | trim(h%source_file)//' (continuation)') endif if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_text ', | 'for source_file global attribute')") call handle_ncerr(istat,char120) endif ! ! Source model time (will be same as start time for this run, ! whether initial or continuation) istat = nf_put_att_int(ncid,NF_GLOBAL,"source_mtime",NF_INT,3, | h%source_mtime) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_text ', | 'for source_mtime global attribute')") call handle_ncerr(istat,char120) endif ! ! Initial file (source file from the initial run) istat = nf_put_att_text(ncid,NF_GLOBAL,"initial_file", | len_trim(h%initial_file),trim(h%initial_file)) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_text ', | 'for initial_file global attribute')") call handle_ncerr(istat,char120) endif ! ! Initial mtime (this is also an attribute of the time coord var) istat = nf_put_att_int(ncid,NF_GLOBAL,"initial_mtime",NF_INT,3, | h%initial_mtime) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_text ', | 'for initial_mtime global attribute')") call handle_ncerr(istat,char120) endif ! ! Method 1 for conversion from lev (ln(p0/p)) to pressure in hPa ! (millibars), using p0: ! write(char80,"('p0*exp(-lev(k))')") istat = nf_put_att_text(ncid,NF_GLOBAL,"lev_to_hPa_method1", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_text ', | 'for lev_to_hPa_method1 global attribute')") call handle_ncerr(istat,char120) endif ! ! Method 2 for conversion from lev (ln(p0/p)) to pressure in hPa ! (millibars), using p0_model: ! write(char80,"('p0_model*1.e-3*exp(-lev(k))')") istat = nf_put_att_text(ncid,NF_GLOBAL,"lev_to_hPa_method2", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_text ', | 'for lev_to_hPa_method2 global attribute')") call handle_ncerr(istat,char120) endif ! ! Number of histories on the file: ivar1(1) = h%ihist ! this will be updated in nc_wrhist istat = nf_put_att_int(ncid,NF_GLOBAL,"nhist",NF_INT,1,ivar1) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_int ', | 'for nhist global attribute')") call handle_ncerr(istat,char120) endif ! ! Delta minutes between histories on the file: ivar1(1) = h%delhmins istat = nf_put_att_int(ncid,NF_GLOBAL,"delhist_mins",NF_INT, | 1,ivar1) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_int ', | 'for delhist_mins global attribute')") call handle_ncerr(istat,char120) endif ! ! Missing value (nf_put_att_double is in util.F): var1(1) = spval istat = nf_put_att_double(ncid,NF_GLOBAL,"missing_value", | NF_DOUBLE,1,var1) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_double ', | 'for missing_value global attribute')") call handle_ncerr(istat,char120) endif ! ! Potential model used (weimer, heelis, or none): istat = nf_put_att_text(ncid,NF_GLOBAL,"potential_model", | len_trim(potential_model),trim(potential_model)) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_text ', | 'for potential_model global attribute')") call handle_ncerr(istat,char120) endif ! ! tuv_lbc_intop==0 (lbc of t,u,v are not stored in top k-slot) ! (only "old" histories stored lbc of t,u,v in top slot) ivar1(1) = h%tuv_lbc_intop ! see init in hist.F istat = nf_put_att_int(ncid,NF_GLOBAL,"tuv_lbc_intop", | NF_INT,1,ivar1) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_int ', | 'for tuv_lbc_intop global attribute')") call handle_ncerr(istat,char120) endif ! ! wacxg5_nlev is number of levels to be nudged by wacxg5: ! (wacxg5_nlev==1 means only 100% at lower boundary, no nudging) ivar1(1) = wacxg5_nlev istat=nf_put_att_int(ncid,NF_GLOBAL,"wacxg5_nlev",NF_INT,1,ivar1) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_int ', | 'for global attribute wacxg5_nlev')") call handle_ncerr(istat,char120) endif ! ! Exit define mode: istat = nf_enddef(ncid) if (istat /= NF_NOERR) call handle_ncerr(istat, + 'Error return from nf_enddef') ! ! Give values to dimension variables: ! ! Geographic horizontal: istat = nf_put_var_double(ncid,idv_lon,glon) if (istat /= NF_NOERR) | call handle_ncerr(istat,'Error giving values to glon coord var') istat = nf_put_var_double(ncid,idv_lat,glat) if (istat /= NF_NOERR) | call handle_ncerr(istat,'Error giving values to glat coord var') ! ! Vertical: istat = nf_put_var_double(ncid,idv_lev,zpmid) ! midpoint levels if (istat /= NF_NOERR) | call handle_ncerr(istat,'Error giving values to zpmid coord var') istat = nf_put_var_double(ncid,idv_ilev,zpint) ! interface levels if (istat /= NF_NOERR) | call handle_ncerr(istat,'Error giving values to zpint coord var') ! ! Bottom interface level LBC for t,u,v lbc: istat = nf_put_var_double(ncid,idv_lbc,zibot) if (istat /= NF_NOERR) | call handle_ncerr(istat,'Error giving value to LBC scalar var') ! ! Magnetic: istat = nf_put_var_double(ncid,idv_mlon,gmlon) istat = nf_put_var_double(ncid,idv_mlat,gmlat) istat = nf_put_var_double(ncid,idv_mlev,zpmid) istat = nf_put_var_double(ncid,idv_imlev,zpint) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_var_double ', + ' to assign values to magnetic dimension vars')") call handle_ncerr(istat,char120) endif ! end subroutine nc_define !----------------------------------------------------------------------- subroutine define_ncfile(ncid,name,idv,att1name,att1val, | att2name,att2val) ! ! Args: integer,intent(in) :: ncid integer,intent(out) :: idv character(len=*),intent(in) :: name,att1name,att1val, | att2name,att2val ! ! Local: integer :: istat,ids2(2) character(len=120) :: char120 ! if (len_trim(name) <= 0) then write(6,"('>>> define_ncfile: empty name! returning..')") return endif ! ids2(1) = id_filelen ids2(2) = id_time ! ! Define time-dependent char variable: istat= nf_def_var(ncid,name,NF_CHAR,2,ids2,idv) if (istat /= NF_NOERR) then write(char120,"('Error defining file var ',a)") trim(name) call handle_ncerr(istat,char120) endif ! ! First attribute: if (len_trim(att1name) > 0.and.len_trim(att1val) > 0) then istat = nf_put_att_text(ncid,idv,att1name,len_trim(att1val), | trim(att1val)) if (istat /= NF_NOERR) then write(char120,"('Error defining ',a,' attribute to var ', | a,': attribute value=',a)") trim(name),trim(att1name), | trim(att1val) call handle_ncerr(istat,char120) endif endif ! ! Second attribute: if (len_trim(att2name) > 0.and.len_trim(att2val) > 0) then istat = nf_put_att_text(ncid,idv,att2name,len_trim(att2val), | trim(att2val)) if (istat /= NF_NOERR) then write(char120,"('Error defining ',a,' attribute to var ', | a,': attribute value=',a)") trim(name),trim(att2name), | trim(att2val) call handle_ncerr(istat,char120) endif endif end subroutine define_ncfile !------------------------------------------------------------------- subroutine defvar_time_dbl(ncid,name,idtime,idvar, | longname,units) ! call defvar_time_dbl(ncid,"ut",ids1,idv_ut, ! | "universal time (from model time hour and minute)", ! | "hours") ! ! Define time-dependent variable "name" on open netcdf file ncid. ! Variable has a single dimension in time (e.g., ut(time)). ! Also define long_name and units attributes. ! ! Args: integer,intent(in) :: ncid,idtime(1) integer,intent(out) :: idvar character(len=*),intent(in) :: name,longname,units ! ! Local: integer :: istat character(len=80) :: char80 ! ! Define the variable: istat = nf_def_var(ncid,name,NF_DOUBLE,1,idtime,idvar) if (istat /= NF_NOERR) then write(char80,"('Error defining double variable ',a)") name call handle_ncerr(istat,trim(char80)) endif ! ! Give long name attribute: if (len_trim(longname) > 0) then istat = nf_put_att_text(ncid,idvar,"long_name", | len_trim(longname),longname) if (istat /= NF_NOERR) then write(char80,"('Error defining long_name of double variable ', | a)") name call handle_ncerr(istat,trim(char80)) endif endif ! ! Give units attribute: if (len_trim(units) > 0) then istat = nf_put_att_text(ncid,idvar,"units",len_trim(units), | units) if (istat /= NF_NOERR) then write(char80,"('Error defining units of variable ',a)") name call handle_ncerr(istat,trim(char80)) endif endif end subroutine defvar_time_dbl !----------------------------------------------------------------------- subroutine deflbc(ncid,vname,idv,type) use input_module,only: sech_nbyte ! ! Define t, u, or v lbc on current history file. Return new variable ! id in idv. (vname should be TLBC, ULBC, or VLBC) ! If defining a secondary history, check sech_nbyte. ! ! Args: integer,intent(in) :: ncid integer,intent(out) :: idv character(len=*),intent(in) :: vname,type ! ! Local: integer :: istat,len character(len=80) :: char80,units ! ids3(1) = id_lon ids3(2) = id_lat ids3(3) = id_time ! ! Define variable: ! Always write doubles to primary files. For secondary files, write ! 4-byte floats if sech_nbyte=4, or 8-byte doubles if sech_nbyte=8 ! if (trim(type)=='primary') then istat = nf_def_var(ncid,trim(vname),NF_DOUBLE,3,ids3,idv) else if (sech_nbyte==4) then istat = nf_def_var(ncid,trim(vname),NF_REAL,3,ids3,idv) else istat = nf_def_var(ncid,trim(vname),NF_DOUBLE,3,ids3,idv) endif endif if (istat /= NF_NOERR) then write(char80,"('Error defining 2-d var ',a)") trim(vname) call handle_ncerr(istat,trim(char80)) endif ! ! Long name attribute: ! (Variable names ending in '_NM' are at previous time n-1) ! write(char80,"('Lower boundary condition of ',a,'N')") vname(1:1) len = len_trim(vname) if (vname(len-2:len)=='_NM') | char80 = trim(char80)//' (TIME N-1)' ! istat = nf_put_att_text(ncid,idv,"long_name", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) then write(char80,"('Error defining long_name attribute of ', | 'variable ',a)") trim(vname) call handle_ncerr(istat,trim(char80)) endif ! ! Units attribute: units = 'cm/s' if (vname(1:1)=='T') units = 'K' istat = nf_put_att_text(ncid,idv,"units",len_trim(units), | trim(units)) if (istat /= NF_NOERR) then write(char80,"('Error defining units attribute of variable ', | a)") trim(vname) call handle_ncerr(istat,trim(char80)) endif ! ! Coordinates attribute: istat = nf_put_att_text(ncid,idv,"coordinates",3,'LBC') if (istat /= NF_NOERR) then write(char80,"('Error defining coordinates attribute of ', | 'variable',a)") trim(vname) call handle_ncerr(istat,trim(char80)) endif ! write(6,"('deflbc: defined var ',a,' idv=',i3)") trim(vname),idv end subroutine deflbc !----------------------------------------------------------------------- subroutine def_fsech(ncid) ! ! Define secondary history fields (geo or mag, 2d+time or 3d+time): ! use fields_module,only: fsechist use hist_module,only: nsech,isecout use input_module,only: sech_nbyte implicit none ! ! Args: integer,intent(in) :: ncid ! ! Local: integer :: ix,istat,iprog,idimids(4),idv,isdef character(len=80) :: char80 real :: var1(1) ! ! External: integer,external :: strloc ! ! write(6,"('Enter def_fsech: isecout=',i3,' nsech=',i3,' istep=', ! | i4)") isecout,nsech,istep if (nsech==2) then ! write(6,"('def_fsech: call nf_redef')") istat = nf_redef(ncid) ! put dataset in define mode if (istat /= NF_NOERR) then write(char80,"('def_fsech: error from nf_redef')") call handle_ncerr(istat,trim(char80)) endif endif ! ! Fields loop: floop: do ix=1,mxfsech if (len_trim(fsechist(ix)%short_name)==0) cycle floop ! ! Define prognostic only if this is 1st history on this file: iprog=strloc(f4d%short_name,size(f4d),fsechist(ix)%short_name) if (iprog > 0.and.nsech==1) then ! ! Define prognostic secondary history field: call defvar_f4d(ncid,f4d(iprog),idv_fsech(ix),sech_nbyte) ! write(6,"('def_fsech: defined 4d prog field ',a,' ix=', ! | i4)") f4d(ix)%short_name,ix ! ! 10/27/05 btf: Assume user-requested 'W' refers to "old hist" W, i.e. OMEGA: ! elseif (nsech==1.and.fsechist(ix)%short_name=='W') then write(6,"('Note: defining requested secondary history ', | 'variable ''W'' as prognostic ''OMEGA''')") iprog = strloc(f4d%short_name,size(f4d),'OMEGA') call defvar_f4d(ncid,f4d(iprog),idv_fsech(ix),sech_nbyte) cycle floop endif ! iprog and nsech ! ! Define diagostics only if not first timestep (because in that case, ! diags have not yet been allocated by addfld): if (istep==0.or.iprog > 0) cycle floop ! if (nsech==1.or.iprog > 0) cycle floop ! ! If data pointer is not associated, addfld was probably never called ! for this field. In this case, do not define the field. (subs def_fsech ! and wrfsech will also check this pointer, see nchist.F, and sub ! mp_gather2root_fsech in mpi.F) ! if (.not.associated(fsechist(ix)%data)) then write(6,"('>>> WARNING: Field ',a,' will not be defined', | ' on secondary histories because',/,4x,'data was not ', | 'allocated. (addfld was probably not called for this ', | 'field)')") trim(fsechist(ix)%short_name) cycle floop endif ! ! Get dim id's for 3 dimensions (geo or mag): if (fsechist(ix)%dimsizes(3) /= 0) then ! 3d var (lon,lat,lev,time) istat=nf_inq_dimid(ncid,fsechist(ix)%dimnames(1),idimids(1)) if (istat /= NF_NOERR) then write(char80,"('def_fsech: Error getting id for dim ', | a,' (1st of 3)')") trim(fsechist(ix)%dimnames(1)) call handle_ncerr(istat,trim(char80)) endif istat=nf_inq_dimid(ncid,fsechist(ix)%dimnames(2),idimids(2)) if (istat /= NF_NOERR) then write(char80,"('def_fsech: Error getting id for dim ', | a,' (2nd of 3)')") trim(fsechist(ix)%dimnames(2)) call handle_ncerr(istat,trim(char80)) endif istat=nf_inq_dimid(ncid,fsechist(ix)%dimnames(3),idimids(3)) if (istat /= NF_NOERR) then write(char80,"('def_fsech: Error getting id for dim ', | a,' (3rd of 3)')") trim(fsechist(ix)%dimnames(3)) call handle_ncerr(istat,trim(char80)) endif idimids(4) = id_time ! ! Define 3d var (lon,lat,lev,time) istat = nf_inq_varid(ncid,trim(fsechist(ix)%short_name),idv) isdef = 0 ! assume var is not yet defined if (istat == 0) isdef = 1 ! var is already defined -- cannot redef if (isdef == 0) then ! safe to define if (sech_nbyte == 4) then ! write 4-byte reals istat = nf_def_var(ncid,fsechist(ix)%short_name,NF_REAL, | 4,idimids,idv_fsech(ix)) ! returns idv_fsech(ix) else ! write 8-byte doubles istat = nf_def_var(ncid,fsechist(ix)%short_name,NF_DOUBLE, | 4,idimids,idv_fsech(ix)) ! returns idv_fsech(ix) endif if (istat /= NF_NOERR) then write(char80,"('def_fsech: Error defining 3d field ',a)") | trim(fsechist(ix)%short_name) call handle_ncerr(istat,trim(char80)) endif ! write(6,"('def_fsech: defined 3d var ',a)") ! | fsechist(ix)%short_name endif ! if var not yet defined ! ! Get dim id's for 2 dimensions (geo or mag): else ! 2d var (lon,lat,time) istat=nf_inq_dimid(ncid,fsechist(ix)%dimnames(1),idimids(1)) if (istat /= NF_NOERR) then write(char80,"('def_fsech: Error getting id for dim ', | a,' (1st of 2)')") trim(fsechist(ix)%dimnames(1)) call handle_ncerr(istat,trim(char80)) endif istat=nf_inq_dimid(ncid,fsechist(ix)%dimnames(2),idimids(2)) if (istat /= NF_NOERR) then write(char80,"('def_fsech: Error getting id for dim ', | a,' (2nd of 2)')") trim(fsechist(ix)%dimnames(2)) call handle_ncerr(istat,trim(char80)) endif idimids(3) = id_time ! ! Define 2d var (lon,lat,time) if (sech_nbyte == 4) then ! write 4-byte reals istat = nf_def_var(ncid,fsechist(ix)%short_name,NF_REAL, | 3,idimids(1:3),idv_fsech(ix)) ! returns idv_fsech(ix) else istat = nf_def_var(ncid,fsechist(ix)%short_name,NF_DOUBLE, | 3,idimids(1:3),idv_fsech(ix)) ! returns idv_fsech(ix) endif if (istat /= NF_NOERR) then write(char80,"('def_fsech: Error defining 2d field ',a)") | fsechist(ix)%short_name call handle_ncerr(istat,trim(char80)) endif endif ! 2d or 3d ! write(6,"('def_fsech: defined sech field ',a,': dimnames=', ! | a,',',a,',',a,' dimsizes=',3i4)") fsechist(ix)%short_name, ! | trim(fsechist(ix)%dimnames(1)),trim(fsechist(ix)%dimnames(2)), ! | trim(fsechist(ix)%dimnames(3)),fsechist(ix)%dimsizes ! ! Attributes: if (len_trim(fsechist(ix)%long_name) > 0) then istat = nf_put_att_text(ncid,idv_fsech(ix),"long_name", | len_trim(fsechist(ix)%long_name),trim(fsechist(ix)%long_name)) endif if (len_trim(fsechist(ix)%units) > 0) then istat = nf_put_att_text(ncid,idv_fsech(ix),"units", | len_trim(fsechist(ix)%units),trim(fsechist(ix)%units)) endif ! ! Missing value: var1(1) = spval istat = nf_put_att_double(ncid,idv_fsech(ix),"missing_value", | NF_DOUBLE,1,var1) if (istat /= NF_NOERR) then write(char80,"('def_fsech: Error return from ', | 'nf_put_att_double for missing_value var attribute')") call handle_ncerr(istat,char80) endif ! ! End field loop: enddo floop ! i=1,mxfsech if (nsech==2) istat = nf_enddef(ncid) ! take out of define mode end subroutine def_fsech !----------------------------------------------------------------------- subroutine defvar_f4d(ncid,f,idvar,nbytes) ! ! Define a prognostic variable on current netcdf history file: ! (not data -- see wrf4d for data write) ! ! Args: integer,intent(in) :: ncid,nbytes integer,intent(out) :: idvar type(fields_4d),intent(in) :: f ! ! Local: integer :: istat,idimids(4) character(len=80) :: char80 character(len=160) :: char160 real :: var1(1) ! ! Prognostics are dimensioned (lon,lat,lev,time): if (.not.fakeflds) then idimids(1) = id_lon idimids(2) = id_lat idimids(3) = id_lev if (f%vcoord(1:3)=='int') idimids(3) = id_ilev ! interfaces idimids(4) = id_time else idimids(1:3) = id_fake idimids(4) = id_time endif ! ! This call returns the variable id, idvar. ! nbytes == 4 means write 4-byte reals (to secondary files), otherwise ! nbytes == 8 means write 8-byte doubles (either secondary or primary files). ! if (nbytes == 4) then ! write 4-byte reals istat = nf_def_var(ncid,f%short_name,NF_REAL,4,idimids,idvar) else ! write 8-byte doubles istat = nf_def_var(ncid,f%short_name,NF_DOUBLE,4,idimids,idvar) endif if (istat /= NF_NOERR) then write(char80,"('defvar_f4d: Error defining 4-d field ',a, | ' nbytes=',i2)") f%short_name,nbytes call handle_ncerr(istat,trim(char80)) endif ! ! Add long_name attribute: istat = nf_put_att_text(ncid,idvar,"long_name", | len_trim(f%long_name),trim(f%long_name)) ! ! Cray apparently crashed here after 1380 steps ! (char80 was used) if (istat /= NF_NOERR) then write(char160,"('Error defining long_name of 4-d', | ' variable ',a,': long_name = ',/,a)") trim(f%short_name), | trim(f%long_name) call handle_ncerr(istat,trim(char160)) endif ! ! Add units attribute: istat = nf_put_att_text(ncid,idvar,"units",len_trim(f%units), | f%units) if (istat /= NF_NOERR) then write(char80,"('Error defining units of 4-d', | ' variable ',a,': units = ',a)") trim(f%short_name), | trim(f%units) call handle_ncerr(istat,trim(char80)) endif ! ! 10/27/05 btf: Add former_name attribute to OMEGA to indicate it ! was previously known as W: if (f%short_name=='OMEGA') then istat = nf_put_att_text(ncid,idvar,"former_name",1,"W") endif ! ! Add missing_value attribute: var1(1) = spval istat = nf_put_att_double(ncid,idvar,"missing_value", | NF_DOUBLE,1,var1) if (istat /= NF_NOERR) then write(char160,"('defvar_f4d: Error return from ', | 'nf_put_att_double for missing_value var attribute')") call handle_ncerr(istat,char160) endif end subroutine defvar_f4d !------------------------------------------------------------------- subroutine nc_wrhist(ncid,diskfile,iprint) ! ! Write current history structure to open netcdf output file: ! use hist_module,only: sh,h,hist_print use input_module,only: start,secstart ! ! Args: integer,intent(in) :: ncid,iprint character(len=*),intent(in) :: diskfile ! ! Local: integer :: mins,istat,ivar1(1),i,j,iprog,imo,ida,mtimeinit(4) character(len=120) :: char120 character(len=240) :: char240 character(len=80) :: char80 integer :: ndims,nvars,ngatts,unlimdimid real :: var1(1),var22(2,2),var2(2),var10(10),rmins,rdays ! ! External: integer,external :: strloc real,external :: mtime_delta ! ! Time: decimal days since start time of initial run: ! h%initial_mtime was set by sub hist_initype (if initial run), ! or was read from source history (if continuation run). ! mtimeinit(1:3) = sh%initial_mtime(1:3) ; mtimeinit(4) = 0 rmins = mtime_delta(mtimeinit,h%modeltime) rdays = rmins/1440. if (rmins < 0.) then write(6,"(/,'>>> nc_wrhist: time from initial start must be', | ' >= 0: mtimeinit=',3i4,' h%modeltime=',3i4,' delta (mins)=', | f10.2)") mtimeinit(1:3),h%modeltime(1:3),rmins call shutdown('delta time from init run') else ! write(6,"('nc_wrhist: mtimeinit=',3i4,' h%modeltime=',3i4, ! | ' delta (mins)=',f10.2,' delta (days)=',f8.2)") ! | mtimeinit(1:3),h%modeltime(1:3),rmins,rdays endif if (idv_time <= 0) istat = nf_inq_varid(ncid,"time",idv_time) istat = nf_put_var1_double(ncid,idv_time,h%ihist,rdays) if (istat /= NF_NOERR) then write(char240,"('Error from nf_put_var1_double defining ', | 'time at h%modeltime=',4i4,' mtimeinit=',3i4, | ' delta (days)=',f8.2)") h%modeltime,rdays,mtimeinit(1:3) call handle_ncerr(istat,char240) endif ! ! Model time (day,hour,minute): start_2d(1) = 1 start_2d(2) = h%ihist count_2d(1) = 3 count_2d(2) = 1 if (idv_mtime <= 0) istat = nf_inq_varid(ncid,"mtime",idv_mtime) ! ! If model day is > 365, write actual day of year, i.e. day-365 ! (e.g., if model day is 366, write model day 1 to the history) ! if (h%modeltime(1) > 365) then write(6,"('Note: Changed model day from ',i4,' to ',i4, | ' before writing to history file.')") | h%modeltime(1), h%modeltime(1)-365 h%modeltime(1) = h%modeltime(1)-365 endif istat = nf_put_vara_int(ncid,idv_mtime,start_2d,count_2d, | h%modeltime(1:3)) if (istat /= NF_NOERR) then write(char120,"('Error defining mtime at h%ihist=',i2, + ' h%modeltime=',4i4,' istat=',i5)") h%ihist,h%modeltime,istat call handle_ncerr(istat,char120) endif ! write(6,"('nc_wrhist: h%ihist=',i4,' h%modeltime=',4i4)") ! | h%ihist,h%modeltime ! ! Calendar year (initially from input, then updated if model is ! advancing calendar time): if (idv_year <= 0) istat = nf_inq_varid(ncid,"year",idv_year) istat = nf_put_var1_int(ncid,idv_year,h%ihist,h%year) if (istat /= NF_NOERR) then write(char120,"('Error from nf_put_var1_int defining ', | 'calendar year at h%ihist=',i3,': h%year=',i5)") | h%ihist,h%year call handle_ncerr(istat,char120) endif ! ! Calendar day (initially from input, then updated if model is ! advancing calendar time): if (idv_day <= 0) istat = nf_inq_varid(ncid,"day",idv_day) istat = nf_put_var1_int(ncid,idv_day,h%ihist,h%day) if (istat /= NF_NOERR) then write(char120,"('Error from nf_put_var1_int defining ', | 'calendar day at h%ihist=',i3,': h%day=',i5)") | h%ihist,h%day call handle_ncerr(istat,char120) endif ! ! Calendar advance: if (idv_calendar_adv <= 0) | istat = nf_inq_varid(ncid,"calendar_advance",idv_calendar_adv) istat = nf_put_var1_int(ncid,idv_calendar_adv,h%ihist, | h%calendar_advance) if (istat /= NF_NOERR) then write(char120,"('Error from nf_put_var1_int defining ', | 'calendar advance flag at h%ihist=',i3, | ': h%calendar_advance=',i5)") h%ihist,h%calendar_advance call handle_ncerr(istat,char120) endif ! ! Date and time history was written: if (idv_writedate <= 0) istat = nf_inq_varid(ncid,"write_date", | idv_writedate) ! write(6,"('nc_wrhist: h%ihist=',i3,' len_trim(h%write_date)=', ! | i3,' h%write_date=$',a,'$')") h%ihist,len_trim(h%write_date), ! | h%write_date ! if (h%write_date == '\0') ! | write(6,"('h%write_date is equal to \0')") if (len_trim(h%write_date)==0.or.h%write_date == '\0') | h%write_date = '[unknown]' start_2d(1) = 1 start_2d(2) = h%ihist count_2d(1) = len(h%write_date) count_2d(2) = 1 istat = nf_put_vara_text(ncid,idv_writedate,start_2d,count_2d, | h%write_date) if (istat /= NF_NOERR) then write(char120,"('Error from nf_put_vara_text writing ', | 'h%write_date=',a)") trim(h%write_date) call handle_ncerr(istat,char120) endif ! ! Iteration number (number of steps from 0,0,0 to current mtime): if (idv_iter <= 0) istat = nf_inq_varid(ncid,"iter",idv_iter) istat = nf_put_var1_int(ncid,idv_iter,h%ihist,h%iter) if (istat /= NF_NOERR) then write(char120,"('Error from nf_put_var1_int defining ', + 'iter at h%ihist=',i3,': h%iter=',i8)") h%ihist,h%iter call handle_ncerr(istat,char120) endif ! ! Universal time (decimal hours): if (idv_ut <= 0) istat = nf_inq_varid(ncid,"ut",idv_ut) istat = nf_put_var1_double(ncid,idv_ut,h%ihist,h%ut) ! ! 1/31/01: tides are now time dependent 2d vars: ! Diurnal tide (dtide(time,2): start_2d(1) = 1 start_2d(2) = h%ihist count_2d(1) = 2 count_2d(2) = 1 if (idv_dtide <= 0) istat = nf_inq_varid(ncid,"dtide",idv_dtide) var2(:) = h%dtide(:) istat = nf_put_vara_double(ncid,idv_dtide,start_2d,count_2d,var2) ! ! Semi-diurnal tide (sdtide(time,10): start_2d(1) = 1 start_2d(2) = h%ihist count_2d(1) = 10 count_2d(2) = 1 if (idv_sdtide <= 0) istat = nf_inq_varid(ncid,"sdtide", | idv_sdtide) var10(:) = h%sdtide(:) istat = nf_put_vara_double(ncid,idv_sdtide,start_2d,count_2d, | var10) ! ! Update gpi_ncfile: call update_ncfile(ncid,"gpi_ncfile",h%gpi_ncfile,idv_gpi_ncfile, | h%ihist) ! ! Update T,Z lbc data files: call update_ncfile(ncid,"ncep_ncfile",h%ncep_ncfile, | idv_ncep_ncfile,h%ihist) call update_ncfile(ncid,"ncep_reanalysis", | h%ncep_reanalysis,idv_ncep_reanalysis,h%ihist) call update_ncfile(ncid,"ecmwf_ncfile", | h%ecmwf_ncfile,idv_ecmwf_ncfile,h%ihist) call update_ncfile(ncid,"wacxg5_ncfile", | h%wacxg5_ncfile,idv_wacxg5_ncfile,h%ihist) ! ! Update gswm files: call update_ncfile(ncid,"gswm_mi_di_ncfile",h%gswm_mi_di_ncfile, | idv_gswm_mi_di_ncfile,h%ihist) call update_ncfile(ncid,"gswm_mi_sdi_ncfile",h%gswm_mi_sdi_ncfile, | idv_gswm_mi_sdi_ncfile,h%ihist) call update_ncfile(ncid,"gswm_nm_di_ncfile",h%gswm_nm_di_ncfile, | idv_gswm_nm_di_ncfile,h%ihist) call update_ncfile(ncid,"gswm_nm_sdi_ncfile",h%gswm_nm_sdi_ncfile, | idv_gswm_nm_sdi_ncfile,h%ihist) ! ! f107d and f107a (daily and average solar fluxes): if (idv_f107d <= 0) istat = nf_inq_varid(ncid,"f107d", | idv_f107d) istat = nf_put_var1_double(ncid,idv_f107d,h%ihist,h%f107d) if (idv_f107a <= 0) istat = nf_inq_varid(ncid,"f107a", | idv_f107a) istat = nf_put_var1_double(ncid,idv_f107a,h%ihist,h%f107a) ! ! Hemispheric power: if (idv_hpower <= 0) istat = nf_inq_varid(ncid,"hpower", | idv_hpower) istat = nf_put_var1_double(ncid,idv_hpower,h%ihist,h%hpower) ! ! Cross-tail potential: if (idv_ctpoten <= 0) istat = nf_inq_varid(ncid,"ctpoten", | idv_ctpoten) istat = nf_put_var1_double(ncid,idv_ctpoten,h%ihist,h%ctpoten) ! ! BY imf: if (idv_byimf <= 0) istat = nf_inq_varid(ncid,"byimf",idv_byimf) istat = nf_put_var1_double(ncid,idv_byimf,h%ihist,h%byimf) ! ! Collision factor: if (idv_colfac <= 0) istat = nf_inq_varid(ncid,"colfac", | idv_colfac) istat = nf_put_var1_double(ncid,idv_colfac,h%ihist,h%colfac) ! ! alfa30: if (idv_alfa30 <= 0) istat = nf_inq_varid(ncid,"alfa30", | idv_alfa30) istat = nf_put_var1_double(ncid,idv_alfa30,h%ihist,h%alfa30) ! ! e30: if (idv_e30 <= 0) istat = nf_inq_varid(ncid,"e30",idv_e30) istat = nf_put_var1_double(ncid,idv_e30,h%ihist,h%e30) ! ! alfad2: if (idv_alfad2 <= 0) istat = nf_inq_varid(ncid,"alfad2", | idv_alfad2) istat = nf_put_var1_double(ncid,idv_alfad2,h%ihist,h%alfad2) ! ! ed2: if (idv_ed2 <= 0) istat = nf_inq_varid(ncid,"ed2",idv_ed2) istat = nf_put_var1_double(ncid,idv_ed2,h%ihist,h%ed2) ! ! Reference pressure used to convert ln(p0/p) to pressure (hPa): ! (h%p0 was set in sub define_hist, output.F) if (idv_p0 <= 0) istat = nf_inq_varid(ncid,"p0",idv_p0) istat = nf_put_var1_double(ncid,idv_p0,h%ihist,h%p0) ! ! Reference pressure used by the model (microbars): ! (h%p0_model was set in sub define_hist, output.F) if (idv_p0_model <= 0) istat = nf_inq_varid(ncid,"p0_model", | idv_p0_model) istat = nf_put_var1_double(ncid,idv_p0_model,h%ihist,h%p0_model) ! ! Gravitational constant used in the model: if (idv_grav <= 0) istat = nf_inq_varid(ncid,"grav",idv_grav) istat = nf_put_var1_double(ncid,idv_grav,h%ihist,h%grav) ! ! Timestep (seconds): if (idv_step <= 0) istat = nf_inq_varid(ncid,"timestep",idv_step) istat = nf_put_var1_int(ncid,idv_step,h%ihist,h%step) if (istat /= NF_NOERR) then write(char120,"('Error from nf_put_var1_int defining ', + 'ut at h%ihist=',i3,': h%step=',i4)") h%ihist,h%step call handle_ncerr(istat,char120) endif ! ! Update run_type (initial or continuation): istat = nf_redef(ncid) ! put dataset in define mode istat = nf_put_att_text(ncid,NF_GLOBAL,"run_type", | len_trim(h%run_type),trim(h%run_type)) if (istat /= NF_NOERR) then write(char120,"('nc_wrhist: Error return from nf_put_att_text ', | 'for run_type global attribute')") call handle_ncerr(istat,char120) endif ! ! Update source_file and source mtime: if (h%run_type(1:4)=='init') then ! initial source file istat = nf_put_att_text(ncid,NF_GLOBAL,"source_file", | len_trim(h%source_file)+10,trim(h%source_file)//' (initial)') elseif (h%run_type(1:4)=='cont') then ! continuation source file istat = nf_put_att_text(ncid,NF_GLOBAL,"source_file", | len_trim(h%source_file)+15, | trim(h%source_file)//' (continuation)') endif if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_text ', | 'for source_file global attribute')") call handle_ncerr(istat,char120) endif istat = nf_put_att_int(ncid,NF_GLOBAL,"source_mtime",NF_INT,3, | h%source_mtime) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_text ', | 'for source_mtime global attribute')") call handle_ncerr(istat,char120) endif istat = nf_enddef(ncid) ! take out of define mode ! ! Update number of histories on the file (global attribute): ivar1(1) = h%ihist istat = nf_put_att_int(ncid,NF_GLOBAL,"nhist",NF_INT,1,ivar1) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_int ', + 'for nhist global attribute')") call handle_ncerr(istat,char120) endif ! ! Add fields data to primary history: ! Currently (2/00), only prognostic fields are written to the ! primary histories: ! idv_tlbc = 0 ; idv_ulbc = 0 ; idv_vlbc = 0 ! forces wrlbc to inquire idv_tlbc_nm = 0 ; idv_ulbc_nm = 0 ; idv_vlbc_nm = 0 ! forces wrlbc to inquire if (h%hist_type(1:3)=='pri') then do i=1,nf4d_hist call wrf4d(ncid,f4d(i)%short_name,h%ihist,fakeflds,i) enddo ! ! Update t,u,v lbc to primary history: ! Global arrays t,u,vlbc_glb were gathered from subdomains for history ! output by sub mp_gather2root_lbc (mpi.F). ! call wrlbc(ncid,"TLBC",tlbc_glb,h%ihist,idv_tlbc,'pri') call wrlbc(ncid,"ULBC",ulbc_glb,h%ihist,idv_ulbc,'pri') call wrlbc(ncid,"VLBC",vlbc_glb,h%ihist,idv_vlbc,'pri') ! call wrlbc(ncid,"TLBC_NM",tlbc_nm_glb,h%ihist,idv_tlbc_nm,'pri') call wrlbc(ncid,"ULBC_NM",ulbc_nm_glb,h%ihist,idv_ulbc_nm,'pri') call wrlbc(ncid,"VLBC_NM",vlbc_nm_glb,h%ihist,idv_vlbc_nm,'pri') ! ! Add fields data to secondary history: ! elseif (h%hist_type(1:3)=='sec') then call wrfsech(ncid,h%ihist) ! write secondary prog and diag fields ! ! Update t,u,v lbc to secondary history: call wrlbc(ncid,"TLBC",tlbc_glb,h%ihist,idv_tlbc,'sec') call wrlbc(ncid,"ULBC",ulbc_glb,h%ihist,idv_ulbc,'sec') call wrlbc(ncid,"VLBC",vlbc_glb,h%ihist,idv_vlbc,'sec') call wrlbc(ncid,"TLBC_NM",tlbc_nm_glb,h%ihist,idv_tlbc_nm,'sec') call wrlbc(ncid,"ULBC_NM",ulbc_nm_glb,h%ihist,idv_ulbc_nm,'sec') call wrlbc(ncid,"VLBC_NM",vlbc_nm_glb,h%ihist,idv_vlbc_nm,'sec') else write(6,"(/,'>>> WARNING nc_wrhist: unknown h%hist_type=',a)") | h%hist_type endif ! ! Report to stdout: if (iprint > 0) call hist_print(h,'WRITE',diskfile) end subroutine nc_wrhist !----------------------------------------------------------------------- subroutine wrfsech(ncid,itime) ! ! Write data values of secondary history fields to ncfile. ! (fields were defined on the file by def_fsech) ! Data was gathered to root task by mp_gather2root_sechist (nchist.F) ! use fields_module,only: fsechist,shortname_len use hist_module,only: nsech,h implicit none ! ! Args: integer,intent(in) :: ncid,itime ! ! Local: integer :: ix,iprog,istat,idvar,iddims(4),itype,ndims,natts,i character(len=80) :: char80 real,allocatable :: f2d(:,:),f3d(:,:,:) ! ! External: integer,external :: strloc ! ! write(6,"('Enter wrfsech: nsech=',i3,' itime=',i3,' istep=', ! | i4)") nsech,itime,istep floop: do ix=1,mxfsech if (len_trim(fsechist(ix)%short_name)==0) cycle floop ! ! Write prognostic: iprog=strloc(f4d%short_name,size(f4d),fsechist(ix)%short_name) if (iprog > 0) then call wrf4d(ncid,fsechist(ix)%short_name,h%ihist,fakeflds, | iprog) cycle floop elseif (fsechist(ix)%short_name=='W') then write(6,"('Note: defining requested secondary history ', | 'variable ''W'' as prognostic ''OMEGA''')") iprog = strloc(f4d%short_name,size(f4d),'OMEGA') call wrf4d(ncid,'OMEGA',h%ihist,fakeflds,iprog) cycle floop endif ! ! Wrote diagnostic only if not first timestep: if (istep==0) cycle floop ! ! If addfld was never called for this field, data will not be ! allocated -- in this case, leave the field defined on the ! history, but do not write data (netcdf will default to _Fillvalue) ! if (.not.associated(fsechist(ix)%data)) then write(6,"('>>> WARNING: sechist field ',a,': data not', | ' allocated.')") trim(fsechist(ix)%short_name) cycle floop endif ! ! Get field id: istat = nf_inq_varid(ncid,fsechist(ix)%short_name,idvar) if (istat /= NF_NOERR) then write(char80,"('wrfsech: Error getting var id of field ',a)") | trim(fsechist(ix)%short_name) call handle_ncerr(istat,char80) endif ! write(6,"('wrfsech write diag ',a,' idvar=',i3, ! | ' fsechist(ix)%dimsizes=',4i4,' fsechist(ix)%geo=',l1)") ! | fsechist(ix)%short_name,idvar,fsechist(ix)%dimsizes, ! | fsechist(ix)%geo ! ! 3d+time field: if (fsechist(ix)%dimsizes(3) > 0) then ! ! Allocate and transfer to local 3d output array: allocate(f3d(fsechist(ix)%dimsizes(1), | fsechist(ix)%dimsizes(2), | fsechist(ix)%dimsizes(3)),stat=istat) if (istat /= 0) then write(6,"(/,'>>> wrfsech: Error allocating f3d: ', | 'dimsizes=',3i4)") fsechist(ix)%dimsizes call shutdown('wrfsech') endif f3d = fsechist(ix)%data ! ! Write values for 3d+time var (lon,lat,lev,time) start_4d(1:3) = 1 start_4d(4) = itime count_4d(2:3) = fsechist(ix)%dimsizes(2:3) count_4d(4) = 1 ! ! Write 3d+time geo field: ! (write only 3:nlonp2 to exclude periodic points): if (fsechist(ix)%geo) then ! 3d geo count_4d(1) = nlon ! write(6,"('wrfsech: writing 3d sech field ',a,'(',a,'=',i3, ! | ',',a,'=',i3,',',a,'=',i3,'): min,max=',2e12.4)") ! | fsechist(ix)%short_name(1:10), ! | trim(fsechist(ix)%dimnames(1)),fsechist(ix)%dimsizes(1), ! | trim(fsechist(ix)%dimnames(2)),fsechist(ix)%dimsizes(2), ! | trim(fsechist(ix)%dimnames(3)),fsechist(ix)%dimsizes(3), ! | minval(f3d(3:nlonp2,:,:)),maxval(f3d(3:nlonp2,:,:)) istat = nf_put_vara_double(ncid,idvar,start_4d,count_4d, | f3d(3:nlonp2,:,:)) ! ! Write 3d+time mag field: else ! 3d mag count_4d(1) = nmlonp1 ! write(6,"('wrfsech: writing 3d mag sech field ',a,'(',a,'=', ! | i3,',',a,'=',i3,',',a,'=',i3,'): min,max=',2e12.4)") ! | trim(fsechist(ix)%short_name), ! | trim(fsechist(ix)%dimnames(1)),fsechist(ix)%dimsizes(1), ! | trim(fsechist(ix)%dimnames(2)),fsechist(ix)%dimsizes(2), ! | trim(fsechist(ix)%dimnames(3)),fsechist(ix)%dimsizes(3), ! | minval(f3d),maxval(f3d) istat = nf_put_vara_double(ncid,idvar,start_4d,count_4d, | f3d) endif ! geo or mag if (istat /= NF_NOERR) then write(char80,"('wrfsech: Error writing 3d+time field ',a)") | trim(fsechist(ix)%short_name) call handle_ncerr(istat,char80) endif ! ! Write 2d+time var (lon,lat,time): else ! ! Allocate and transfer to local 2d array: allocate(f2d(fsechist(ix)%dimsizes(1), | fsechist(ix)%dimsizes(2)),stat=istat) if (istat /= 0) then write(6,"(/,'>>> wrfsech: Error allocating f2d: ', | 'dimsizes=',3i4)") fsechist(ix)%dimsizes(1:2) call shutdown('wrfsech') endif f2d = fsechist(ix)%data(:,:,1) ! ! Write values for 2d+time var (lon,lat,time) start_3d(1:2) = 1 start_3d(3) = itime count_3d(2) = fsechist(ix)%dimsizes(2) count_3d(3) = 1 ! ! Write 2d+time geo field: ! (write only 3:nlonp2 to exclude periodic points): if (fsechist(ix)%geo) then ! 2d geo count_3d(1) = nlon ! write(6,"('wrfsech: writing 2d geo sech field ',a,'(',a,'=', ! | i3,',',a,'=',i3,' itime=',i3,'): min,max=',2e12.4)") ! | fsechist(ix)%short_name(1:10), ! | trim(fsechist(ix)%dimnames(1)),fsechist(ix)%dimsizes(1), ! | trim(fsechist(ix)%dimnames(2)),fsechist(ix)%dimsizes(2), ! | itime,minval(f2d(3:nlonp2,:)),maxval(f2d(3:nlonp2,:)) istat = nf_put_vara_double(ncid,idvar,start_3d,count_3d, | f2d(3:nlonp2,:)) ! ! Write 2d+time mag field: else ! 2d mag count_3d(1) = nmlonp1 ! write(6,"('wrfsech: writing 2d mag sech field ',a,'(',a,'=', ! | i3,',',a,'=',i3,'): min,max=',2e12.4)") ! | trim(fsechist(ix)%short_name), ! | trim(fsechist(ix)%dimnames(1)),fsechist(ix)%dimsizes(1), ! | trim(fsechist(ix)%dimnames(2)),fsechist(ix)%dimsizes(2), ! | minval(f2d),maxval(f2d) istat = nf_put_vara_double(ncid,idvar,start_3d,count_3d, | f2d) endif if (istat /= NF_NOERR) then write(char80,"('wrfsech: Error writing 2d+time field ',a)") | trim(fsechist(ix)%short_name) call handle_ncerr(istat,char80) endif endif if (allocated(f2d)) deallocate(f2d) if (allocated(f3d)) deallocate(f3d) enddo floop ! ix=1,mxfsech end subroutine wrfsech !------------------------------------------------------------------- subroutine wrf4d(ncid,name,itime,fake,ix) ! ! Write prognostic field type f%data to current netcdf output ! history file attached to ncid: ! use input_module,only: dynamo use init_module,only: istep use mpi_module,only: lon0,lon1,lat0,lat1 use fields_module,only: dynpot ! ! Args: integer,intent(in) :: ncid,itime,ix character(len=*),intent(in) :: name logical,intent(in) :: fake ! ! Local: integer :: nxk,k,i,j,idvar,istat,lonbeg,lonend,nnans,nanfatal real :: f3d(nlon,nlat,nlevp1) ! note i,j,k for history character(len=120) :: char120 real :: fakevar(1,1,1),fmin,fmax ! ! Fake means 3d fields are dimensioned (1,1,1): ! (see fakeflds in flds_mod.f) if (fake) then istat = nf_inq_varid(ncid,name,idvar) if (istat /= NF_NOERR) then write(char120,"('wrf4d: Error getting id of field var ',a)") | trim(name) call handle_ncerr(istat,char120) endif start_4d(1:3) = 1 start_4d(4) = itime count_4d(1) = 1 count_4d(2) = 1 count_4d(3) = 1 count_4d(4) = 1 fakevar(1,1,1) = 0. istat = nf_put_vara_double(ncid,idvar,start_4d,count_4d,fakevar) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_vara_double', + ' for fake field var ',a,' itime=',i2)") name,itime call handle_ncerr(istat,char120) endif return endif ! ! Transform from (k,i,j) to (i,j,k) for the history. ! This routine is executed by root task only. The root task has ! gathered all subdomain data from other tasks into foutput ! (see call mp_gather2root in advance). Foutput is a pointer to ! all fields at 3d grid, declared in fields.F and allocated in ! allocdata.F). ! ! f3d(nlon,nlat,nlevp1) ! note i,j,k for history ! dynpot(nlonp1,0:nlatp1,nlevp1) ! 3d electric potential geographic ! if (trim(name) /= "POTEN") then do j=1,nlat do k=1,nlevp1 f3d(:,j,k) = foutput(k,3:nlon+2,j,ix) enddo enddo else ! ! Save geographic potential dynpot from dynamo. ! do j=1,nlat do k=1,nlevp1 f3d(:,j,k) = dynpot(1:nlon,j,k) enddo enddo endif ! ! Get field id: istat = nf_inq_varid(ncid,name,idvar) if (istat /= NF_NOERR) then write(char120,"('wrf4d: Error getting id of field var ',a)") | trim(name) call handle_ncerr(istat,char120) endif ! ! Put data onto netcdf file (this is where the majority of output ! i/o happens). ! start_4d(1:3) = 1 start_4d(4) = itime count_4d(1) = nlon count_4d(2) = nlat count_4d(3) = nlevp1 count_4d(4) = 1 istat = nf_put_vara_double(ncid,idvar,start_4d,count_4d,f3d) if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_vara_double', + ' for field var ',a,' itime=',i2)") name,itime call handle_ncerr(istat,char120) endif ! ! Write min,max of each field if not echoing source history: ! if (istep /= 0) then ! call fminmaxspv(f3d,nlon*nlev*nlat,fmin,fmax,spval) ! write(6,"('wrf4d: Wrote field ',a,' istep=',i3, ! | ' 3d min,max=',2e12.4)") name(1:8),istep,fmin,fmax ! endif ! ! Check for NaNs: ! nanfatal = 1 ! if > 0, then abort if nans are found nanfatal = 0 ! if = 0, then do not abort if nans are found if (check_nan) | call check_nans(f3d,nlon,nlat,nlevp1,name(1:8),nnans, | 0,0.,1,nanfatal) end subroutine wrf4d !----------------------------------------------------------------------- subroutine update_ncfile(ncid,vname,ncfile,idv,ihist) ! ! Update file text var vname for history ihist. ! Value of the variable to write is ncfile. ! implicit none ! ! Args: integer,intent(in) :: ncid,ihist integer,intent(inout) :: idv character(len=*),intent(in) :: vname,ncfile ! ! Local: integer :: istat,start(2),count(2) character(len=len(ncfile)) :: file character(len=120) :: char120 ! if (idv <= 0) | istat = nf_inq_varid(ncid,vname,idv) start(1) = 1 start(2) = ihist count(2) = 1 file = ncfile if (len_trim(ncfile) <= 0) file = '[none]' count(1) = len(ncfile) istat = nf_put_vara_text(ncid,idv,start,count,file) if (istat /= NF_NOERR) then write(char120,"('Error from nf_put_vara_text to update ncfile ', | a)") trim(file) call handle_ncerr(istat,char120) else ! write(6,"('Updated var ',a,' for ihist=',i3,' file=',a)") ! | trim(vname),ihist,trim(file) endif end subroutine update_ncfile !------------------------------------------------------------------- subroutine wrlbc(ncid,vname,flbc,itime,idv,type) implicit none ! ! Update global lbc of t,u,v (TLBC,ULBC,VLBC (lon,lat,time)) to ! current history (may be primary or secondary history). ! This is executed by master task only (subdomains were gathered ! into flbc by mp_gather2root_lbc) ! ! Args: integer,intent(in) :: ncid,itime integer,intent(inout) :: idv character(len=*),intent(in) :: vname,type real,intent(in) :: flbc(nlonp4,nlat) ! data to save on file ! ! Local: integer :: istat,j character(len=80) :: char80 character(len=120) :: char120 real :: fglb(nlon,nlat) ! local for write to history ! start_3d(1:2) = 1 start_3d(3) = itime count_3d(1) = nlon count_3d(2) = nlat count_3d(3) = 1 ! ! idv may be 0 for primary continuation file: if (idv <= 0) then istat = nf_inq_varid(ncid,trim(vname),idv) if (istat /= NF_NOERR) then write(char80,"('Error getting var id for ',a, | ' (ncid=',i3,')')") trim(vname),ncid call handle_ncerr(istat,char80) else ! write(6,"('wrlbc: vname=',a,' idv=',i3,' (idv was 0 on ', ! | 'input)')") trim(vname),idv endif else ! write(6,"('wrlbc: vname=',a,' idv on input=',i3)") ! | trim(vname),idv endif ! fglb(:,:) = flbc(3:nlon+2,:) istat = nf_put_vara_double(ncid,idv,start_3d,count_3d,fglb) if (istat /= NF_NOERR) then write(6,"('>>> wrlbc: error writing ',a,': itime=',i4, | ' type=',a,' global (nlon,nlat) min,max=',2f8.2)") | trim(vname),itime,trim(type),minval(fglb),maxval(fglb) write(char80,"('Error return from nf_put_vara_double for ', | 'var ',a,': itime=',i2,' type=',a)") trim(vname),itime, | trim(type) call handle_ncerr(istat,char80) else ! write(6,"('wrlbc: wrote ',a,': itime=',i4,' type=',a, ! | ' global (nlon,nlat) min,max=',2f8.2)") trim(vname), ! | itime,trim(type),minval(fglb),maxval(fglb) endif ! if (trim(vname)=='TLBC') then ! do j=1,nlat ! write(6,"('wrlbc: j=',i3,' tlbc(:,j)=',/,(6e12.4))") ! | j,fglb(:,j) ! enddo ! endif end subroutine wrlbc !----------------------------------------------------------------------- subroutine nc_fileinfo(ncid,fileinfo) ! ! A netcdf file on ncid has been written and is ready to be ! closed and disposed. This routine builds a fileinfo string ! giving beginning and ending history dates, model version string, ! and "primary" or "secondary". This is added as a global ! file attribute, and is included in the mss comment field ! if the file is disposed to the hpss, and is therefore accessible ! via "hsi ls -lA". Example fileinfo string: ! ! 96080 80, 0, 0 to 96084 84, 0, 0 by 720 tgcm22 primary ! ! Args: integer,intent(in) :: ncid character(len=*),intent(out) :: fileinfo ! ! Local: integer :: istat,idvmtime,idtime,ntimes,idvyear,idvday, | mtime_beg(3),mtime_end(3),istart2(2),icount2(2), | iyear_beg,iyear_end,iday_beg,iday_end,idelhist, | ibegiyd,iendiyd character(len=16) :: mversion,htype character(len=80) :: contents_desc ! ! Init: fileinfo = ' ' ! ! Get number of histories (times) on the file: istat = nf_inq_unlimdim(ncid,idtime) ! id of unlimited record var if (istat /= NF_NOERR) call handle_ncerr(istat, | 'nc_fileinfo: Error getting unlimited dimension id') istat = nf_inq_dimlen(ncid,idtime,ntimes) ! length of time var if (istat /= NF_NOERR) call handle_ncerr(istat, | 'nc_fileinfo: Error getting length of time record variable') ! ! Get id for mtime variable: istat = nf_inq_varid(ncid,"mtime",idvmtime) if (istat /= NF_NOERR) | call handle_ncerr(istat, | 'nc_fileinfo: Error getting id of variable mtime') ! ! Get beginning and ending mtimes: istart2(1) = 1 istart2(2) = 1 icount2(1) = 3 icount2(2) = 1 istat = nf_get_vara_int(ncid,idvmtime,istart2,icount2, | mtime_beg) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'nc_fileinfo: Error getting beginning mtime') istart2(2) = ntimes istat = nf_get_vara_int(ncid,idvmtime,istart2,icount2, | mtime_end) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'nc_fileinfo: Error getting ending mtime') ! ! Get beginning and ending year: istat = nf_inq_varid(ncid,"year",idvyear) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'nc_fileinfo: Error getting id of variable year') istat = nf_get_var1_int(ncid,idvyear,1,iyear_beg) istat = nf_get_var1_int(ncid,idvyear,ntimes,iyear_end) ! ! Get beginning and ending day: istat = nf_inq_varid(ncid,"day",idvday) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'nc_fileinfo: Error getting id of variable day') istat = nf_get_var1_int(ncid,idvday,1,iday_beg) istat = nf_get_var1_int(ncid,idvday,ntimes,iday_end) ! ! Get delta time between histories (minutes), model version, ! and history type: istat = nf_get_att_int (ncid,NF_GLOBAL,"delhist_mins",idelhist) mversion = ' ' istat = nf_get_att_text(ncid,NF_GLOBAL,"model_version",mversion) htype = ' ' istat = nf_get_att_text(ncid,NF_GLOBAL,"history_type",htype) ! ! Construct file info string: ibegiyd = (iyear_beg-iyear_beg/100*100)*1000+iday_beg iendiyd = (iyear_end-iyear_end/100*100)*1000+iday_end write(fileinfo,"(i5.5,' ',i3,' ',i2,' ',i2,' to ', | i5.5,' ',i3,' ',i2,' ',i2,' by ',i4, | ' ',a,' ',a)") ibegiyd, mtime_beg,iendiyd, mtime_end, | idelhist,trim(mversion),trim(htype) ! ! Add file info string and description of file info global attributes: istat = nf_redef(ncid) istat = nf_put_att_text(ncid,NF_GLOBAL,"contents", | len_trim(fileinfo),trim(fileinfo)) write(contents_desc,"('yyddd day hour min', | ' to yyddd day hour min by delta_mins')") istat = nf_put_att_text(ncid,NF_GLOBAL,"contents_desc", | len_trim(contents_desc),trim(contents_desc)) istat = nf_enddef(ncid) end subroutine nc_fileinfo !------------------------------------------------------------------- subroutine handle_ncerr(istat,msg) ! ! 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 !------------------------------------------------------------------- end module nchist_module