MODULE process_domain_module USE input_module USE output_module USE module_model_basics USE module_arrays integer :: next_file logical :: run_out_of_files logical :: file_is_open integer :: vertical !!! needed for v5dcreation CONTAINS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: process_domain ! Purpose: Process the input data !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE process_domain () USE date_pack USE gridinfo_module USE misc_definitions_module USE module_debug USE v5d_module IMPLICIT NONE ! ------ some vis5d stuff, contain MAXVARS,MAXTIMES,MAXROWS,MAXCOLUMNS,MAXLEVELS include 'v5df.h' ! Local variables integer :: idiff, n_times, i character (len=19) :: valid_date, temp_date, datestr character (len=128) :: cname, stagger, cunits, cdesc integer :: istatus, reclength integer :: ret !!! v5d stuff integer, dimension(MAXLEVELS) :: v5d_nl real, dimension(MAXLEVELS) :: v5d_levels integer, dimension(MAXTIMES) :: timestamp, datestamp integer :: projection real, dimension(100) :: proj_args ! Compute number of times that we will process CALL geth_idts(end_date, start_date, idiff) CALL mprintf((idiff < 0),ERROR,'Ending date is earlier than starting date in namelist.') ! Check that the interval evenly divides the range of times to process n_times = idiff / interval_seconds CALL mprintf((mod(idiff, interval_seconds) /= 0),WARN, & 'In namelist, interval_seconds does not evenly divide '// & '(end_date - start_date). Only %i time periods '// & 'will be processed.', i2=n_times) ! ! DO TIME_INDEPENDANT PROCESSING ! ! Initialize the input module to read input data CALL input_init(1, istatus) CALL mprintf((istatus /= 0),ERROR, 'input_init(): Error opening input.') ! Read global attributes from the input file CALL read_global_attrs () IF ( iprogram .le. 1 ) n_times = 0 IF ( output_type == 'grads' ) THEN !!! We are making grads type files ! Open .dat and .ctl files DO cunit=10,100 INQUIRE(unit=cunit, opened=is_used) if (.not. is_used) exit END DO ctlfile = trim(output_root_name)//'.ctl' OPEN(cunit,file=ctlfile) DO dunit=10,100 INQUIRE(unit=dunit, opened=is_used) IF (.not. is_used) EXIT END DO IF ( west_east_dim .gt. 2 .AND. south_north_dim .gt. 2 ) THEN reclength = (west_east_dim)*(south_north_dim) ELSE IF ( west_east_dim .gt. 2 .AND. south_north_dim .eq. 2 ) THEN reclength = west_east_dim ELSE IF ( west_east_dim .eq. 2 .AND. south_north_dim .gt. 2 ) THEN reclength = south_north_dim ELSE IF ( west_east_dim .eq. 2 .AND. south_north_dim .eq. 2 ) THEN reclength = 1 END IF #ifdef RECL4 reclength = reclength*4 #endif datfile = trim(output_root_name)//'.dat' ! OPEN(dunit,file=datfile,form='unformatted',access="direct", recl=reclength) OPEN(dunit,file=datfile) ELSE IF ( output_type == 'v5d' ) THEN !!! We are making grads type files v5dfile = trim(output_root_name)//'.v5d' END IF CALL arw_get_next_time(handle, datestr, istatus) !!! !!! !!! MARS: we go to next time to init, because first step of PHTOT and PTOT are not filled !!! MARS: caca caca caca. il vaut mieux le geop init pour calculer les hauteurs !!! !CALL arw_get_next_time(handle, datestr, istatus) vertical = 0 !!! needed for v5d creation (equally spaced levels in generic units) ! If interp_method /= 0, then we need to know if this can be done, for now only doing it for wrfout data IF ( iprogram .ge. 6 .AND. interp_method /= 0 ) THEN CALL get_interp_info (datestr) ELSE vertical_type = 'n' !! In case user set this for other programs number_of_zlevs = bottom_top_dim END IF #ifdef V5D IF ( output_type == 'v5d' ) THEN !!! Create the vis5d file v5d_nl = IMISSING CALL v5d_varnames (v5d_nl) CALL v5d_times (n_times+1, timestamp, datestamp) CALL v5d_proj (datestr, projection, proj_args) v5d_levels = MISSING IF (vertical == 0 ) THEN v5d_levels(1) = 0.0 !!! Same basics set up for sigma level data v5d_levels(2) = 1.0 ELSE DO i=1,number_of_zlevs v5d_levels(i) = interp_levels(i) END DO END IF ret = v5dCreate( v5dfile, n_times+1, numvars, & !!west_east_dim, south_north_dim, v5d_nl, & south_north_dim, west_east_dim, v5d_nl, & varnames, timestamp, datestamp, 1, projection, & proj_args, vertical, v5d_levels ) IF ( ret == 0 ) THEN CALL mprintf(.true.,ERROR, 'v5dCreate') ELSE CALL mprintf(.true.,STDOUT, ' ') CALL mprintf(.true.,STDOUT, ' SUCCESS: v5d file has been Created') END IF END IF #endif CALL input_close() !!! Will be opened again by next process if needed ! ! BEGIN TIME-DEPENDANT PROCESSING next_file = 1 ivars = 0 ! Loop over all times to be processed could_not_find = plot_these_fields !!! Keep a list of requested fields run_out_of_files = .FALSE. file_is_open = .FALSE. rec = 0 all_files : DO time=0,n_times IF ( iprogram == 1 ) then temp_date = "0000-00-00_00:00:00" ELSE CALL geth_newdate(valid_date, trim(start_date), time*interval_seconds) temp_date = ' ' write(temp_date,'(a19)') valid_date(1:10)//'_'//valid_date(12:19) ENDIF IF ( time == 0 ) tdef_date = temp_date CALL get_fields(temp_date) IF (run_out_of_files) EXIT all_files ENDDO all_files ! Loop over n_times CALL mprintf(.true.,STDOUT, ' ') CALL mprintf(.true.,STDOUT, 'DONE Processing Data') IF ( output_type == 'grads' ) THEN !!! We are making grads type files ! Create the .ctl file CALL mprintf(.true.,STDOUT, ' ') CALL mprintf(.true.,STDOUT, 'CREATING .ctl file') CALL create_ctl () END IF #ifdef V5D IF ( output_type == 'v5d' ) THEN !!! Create the vis5d file ret = v5dclose() END IF #endif END SUBROUTINE process_domain !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: get_interp_info ! Purpose: Get extra info if we need to interpolate to pressure / height !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE get_interp_info (valid_date) IMPLICIT NONE ! Arguments character (len=19) :: valid_date ! Local variables integer :: istatus, wrftype real, pointer, dimension(:,:,:) :: real_array character (len=3) :: memorder character (len=128) :: cname, stagger character (len=128), dimension(3) :: dimnames integer, dimension(3) :: domain_start, domain_end integer, dimension(2) :: loc_of_min_z real, allocatable, dimension(:,:,:) :: ph_tmp, phb_tmp integer :: found found = 0 wrftype = 104 memorder = "XYZ" domain_start = 1 domain_end(1) = west_east_dim domain_end(2) = south_north_dim IF (iprogram == 6 .AND. vertical_type == 'p') THEN domain_end(3) = bottom_top_dim stagger = '-' istatus = 0 CALL read_spec_field(domain_start, domain_end, 'QVAPOR', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) found = found + istatus memorder = "XY" domain_end(3) = 1 istatus = 0 CALL read_spec_field(domain_start, domain_end, 'MU', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) found = found + istatus istatus = 0 CALL read_spec_field(domain_start, domain_end, 'MUB', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) found = found + istatus memorder = "0" domain_end(1) = 1 domain_end(2) = 1 istatus = 0 CALL read_spec_field(domain_start, domain_end, 'P_TOP', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) found = found + istatus memorder = "Z" domain_end(1) = bottom_top_dim istatus = 0 CALL read_spec_field(domain_start, domain_end, 'ZNU', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) found = found + istatus domain_end(1) = bottom_top_dim + 1 stagger = 'Z' istatus = 0 CALL read_spec_field(domain_start, domain_end, 'ZNW', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) found = found + istatus IF ( found == 0 ) THEN vertical = 3 !!! (unequally spaced levels in mb - v5d) CALL mprintf(.true.,STDOUT, ' Interpolating to PRESSURE levels ') CALL mprintf(.true.,STDOUT, ' ') ELSE CALL mprintf(.true.,STDOUT, ' WARNING: Asked to interpolate to PRESSURE, but we do not have enough information ' ) CALL mprintf(.true.,STDOUT, ' Will output data on MODEL LEVELS ' ) vertical_type = 'n' ENDIF ENDIF IF (iprogram == 8 .AND. vertical_type == 'p') THEN domain_end(3) = bottom_top_dim stagger = '-' istatus = 0 CALL read_spec_field(domain_start, domain_end, 'P', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) found = found + istatus istatus = 0 CALL read_spec_field(domain_start, domain_end, 'PB', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) found = found + istatus IF ( found == 0 ) THEN vertical = 3 !!! (unequally spaced levels in mb - v5d) CALL mprintf(.true.,STDOUT, ' Interpolating to PRESSURE levels ') CALL mprintf(.true.,STDOUT, ' ') ELSE CALL mprintf(.true.,STDOUT, ' WARNING: Asked to interpolate to PRESSURE, but we do not have enough information ' ) CALL mprintf(.true.,STDOUT, ' Will output data on MODEL LEVELS ' ) vertical_type = 'n' ENDIF ENDIF IF (vertical_type == 'z' .AND. interp_method == 1 ) THEN !!! !!! MARS MARS: caca caca caca: corriger ce probleme de staggering ? !!! domain_end(3) = bottom_top_dim + 1 stagger = 'Z' ! domain_end(3) = bottom_top_dim ! stagger = '' !!! !!! MARS MARS: caca caca caca !!! istatus = 0 CALL read_spec_field(domain_start, domain_end, 'PH', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) found = found + istatus istatus = 0 CALL read_spec_field(domain_start, domain_end, 'PHB', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) found = found + istatus !!!print *, 'get_interp_info' ! ! istatus = 0 ! CALL read_spec_field(domain_start, domain_end, 'PHTOT', wrftype, & ! memorder, stagger, dimnames, real_array, valid_date, istatus) ! found = found + istatus IF ( found == 0 ) THEN vertical = 2 !!! (unequally spaced levels in km - v5d) CALL mprintf(.true.,STDOUT, ' Interpolating to USER SPECIFIED HEIGHT levels ') CALL mprintf(.true.,STDOUT, ' ') ELSE CALL mprintf(.true.,STDOUT, ' WARNING: Asked to interpolate& & to USER SPECIFIED HEIGHT, but we do not have enough information ' ) CALL mprintf(.true.,STDOUT, ' Will output data on MODEL LEVELS ' ) vertical_type = 'n' ENDIF ENDIF IF (vertical_type == 'z' .AND. interp_method == -1 ) THEN domain_end(3) = bottom_top_dim + 1 stagger = 'Z' istatus = 0 CALL read_spec_field(domain_start, domain_end, 'PH', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) found = found + istatus IF ( istatus == 0 ) THEN allocate(ph_tmp(domain_end(1),domain_end(2),domain_end(3)-1)) ph_tmp = 0.5*(real_array(:,:,1:domain_end(3)-1)+real_array(:,:,2:domain_end(3))) ENDIF istatus = 0 CALL read_spec_field(domain_start, domain_end, 'PHB', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) found = found + istatus IF ( istatus == 0 ) THEN allocate(phb_tmp(domain_end(1),domain_end(2),domain_end(3)-1)) phb_tmp = 0.5*(real_array(:,:,1:domain_end(3)-1)+real_array(:,:,2:domain_end(3))) ENDIF IF ( found == 0 ) THEN vertical = 2 !!! (unequally spaced levels in km - v5d) !! ph_tmp = ( (ph_tmp+phb_tmp) / 9.81)/1000. !!! convert to height in km ph_tmp = ( (ph_tmp+phb_tmp) / 3.72)/1000. !!! convert to height in km !!! Generate nice heights. !!! Get location of lowest height on model grid !!! Use column above this point as our heights to interpolate to !!! Adjust lowerst and heighest heights a bit number_of_zlevs = bottom_top_dim loc_of_min_z = minloc(ph_tmp(:,:,1)) interp_levels(1:number_of_zlevs) = & ph_tmp(loc_of_min_z(1),loc_of_min_z(2),1:number_of_zlevs) !!--- .002 problematic on Mars with low pressure !!interp_levels(1) = interp_levels(1) + 0.002 !!interp_levels(1) = MAX(interp_levels(1), interp_levels(2)/2.0) !! no neg values !!interp_levels(number_of_zlevs) = interp_levels(number_of_zlevs) - 0.002 vertical = 2 CALL mprintf(.true.,STDOUT, ' Interpolating to GENERATED HEIGHT levels ') CALL mprintf(.true.,STDOUT, ' ') ELSE CALL mprintf(.true.,STDOUT, ' WARNING: Asked to interpolate to GENERATED & &HEIGHT, but we do not have enough information ' ) CALL mprintf(.true.,STDOUT, ' Will output data on MODEL LEVELS ' ) vertical_type = 'n' ENDIF IF (ALLOCATED(ph_tmp)) DEALLOCATE(ph_tmp) IF (ALLOCATED(phb_tmp)) DEALLOCATE(phb_tmp) ENDIF END SUBROUTINE get_interp_info !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: get_interp_array ! Purpose: Get array we will use when interpolting !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE get_interp_array (valid_date) USE module_pressure IMPLICIT NONE ! Arguments character (len=19) :: valid_date ! Local variables integer :: istatus, wrftype real, pointer, dimension(:,:,:) :: real_array character (len=3) :: memorder character (len=128) :: cname, stagger character (len=128), dimension(3) :: dimnames integer, dimension(3) :: domain_start, domain_end integer, dimension(2) :: loc_of_min_z real, allocatable, dimension(:,:,:) :: p_tmp, pb_tmp, ph_tmp, phb_tmp real, pointer, dimension(:,:,:) :: SCR3 wrftype = 104 memorder = "XYZ" domain_start = 1 domain_end(1) = west_east_dim domain_end(2) = south_north_dim IF (iprogram == 6 .AND. vertical_type == 'p') THEN domain_end(3) = bottom_top_dim stagger = '-' CALL read_spec_field(domain_start, domain_end, 'QVAPOR', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) allocate(QV(west_east_dim,south_north_dim,bottom_top_dim)) QV = real_array memorder = "XY" domain_end(3) = 1 stagger = '-' CALL read_spec_field(domain_start, domain_end, 'MU', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) allocate(MU(west_east_dim,south_north_dim)) MU = real_array(:,:,1) CALL read_spec_field(domain_start, domain_end, 'MUB', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) allocate(MUB(west_east_dim,south_north_dim)) MUB = real_array(:,:,1) memorder = "0" domain_end(1) = 1 domain_end(2) = 1 CALL read_spec_field(domain_start, domain_end, 'P_TOP', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) PTOP = real_array(1,1,1) memorder = "Z" domain_end(1) = bottom_top_dim CALL read_spec_field(domain_start, domain_end, 'ZNU', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) allocate(ZNU(bottom_top_dim)) ZNU = real_array(:,1,1) domain_end(1) = bottom_top_dim + 1 stagger = 'Z' CALL read_spec_field(domain_start, domain_end, 'ZNW', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) allocate(ZNW(bottom_top_dim+1)) ZNW = real_array(:,1,1) CALL pressure(SCR3) ALLOCATE(vert_array(west_east_dim,south_north_dim,bottom_top_dim)) vert_array = 0.01 * SCR3 !! pressure array in HPa IF (ALLOCATED(QV)) DEALLOCATE(QV) IF (ALLOCATED(MU)) DEALLOCATE(MU) IF (ALLOCATED(MUB)) DEALLOCATE(MUB) IF (ALLOCATED(ZNU)) DEALLOCATE(ZNU) IF (ALLOCATED(ZNW)) DEALLOCATE(ZNW) ENDIF IF (iprogram == 8 .AND. vertical_type == 'p') THEN domain_end(3) = bottom_top_dim stagger = '-' CALL read_spec_field(domain_start, domain_end, 'P', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) allocate(p_tmp(west_east_dim,south_north_dim,bottom_top_dim)) p_tmp = real_array CALL read_spec_field(domain_start, domain_end, 'PB', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) allocate(pb_tmp(west_east_dim,south_north_dim,bottom_top_dim)) pb_tmp = real_array IF (ALLOCATED(vert_array)) DEALLOCATE(vert_array) ALLOCATE(vert_array(west_east_dim,south_north_dim,bottom_top_dim)) vert_array = (p_tmp+pb_tmp)/100. !! pressure array in HPa IF (ALLOCATED(p_tmp)) DEALLOCATE(p_tmp) IF (ALLOCATED(pb_tmp)) DEALLOCATE(pb_tmp) ENDIF IF (vertical_type == 'z' ) THEN !!! !!! MARS MARS: caca caca caca: corriger ce probleme de staggering ? !!! domain_end(3) = bottom_top_dim + 1 stagger = 'Z' ! domain_end(3) = bottom_top_dim ! stagger = '' !!! !!! MARS MARS: caca caca caca !!! CALL read_spec_field(domain_start, domain_end, 'PH', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) allocate(ph_tmp(west_east_dim,south_north_dim,bottom_top_dim)) ph_tmp = 0.5*(real_array(:,:,1:domain_end(3)-1)+real_array(:,:,2:domain_end(3))) CALL read_spec_field(domain_start, domain_end, 'PHB', wrftype, & memorder, stagger, dimnames, real_array, valid_date, istatus) allocate(phb_tmp(west_east_dim,south_north_dim,bottom_top_dim)) phb_tmp = 0.5*(real_array(:,:,1:domain_end(3)-1)+real_array(:,:,2:domain_end(3))) ! CALL read_spec_field(domain_start, domain_end, 'PHTOT', wrftype, & ! memorder, stagger, dimnames, real_array, valid_date, istatus) ! !!print *, 'get_interp_array' !!!! !!!! MARS MARS: caca caca caca !!!! ! allocate(ph_tmp(west_east_dim,south_north_dim,bottom_top_dim)) ! !!ph_tmp = 0.5*0.5*(real_array(:,:,1:domain_end(3)-1)+real_array(:,:,2:domain_end(3))) !ph_tmp = 0.5*real_array ! ! allocate(phb_tmp(west_east_dim,south_north_dim,bottom_top_dim)) ! !!phb_tmp = 0.5*0.5*(real_array(:,:,1:domain_end(3)-1)+real_array(:,:,2:domain_end(3))) !phb_tmp = 0.5*real_array !!!! !!!! MARS MARS: caca caca caca !!!! IF (ALLOCATED(vert_array)) DEALLOCATE(vert_array) ALLOCATE(vert_array(west_east_dim,south_north_dim,bottom_top_dim)) !! vert_array = ( (ph_tmp+phb_tmp) / 9.81)/1000. !! height array in km vert_array = ( (ph_tmp+phb_tmp) / 3.72)/1000. !! height array in km IF (ALLOCATED(ph_tmp)) DEALLOCATE(ph_tmp) IF (ALLOCATED(phb_tmp)) DEALLOCATE(phb_tmp) ENDIF END SUBROUTINE get_interp_array !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: get_fields ! Purpose: Read all fields in input file and process required output fields !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE get_fields (valid_date) USE module_get_file_names USE module_diagnostics USE module_interp IMPLICIT NONE ! Local variables integer :: istatus, i, j, k, ii, jj, kk real, pointer, dimension(:,:,:) :: real_array character (len=3) :: memorder character (len=19) :: datestr, valid_date character (len=128) :: cname, stagger, cunits, cdesc character (len=128), dimension(3) :: dimnames integer, dimension(3) :: domain_start, domain_end real, pointer, dimension(:,:,:) :: data_out real, allocatable, dimension(:,:,:):: SCR integer :: nxout, nyout, nzout character (len=20) :: dummy integer :: is_there ! Initialize the input module to read static input data for this domain CALL mprintf(.true.,STDOUT, ' ') CALL mprintf(.true.,STDOUT, ' Processing time --- %s', s1=trim(valid_date)) get_right_time : DO IF ( next_file > number_of_input_files ) THEN CALL mprintf(.true.,STDOUT, ' --- WE HAVE RUN OUT OF INPUT FILES ---') run_out_of_files = .TRUE. RETURN END IF IF ( .not. file_is_open) THEN CALL input_init(next_file, istatus) CALL mprintf((istatus /= 0),ERROR, 'input_init(): Error opening input.') file_is_open = .TRUE. END IF CALL arw_get_next_time(handle, datestr, istatus) IF ( istatus /= 0 ) THEN ! might be in a next file CALL input_close() file_is_open = .FALSE. CALL mprintf(.true.,STDOUT, ' Date not in this file - see if there are more files ') next_file = next_file + 1 CYCLE get_right_time END IF IF ( TRIM(datestr) .LT. TRIM(valid_date) ) THEN CALL mprintf(.true.,STDOUT, ' Not looking for %s - skipping to next time' , s1=datestr ) CYCLE get_right_time ELSEIF ( TRIM(datestr) .EQ. TRIM(valid_date) ) THEN CALL mprintf(.true.,STDOUT, ' Found the right date - continue ' ) EXIT get_right_time ELSEIF ( TRIM(datestr) .GT. TRIM(valid_date) ) THEN CALL mprintf(.true.,STDOUT, ' Found %s before %s', s1=trim(datestr),s2=trim(valid_date)) run_out_of_files = .TRUE. RETURN ENDIF ENDDO get_right_time CALL clobber_arrays #ifdef IO_GRIB1 IF (io_form_input == GRIB1) THEN REWIND (13) END IF #endif !! IF we are interpolting we need the pressure/height array to interpolate to IF ( vertical_type /= 'n' ) THEN CALL get_interp_array (valid_date) ENDIF ! Read fields using the input module; we know that there are no more ! fields to be read when read_next_field() returns a non-zero status. istatus = 0 process_all_fields : DO WHILE (istatus == 0) CALL read_next_field(domain_start, domain_end, cname, cunits, cdesc, & memorder, stagger, dimnames, real_array, valid_date, istatus) #ifdef IO_GRIB1 IF (io_form_input == GRIB1) THEN IF (istatus == -5) THEN istatus = 0 CYCLE process_all_fields END IF END IF #endif IF (istatus == 0) THEN IF (cname(len_trim(cname)-1:len_trim(cname)) == "_U" ) CYCLE process_all_fields IF (cname(len_trim(cname)-1:len_trim(cname)) == "_V" ) CYCLE process_all_fields CALL keep_arrays(cname, real_array) !!! Keep some fields around !! !! DECIDE IF WE WANT THE FIELD AND IF YES, INTERPOLATE TO CORRECT GRID IF NEEEDED !! assume we got the grid info up at the top somewhere !! dummy = ","//trim(cname)//"," is_there = INDEX(plot_these_fields,trim(dummy)) IF ( ( INDEX(plot,'all') /= 0 .OR. is_there /= 0) .AND. domain_end(2) /= 1) THEN IF (ALLOCATED(SCR)) DEALLOCATE(SCR) ALLOCATE(SCR(domain_end(1),domain_end(2),domain_end(3))) SCR = real_array CALL interp( SCR, domain_end(1), domain_end(2), domain_end(3), & data_out, nxout, nyout, nzout, & vert_array, interp_levels, number_of_zlevs) ! Write the fields we want out to the .dat file, also keeps a list of what is written out CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits) ENDIF ENDIF !! end "istatus==0" END DO process_all_fields IF ( INDEX(plot,'list') /= 0 .OR. INDEX(plot,'file') /= 0 ) THEN !! Do we have any DIAGNOSTICS to process? CALL process_diagnostics () END IF !! We are DONE for this time !! print a list of the requested fields we could not find IF ( len_trim(could_not_find) > 1 ) THEN CALL mprintf(.true.,STDOUT, ' ') CALL mprintf(.true.,STDOUT, ' WARNING: The following requested fields could not be found ' ) could_not_find = could_not_find(2:len_trim(could_not_find)) is_there = INDEX(could_not_find,",") #ifdef V5D IF ( output_type == 'v5d' .AND. is_there > 1 ) THEN !!! Warning for possible bad data file CALL mprintf(.true.,STDOUT, ' WARNING: This will create a problem for VIS5D, please ' ) CALL mprintf(.true.,STDOUT, ' remove these fields from your input list and ' ) CALL mprintf(.true.,STDOUT, ' run the code again. ' ) CALL mprintf(.true.,STDOUT, ' ' ) END IF #endif DO WHILE ( is_there > 1 ) PRINT*," ", could_not_find(1:is_there-1) could_not_find = could_not_find(is_there+1:len_trim(could_not_find)) is_there = INDEX(could_not_find,",") END DO CALL mprintf(.true.,STDOUT, ' ') ENDIF END SUBROUTINE get_fields END MODULE process_domain_module