MODULE output_module USE input_module USE module_model_basics USE module_arrays integer :: time, rec, ivars !!! process time ; rec in .dat file ; variables in .dat file integer :: cunit, dunit !!! .ctl and .dat output files character (len=128) :: ctlfile, datfile, v5dfile character (len=200), dimension(200) :: ctl_var_string !!! List of fields foe .ctl file character (len=2000) :: could_not_find !!! Keep list of what we wrote out character (len=19) :: tdef_date !!! Output start date character (len=40) :: tdef !!! tdef string in .ctl file character (len=10), dimension(100) :: varnames integer :: numvars CONTAINS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: create_ctl ! Purpose: Create the .ctl file !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE create_ctl () IMPLICIT NONE ! Local variables real, dimension(4) :: xlat_a, xlon_a real :: dxll, slack, temp real :: abslatmin, abslonmax, abslonmin, abslatmax integer :: i, k, ipoints, jpoints, ilon integer :: ipole real :: xi, r, re, xj, dx_ps, earthr, radpd real :: rlat, rlong, ref_lon, ref_lat integer :: file_start, file_end integer :: we_dim, sn_dim file_start = INDEX(datfile,"/",BACK=.TRUE.) file_end = len_trim(datfile) write(cunit,'("dset ^",A)') datfile(file_start+1:file_end) #ifdef bytesw write (cunit, '("options byteswapped")') #endif write(cunit, '("undef 1.e37")') IF ( output_title /= ' ') THEN write(cunit,'("title ",A)') trim(output_title) ELSE write(cunit,'("title ",A)') trim(title) END IF IF (map_proj == 1) THEN !! Lambert Projection !!!!!!! !!!!!!! TEMP !!!!!!! IF ( mercator_defs ) THEN WRITE(cunit,'("xdef ",i4," levels ")') west_east_dim DO i = 1,west_east_dim WRITE(cunit,*) XLONG(i,int(west_east_dim/2)) ENDDO WRITE(cunit,'("ydef ",i4," levels ")') south_north_dim DO i = 1,south_north_dim WRITE(cunit,*) XLAT(int(south_north_dim/2),i) END DO ELSE ! Make sure truelat1 is the larger number IF (truelat1 < truelat2) THEN temp = truelat1 truelat1 = truelat2 truelat2 = temp END IF !!WRITE(cunit,'("pdef ",i4," ",i3," lcc ",f9.5," ",f10.5," ", & !! rounding off WRITE(cunit,'("pdef ",i4," ",i3," lcc ",f7.3," ",f8.3," ", & & f8.3," ",f8.3," ",f9.5," ",f9.5," ",f10.5," ", & & f10.3," ",f10.3)') & west_east_dim,south_north_dim,cen_lat,cen_lon, & (west_east_dim+1)/2.,(south_north_dim+1)/2., & truelat1,truelat2,stand_lon,dx,dy xlon_a(1) = XLONG(1,1) xlon_a(2) = XLONG(1,south_north_dim) xlon_a(3) = XLONG(west_east_dim,1) xlon_a(4) = XLONG(west_east_dim,south_north_dim) !check for dateline ilon = 0 if ( abs(xlon_a(1) - xlon_a(3)) .GT. 180. ) ilon = 1 if ( abs(xlon_a(2) - xlon_a(4)) .GT. 180. ) ilon = 1 abslatmin = minval(XLAT) abslatmax = maxval(XLAT) abslonmin = 99999. abslonmax = -99999. DO i=1,4 IF ( xlon_a(i) .lt. 0.0 .AND. ilon .eq. 1 ) THEN abslonmin=min(abslonmin,360.+xlon_a(i)) abslonmax=max(abslonmax,360.+xlon_a(i)) ELSE abslonmin=min(abslonmin,xlon_a(i)) abslonmax=max(abslonmax,xlon_a(i)) ENDIF END DO ! dxll = (dx/1000.)/111./2. dxll = (dx/1000.)/59./2. slack = ((dx/1000.)*3.)/100. ipoints = int(((abslonmax-abslonmin)+(3*slack))/dxll) jpoints = int(((abslatmax-abslatmin)+(3*slack))/dxll) WRITE(cunit,'("xdef ",i4," linear ",f10.5," ",f12.8)') ipoints, & (abslonmin-1.5*slack),dxll WRITE(cunit,'("ydef ",i4," linear ",f10.5," ",f12.8)') jpoints, & (abslatmin-1.5*slack),dxll !!!!!!!!!!!!! !!!!!!!!!!!!! TEMP !!!!!!!!!!!!! ENDIF ENDIF IF (map_proj == 2) THEN !! Polar Stereo Projection ipole = 1 IF (truelat1 .lt. 0.0) ipole = -1 radpd = .01745329 ! earthr = 6.3712E6 earthr = 3.3972E6 dx_ps = dx/( (1.0+sin(ipole*TRUELAT1*radpd))/(1.0+sin(60*radpd)) ) !! true dx (meter) at 60degrees re = ( earthr * (1.0 + sin(60*radpd)) ) / dx_ps !! IF ( ipole == 1 ) THEN rlong = (XLONG(1,1) - stand_lon) * radpd rlat = XLAT(1,1) * radpd r = (re * cos(rlat)) / (1.0 + sin(rlat)) xi = (1. - r * sin(rlong)) xj = (1. + r * cos(rlong)) write(cunit,'("pdef ",i3," ",i3," nps ",f8.3," ",f8.3," ", & & f12.4," ",f12.7)') west_east_dim, south_north_dim, & xi, xj, stand_lon, dx_ps*0.001 ELSE IF ( ipole == -1 ) THEN rlong = (180.0 + XLONG(1,1) - stand_lon) * radpd rlat = ipole*XLAT(1,1) * radpd r = (re * cos(rlat)) / (1.0 + sin(rlat)) xi = (1. - ipole * r * sin(rlong)) xj = (1. + r * cos(rlong)) write(cunit,'("pdef ",i3," ",i3," sps ",f8.3," ",f8.3," ", & & f12.4," ",f12.7)') west_east_dim, south_north_dim, & xi, xj, (180+stand_lon), -0.001*dx_ps END IF abslonmin = minval(XLONG) abslonmax = maxval(XLONG) abslatmin = minval(XLAT) abslatmax = maxval(XLAT) !! dxll = (dx_ps/1000.)/111./2. dxll = (dx_ps/1000.)/59./2. ipoints = int(((abslonmax-abslonmin)+2.)/dxll)-2 jpoints = int(((abslatmax-abslatmin)+2.)/dxll)-3 ref_lon = abslonmin-1. IF ( abslonmax .ge. 175. .AND. abslonmax .le. 180.) THEN IF ( abslonmin .ge. -180. .AND. abslonmin .le. -175.) THEN ref_lon = 0.0 END IF END IF ref_lat = abslatmin-1. !!IF ( ipole == -1 ) ref_lat = abslatmax+1. IF ( ipole == -1 ) ref_lat = abslatmin WRITE(cunit,'("xdef ",i4," linear ",f6.1," ",f12.8)') ipoints, & ref_lon,dxll WRITE(cunit,'("ydef ",i4," linear ",f6.1," ",f12.8)') jpoints, & ref_lat,dxll ENDIF IF (map_proj == 3) THEN !! Mercator Projection !check for dateline ilon = 0 IF ( abs(XLONG(1,1) - XLONG(west_east_dim,1)) .GT. 180. ) ilon = 1 IF ( abs(XLONG(1,south_north_dim) - XLONG(west_east_dim,south_north_dim)) .GT. 180. ) ilon = 1 IF ( mercator_defs .AND. map_proj == 3 ) THEN WRITE(cunit,'("xdef ",i4," levels ")') west_east_dim DO i = 1,west_east_dim WRITE(cunit,*) XLONG(i,1) !!IF ( XLONG(i,1) .lt. 0.0 ) THEN !!WRITE(cunit,*) XLONG(i,1) + 180.0 !!ELSE !!WRITE(cunit,*) XLONG(i,1) !!END IF END DO ELSE IF ( ilon == 1 ) THEN WRITE(cunit,'("xdef ",i4," linear ",f9.4," ",f8.4)') & west_east_dim,XLONG(1,1), & (abs(XLONG(1,1)-(360.+XLONG(west_east_dim,south_north_dim)))/(west_east_dim-1)) ELSE WRITE(cunit,'("xdef ",i4," linear ",f9.4," ",f8.4)') & west_east_dim,XLONG(1,1), & (abs(XLONG(1,1)-XLONG(west_east_dim,south_north_dim))/(west_east_dim-1)) ENDIF ENDIF IF ( mercator_defs .AND. map_proj == 3 ) THEN WRITE(cunit,'("ydef ",i4," levels ")') south_north_dim DO i = 1,south_north_dim WRITE(cunit,*) XLAT(1,i) END DO ELSE WRITE(cunit,'("ydef ",i4," linear ",f9.4," ",f8.4)') & south_north_dim,XLAT(1,1),(abs(XLAT(1,1)-XLAT(west_east_dim,south_north_dim))/(south_north_dim-1)) ENDIF ENDIF IF (map_proj == 0) THEN !! Ideal we_dim = west_east_dim sn_dim = south_north_dim IF (we_dim == 2) we_dim = 1 IF (sn_dim == 2) sn_dim = 1 write (cunit, '("xdef ",i4, " linear 0 0.0001")') we_dim write (cunit, '("ydef ",i4, " linear 0 0.0001")') sn_dim ENDIF !! END of projection specific section IF (iprogram .ge. 6 .and. (vertical_type=='p'.or.vertical_type=='z')) THEN WRITE(cunit,'("zdef ",i3, " levels ")') number_of_zlevs DO k = 1,number_of_zlevs WRITE(cunit,'(f10.5)') interp_levels(k) END DO ELSE WRITE(cunit,'("zdef ",i4," linear 1 1 ")') bottom_top_dim ENDIF IF (tdef_date == '0000-00-00_00:00:00') THEN tdef = ' 1 linear 00z01jan2000 1hr' ELSE IF ( output_type == 'grads' ) THEN CALL tdef_calc_grads () ELSE CALL tdef_calc () ENDIF END IF WRITE(cunit,'("tdef ",A)') tdef WRITE(cunit,'("VARS ",i3)') ivars DO i = 1,ivars WRITE(cunit,'(A)') trim(ctl_var_string(i)) ENDDO WRITE(cunit,'("ENDVARS")') DO i=1,iatts WRITE(cunit,'(A)') trim(catts(i)) ENDDO END SUBROUTINE create_ctl !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: tdef_calc ! Purpose: Make the tdef line for the .ctl file !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE tdef_calc () IMPLICIT NONE character (len=3) :: cmonth integer :: year,month,day,hour tdef = ' ' !! Time comes in as YYYY-MM-DD_HH:MM:SS !! 1234 67 90 23 45 89 READ(tdef_date(1:4),*) year READ(tdef_date(6:7),*) month READ(tdef_date(9:10),*) day READ(tdef_date(12:13),*) hour IF (month == 1) cmonth = 'M1' IF (month == 2) cmonth = 'M2' IF (month == 3) cmonth = 'M3' IF (month == 4) cmonth = 'M4' IF (month == 5) cmonth = 'M5' IF (month == 6) cmonth = 'M6' IF (month == 7) cmonth = 'M7' IF (month == 8) cmonth = 'M8' IF (month == 9) cmonth = 'M9' IF (month ==10) cmonth = 'M10' IF (month ==11) cmonth = 'M11' IF (month ==12) cmonth = 'M12' IF ( year < 100 ) year = 2000 WRITE (tdef,'(i4," linear ",i2.2,"Z",i2.2,A,i4," ",i6,"MN")') time, hour, day, cmonth, year, interval_seconds END SUBROUTINE tdef_calc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: tdef_calc_grads ! Purpose: Make the tdef line for the .ctl file !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE tdef_calc_grads () IMPLICIT NONE character (len=3) :: cmonth integer :: year,month,day,hour tdef = ' ' !! Time comes in as YYYY-MM-DD_HH:MM:SS !! 1234 67 90 23 45 89 READ(tdef_date(1:4),*) year READ(tdef_date(6:7),*) month READ(tdef_date(9:10),*) day READ(tdef_date(12:13),*) hour IF (month == 1) cmonth = 'jan' IF (month == 2) cmonth = 'fev' IF (month == 3) cmonth = 'mar' IF (month == 4) cmonth = 'avr' IF (month == 5) cmonth = 'may' IF (month == 6) cmonth = 'jun' IF (month == 7) cmonth = 'jul' IF (month == 8) cmonth = 'aug' IF (month == 9) cmonth = 'sep' IF (month ==10) cmonth = 'oct' IF (month ==11) cmonth = 'nov' IF (month ==12) cmonth = 'dec' IF ( year < 100 ) year = 2000 WRITE (tdef,'(i4," linear ",i2.2,"Z",i2.2,A,i4," ",i6,"MN")') time, hour, day, cmonth, year, interval_seconds END SUBROUTINE tdef_calc_grads !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: write_dat ! Purpose: Write the data to the .dat file. Also keep record of what has been ! written for creation of the .ctl file !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE write_dat (data_out, nx, ny, nz, cname, cdesc, cunits) IMPLICIT NONE ! ------ some vis5d stuff, contain MAXVARS,MAXTIMES,MAXROWS,MAXCOLUMNS,MAXLEVELS include 'v5df.h' ! Arguments integer :: nx, ny, nz real, dimension(nx,ny,nz) :: data_out character (len=128) :: cname, cunits, cdesc ! Local variables integer :: ii, jj, kk integer :: is_there, ret, ivar character (len=20) :: dummy real, dimension(ny,nx,nz) :: v5d_data_out ! IF ( output_type == 'ncdf' ) THEN !!! Create the netcdf file ! IF ( nx == 2 ) nx = 1 ! IF ( ny == 2 ) ny = 1 ! !DO kk=1,nz ! ! rec = rec + 1 ! ! WRITE(dunit,rec=rec) ((data_out(ii,jj,kk),ii=1,nx),jj=1,ny) ! !END DO ! CALL writencdf(unite nx,ny,nz,"u","W.m-2",3,data_out) ! END IF IF ( output_type == 'grads' ) THEN !!! Create the grads file IF ( nx == 2 ) nx = 1 IF ( ny == 2 ) ny = 1 DO kk=1,nz rec = rec + 1 WRITE(dunit,rec=rec) ((data_out(ii,jj,kk),ii=1,nx),jj=1,ny) END DO END IF !!!!****MARS IDL IF ( output_type == 'idl' ) THEN !!! Create the IDL ascii file IF ( nx == 2 ) nx = 1 IF ( ny == 2 ) ny = 1 DO kk=1,nz rec = rec + 1 !!!****MARS ! WRITE(dunit,rec=rec) ((data_out(ii,jj,kk),ii=1,nx),jj=1,ny) WRITE(dunit,FMT='(F20.8)') ((data_out(ii,jj,kk),ii=1,nx),jj=1,ny) END DO END IF #ifdef V5D IF ( output_type == 'v5d' ) THEN !!! Create the vis5d file DO kk = 1,nz DO jj = 1,nx DO ii = 1,ny v5d_data_out(ii,jj,kk) = data_out(jj,ny-ii+1,kk) END DO END DO END DO ivar = 0 DO ii = 1,numvars IF ( trim(varnames(ii)) == trim(cname) ) THEN ivar = ii END IF END DO ret = v5dwrite ( time+1, ivar, v5d_data_out ) END IF #endif IF (time == 0) THEN IF ( cname(1:3) == 'clf' ) then dummy = ",clfr," ELSE dummy = ","//trim(cname)//"," END IF is_there = INDEX(could_not_find,trim(dummy)) DO WHILE ( is_there /= 0 ) could_not_find = could_not_find(1:is_there)//could_not_find(is_there+len_trim(dummy):len_trim(could_not_find)) is_there = INDEX(could_not_find,trim(dummy)) END DO ivars = ivars + 1 WRITE(ctl_var_string(ivars),'(A," ",i4," 0 ",A," (",A,")")') cname(1:10),nz,trim(cdesc),trim(cunits) ENDIF END SUBROUTINE write_dat ! SUBROUTINE writencdf(unite,nx,ny,nz,"u","W.m-2",3,data_out) ! ! ! ! NETCDF file ! fichnom="diagfi.nc" ! if (firstnom.eq.'1234567890') then ! firstnom = nom ! ierr = NF_CREATE (fichnom, NF_CLOBBER, nid) ! ! Build dimensions ! ierr = NF_DEF_DIM (nid, "Time", NF_UNLIMITED, idim) ! ierr = NF_DEF_DIM (nid, "latitude", ny, idim_lat) ! ierr = NF_DEF_DIM (nid, "longitude", nx, idim_lon) ! ierr = NF_DEF_DIM (nid, "altitude", nz, idim_alt) ! ! Write dimensions attributes ! ierr = NF_DEF_VAR (nid, "Time", NF_FLOAT, 1, idim, varid) ! ierr = NF_PUT_ATT_TEXT (nid, varid, "long_name", 4, "Time") ! ierr = NF_PUT_ATT_TEXT (nid, varid, 'units', 29, "Time") ! ierr = NF_ENDDEF(nid) ! ! ! ierr = NF_DEF_VAR (nid, "latitude", NF_FLOAT, 1, idim_lat, varid) ! ierr = NF_PUT_ATT_TEXT (nid, varid, 'units', 13, "degrees_north") ! ierr = NF_PUT_ATT_TEXT (nid, varid, "long_name", 14, "North latitude") ! ierr = NF_ENDDEF(nid) ! ierr = NF_DEF_VAR (nid, "longitude", NF_FLOAT, 1, idim_lon, varid) ! ierr = NF_PUT_ATT_TEXT (nid, varid, 'units', 12, "degrees_east") ! ierr = NF_PUT_ATT_TEXT (nid, varid, "long_name", 14, "East longitude") ! ierr = NF_ENDDEF(nid) ! ierr = NF_DEF_VAR (nid, "altitude", NF_FLOAT, 1, idim_alt, varid) ! ierr = NF_PUT_ATT_TEXT (nid, varid, 'units', 2, "km") ! ierr = NF_PUT_ATT_TEXT (nid, varid, "long_name", 10, "pseudo-alt") ! ierr = NF_PUT_ATT_TEXT (nid, varid, 'positive', 2, "up") ! ierr = NF_ENDDEF(nid) ! else ! ierr = NF_OPEN(fichnom,NF_WRITE,nid) ! endif ! ! Advance time ! ierr = NF_INQ_VARID (nid, "Time", varid) ! ierr = NF_PUT_VARA_REAL (nid, varid, ntime, 1, date) ! ! Write field ! ierr = NF_INQ_VARID(nid,nom,varid) ! if (ierr /= NF_NOERR) then ! ! Create variable if does not exist ! ierr = NF_INQ_DIMID(nid,"longitude",id(1)) ! ierr = NF_INQ_DIMID(nid,"latitude",id(2)) ! ierr = NF_INQ_DIMID(nid,"altitude",id(3)) ! ierr = NF_INQ_DIMID(nid,"Time",id(4)) ! write (*,*) "=====================================" ! write (*,*) "NETCDF: create ",nom ! call def_var (nid,nom,titre,unite,4,id,varid,ierr) ! endif ! END SUBROUTINE writencdf !!================================================================================== !!================================================================================== ! !subroutine def_var(nid,name,title,units,nbdim,dim,nvarid,ierr) ! !implicit none ! !include "netcdf.inc" ! !character (len=*) :: title,units,name !integer :: nid,nbdim,nvarid,ierr !integer, dimension(nbdim) :: dim ! !ierr=NF_REDEF(nid) !#ifdef NC_DOUBLE !ierr = NF_DEF_VAR (nid,adjustl(name),NF_DOUBLE,nbdim,dim,nvarid) !#else !ierr = NF_DEF_VAR (nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid) !#endif !if(ierr/=NF_NOERR) then ! write(*,*) NF_STRERROR(ierr) ! stop "in def_var" !endif !ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", !len_trim(adjustl(title)),adjustl(title)) !if(ierr/=NF_NOERR) then ! write(*,*) NF_STRERROR(ierr) ! stop "in def_var" !endif !ierr = NF_PUT_ATT_TEXT (nid, nvarid, "units", !len_trim(adjustl(units)),adjustl(units)) !if(ierr/=NF_NOERR) then ! write(*,*) NF_STRERROR(ierr) ! stop "in def_var" !endif !ierr = NF_ENDDEF(nid) ! !end subroutine def_var END MODULE output_module