|
XBeach
|
00001 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00002 !!!!!!!!!!!!!!!!!!!!!!! MODULE OUTPUT !!!!!!!!!!!!!!!!!!!!!!!!!!! 00003 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00004 ! this module will provide for netcdf output in the way varoutput provides for binary output. 00005 ! It should be called optionally through one of the params.txt settings 00006 ! It will also be compiled conditionally like mpi. Only if it proves usefull will it be added by default. 00007 ! it will add dependencies on the netcdf fortran library (http://www.unidata.ucar.edu/software/netcdf/) 00008 ! 00009 ! With contributions from Uwe Rosebrock (CSIRO). 00010 ! 00011 ! if NCSINGLE is defined, output of all real variables in single precision 00012 ! else double precision 00013 ! 00014 ! P+R: tijdelijk uit voor testbed 00015 !#define NCSINGLE 00016 !#ifdef NCSINGLE 00017 !#define NCREAL NF90_REAL 00018 !#define CONVREAL sngl 00019 !#define CONVREALTYPE real*4 00020 !#else 00021 !#define NCREAL NF90_DOUBLE 00022 #define CONVREAL 00023 #define CONVREALTYPE real*8 00024 !#endif 00025 00026 ! NF90: macro to call nf90 function: if return value .ne. nf90_noerr, 00027 ! an error message is produced, including file and lineno of the error, 00028 ! and the program is halted 00029 ! example: 00030 ! replace these two lines: 00031 ! status = nf90_create(path = par%ncfilename, cmode=ior(NF90_CLOBBER,NF90_64BIT_OFFSET), ncid = ncid) 00032 ! if (status /= nf90_noerr) call handle_err(status,__FILE__,__LINE__) 00033 ! with this one line: 00034 ! NF90(nf90_create(path = par%ncfilename, cmode=ior(NF90_CLOBBER,NF90_64BIT_OFFSET), ncid = ncid)) 00035 ! 00036 ! NOTE: the NF90 call must be on one line because of preprocessor restrictions 00037 ! 00038 #ifdef USENETCDF 00039 #define NF90(nf90call) call handle_err(nf90call,__FILE__,__LINE__) 00040 #endif 00041 00042 module ncoutput_module 00043 00044 #ifdef HAVE_CONFIG_H 00045 #include "config.h" 00046 #endif 00047 use xmpi_module 00048 #ifdef USENETCDF 00049 use netcdf 00050 #endif 00051 use mnemmodule, only: mnemonics, numvars 00052 use typesandkinds, only: slen 00053 implicit none 00054 save 00055 private 00056 public ncoutput, fortoutput_init, points_output_init 00057 #ifdef USENETCDF 00058 public ncoutput_init 00059 #endif 00060 00061 ! wwvv todo why the save's? 00062 ! see http://stackoverflow.com/questions/2893097/fortran-save-statement 00063 ! so, we should use save everywhere in modules, or use them in main program 00064 ! we could also add a line 00065 ! save 00066 ! after 'implicit none' 00067 00068 integer :: ncid 00069 00070 ! parameters 00071 integer :: parvarid 00072 00073 ! grid 00074 integer :: xdimid, ydimid 00075 integer :: xvarid, yvarid 00076 00077 ! Wave angle 00078 integer :: thetadimid 00079 ! Sediment 00080 integer :: sedimentclassesdimid, bedlayersdimid 00081 00082 ! Drifters 00083 integer :: driftersdimid, driftersdimid2, drifterstimedimid, drifterstimevarid 00084 integer, allocatable :: driftersvarids(:) 00085 real*8 :: drift(3) 00086 00087 ! Ships 00088 integer :: shipdimid 00089 00090 ! Vegetation 00091 integer :: vegdimid 00092 00093 ! Q3D 00094 integer :: Q3Ddimid 00095 00096 ! global 00097 integer, dimension(:), allocatable :: globalvarids 00098 ! default output (fixed length) 00099 00100 ! points 00101 integer :: pointsdimid, pointnamelengthdimid 00102 integer :: xpointsvarid, ypointsvarid, pointtypesvarid, xpointindexvarid, ypointindexvarid, stationidvarid 00103 integer, dimension(:), allocatable :: pointsvarids 00104 integer, dimension(:),allocatable :: xpoints ! model x-coordinate of output points 00105 integer, dimension(:),allocatable :: ypoints ! model y-coordinate of output points 00106 00107 ! mean 00108 ! number of variables by number of parameters per variable (mean, sigma^2, min, max) 00109 integer, dimension(:,:), allocatable :: meanvarids 00110 character(slen), dimension(:), allocatable :: meanvartypes 00111 integer*4 :: nmeanvartypes = 4 ! number of time-average variable types 00112 integer*4,dimension(:),allocatable :: rugrowindex ! Array with row index where runup gauge can be found 00113 00114 00115 ! time 00116 integer :: globaltimedimid, pointtimedimid, meantimedimid 00117 integer :: globaltimevarid, pointtimevarid, meantimevarid 00118 00119 00120 ! TODO: check out why these are sometimes used.... 00121 integer :: tidetimedimid, windtimedimid 00122 integer :: inoutdimid, tidecornersdimid 00123 00124 ! local variables 00125 integer :: npointstotal 00126 integer :: itg,itp,itc,itm,itd 00127 ! Only alive at xmaster 00128 integer :: stpm ! size of tpm 00129 integer :: wordsize ! size of word in bytes 00130 00131 integer :: noutnumbers = 0 ! the number of outnumbers 00132 integer, dimension(numvars) :: outnumbers ! numbers, corrsponding to mnemonics, which are to be output 00133 00134 ! Output type 00135 integer :: NCREAL ! can be NF90_double or NF90_float 00136 00137 ! Fill numbers 00138 integer,parameter :: iFill = -huge(0) 00139 real,parameter :: sFill = -huge(0.0) 00140 real*8,parameter :: dFill = -dble(huge(0.0)) ! Robert: easier this way than to catch dFill<sFill later on 00141 contains 00142 00143 #ifdef USENETCDF 00144 ! Error handling of netcdf errors 00145 subroutine handle_err(status,file,line) 00146 use netcdf 00147 implicit none 00148 00149 integer, intent ( in) :: status 00150 character(*), intent(in) :: file 00151 integer, intent ( in) :: line 00152 integer :: status2 00153 00154 if(status /= nf90_noerr) then 00155 !UNIT=6 for stdout and UNIT=0 for stderr. 00156 write(0,'("NETCDF ERROR: ",a,i6,":",a)') file,line,trim(nf90_strerror(status)) 00157 write(0,*) 'closing file' 00158 status2 = nf90_close(ncid) 00159 if (status2 /= nf90_noerr) then 00160 write(0,*) 'NETCDF ERROR: ', __FILE__,__LINE__,trim(nf90_strerror(status2)) 00161 end if 00162 call halt_program 00163 end if 00164 end subroutine handle_err 00165 #endif 00166 00167 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00168 !!!!!!!!!!!!!!!!!!! INITIALISE OUTPUT !!!!!!!!!!!!!!!!!!!!!!!!!!!! 00169 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00170 #ifdef USENETCDF 00171 subroutine ncoutput_init(s, sl, par, tpar) 00172 use xmpi_module 00173 use params, only: parameters 00174 use spaceparams, only: indextos 00175 use spaceparamsdef 00176 use mnemmodule, only: arraytype, chartoindex 00177 use timestep_module, only: timepars 00178 use means_module, only: meanspars, meansparsglobal 00179 use logging_module, only: writelog 00180 ! This module is awaiting comments from Robert. 00181 use getkey_module, only: parameter, getkey, getkeys 00182 use paramsconst 00183 00184 00185 00186 implicit none 00187 00188 type(spacepars), intent(inout) :: s ! Use s => global data and sl => local data 00189 type(spacepars), intent(in) :: sl ! Use s => global data and sl => local data 00190 type(parameters), intent(in) :: par 00191 type(timepars), intent(in) :: tpar 00192 00193 ! Part of the getkey 00194 type(parameter) :: val 00195 type(arraytype) :: t 00196 type(meanspars) :: meanvar 00197 integer :: i,j 00198 integer :: rc ! return code 00199 character(slen) :: mnem 00200 00201 !integer :: npointstotal 00202 logical :: outputp, outputg, outputm 00203 integer, dimension(:), allocatable :: dimids ! store the dimids in a vector 00204 character(slen) :: coordinates 00205 character(slen) :: cellmethod 00206 00207 character(slen), dimension(:), allocatable :: keys 00208 logical :: dofortran, donetcdf 00209 ! subversion information 00210 include 'version.def' 00211 include 'version.dat' 00212 00213 00214 if (.not. xomaster) return 00215 00216 dofortran = par%outputformat .eq. OUTPUTFORMAT_FORTRAN .or. & 00217 & par%outputformat .eq. OUTPUTFORMAT_DEBUG 00218 donetcdf = par%outputformat .eq. OUTPUTFORMAT_NETCDF .or. & 00219 & par%outputformat .eq. OUTPUTFORMAT_DEBUG 00220 00221 outputp = .false. 00222 00223 ! initialize values 00224 00225 ! set output precision for NetCDF 00226 if(par%outputprecision == OUTPUTPRECISION_SINGLE) then 00227 NCREAL = NF90_REAL 00228 else 00229 NCREAL = NF90_DOUBLE 00230 endif 00231 00232 ! global 00233 00234 ! store netcdf variable ids for each variable 00235 allocate(globalvarids(par%nglobalvar)) 00236 globalvarids = -1 ! initialize to -1, so an error is raised when we miss something... 00237 outputg = .true. 00238 00239 npointstotal = par%npoints+par%nrugauge 00240 outputp = (npointstotal .gt. 0) .and. (size(tpar%tpp) .gt. 0) 00241 allocate(pointsvarids(par%npointvar)) 00242 00243 allocate(driftersvarids(par%ndrifter)) 00244 00245 allocate(meanvarids(par%nmeanvar,nmeanvartypes)) 00246 meanvarids = -1 00247 outputm = (par%nmeanvar .gt. 0) 00248 allocate(meanvartypes(nmeanvartypes)) 00249 meanvartypes = (/ 'mean ', 'var ', 'min ', 'max ' /) 00250 00251 ! create a file 00252 00253 NF90(nf90_create(path = par%ncfilename, cmode=ior(NF90_CLOBBER,NF90_64BIT_OFFSET), ncid = ncid)) 00254 00255 ! dimensions TODO: only output dimensions that are used 00256 ! grid 00257 NF90(nf90_def_dim(ncid, 'nx', s%nx+1, xdimid)) 00258 NF90(nf90_def_dim(ncid, 'ny', s%ny+1, ydimid)) 00259 00260 ! wave angles 00261 NF90(nf90_def_dim(ncid, 'wave_angle', s%ntheta, thetadimid)) 00262 00263 ! computational layers in bed ... 00264 NF90(nf90_def_dim(ncid, 'bed_layers', par%nd, bedlayersdimid)) 00265 00266 ! sediment classes 00267 NF90(nf90_def_dim(ncid, 'sediment_classes', par%ngd, sedimentclassesdimid)) 00268 00269 ! dimensions of length 2.... what is this.... TODO: find out what this is 00270 NF90(nf90_def_dim(ncid, 'inout', 2, inoutdimid)) 00271 00272 ! write(*,*) 'Writing ndrifter', par%ndrifter 00273 if (par%ndrifter .gt. 0) then 00274 ! create dimensions for drifters and drifterstime: 00275 NF90(nf90_def_dim(ncid, 'ndrifter', par%ndrifter, driftersdimid)) 00276 NF90(nf90_def_dim(ncid, 'drifterstime', size(tpar%tpp), drifterstimedimid)) 00277 end if 00278 00279 ! write(*,*) 'Writing nship', par%nship 00280 if (par%nship .gt. 0) then 00281 NF90(nf90_def_dim(ncid, 'nship', par%nship, shipdimid)) 00282 end if 00283 00284 ! vegetation sections 00285 if (par%vegetation .gt. 0) then 00286 NF90(nf90_def_dim(ncid, 'nsecvegmax', s%nsecvegmax, vegdimid)) 00287 endif 00288 00289 ! write(*,*) 'Writing nz', par%nz 00290 if (par%nz .gt. 1) then 00291 NF90(nf90_def_dim(ncid, 'nz', par%nz, Q3Ddimid)) 00292 end if 00293 00294 ! time dimensions are fixed, only defined if there are points 00295 if (outputg) then 00296 NF90(nf90_def_dim(ncid, 'globaltime', NF90_unlimited, globaltimedimid)) 00297 end if 00298 if (outputp) then 00299 ! points 00300 NF90(nf90_def_dim(ncid, 'points', npointstotal, pointsdimid)) 00301 NF90(nf90_def_dim(ncid, 'pointtime', size(tpar%tpp), pointtimedimid)) 00302 NF90(nf90_def_dim(ncid, 'pointnamelength', 64, pointnamelengthdimid)) 00303 end if 00304 if (outputm) then 00305 NF90(nf90_def_dim(ncid, 'meantime', size(tpar%tpm)-1, meantimedimid)) 00306 end if 00307 00308 if (s%tidelen > 0) then 00309 NF90(nf90_def_dim(ncid, 'tidetime', s%tidelen, tidetimedimid)) 00310 endif 00311 00312 if (par%tideloc > 0) then 00313 NF90(nf90_def_dim(ncid, 'tidecorners', par%tideloc, tidecornersdimid)) 00314 endif 00315 00316 if (s%windlen > 0) then 00317 NF90(nf90_def_dim(ncid, 'windtime', s%windlen, windtimedimid)) 00318 endif 00319 00320 ! define empty parameter variable 00321 NF90(nf90_def_var(ncid, '_parameters', NCREAL, varid=parvarid)) 00322 ! define space & time variables 00323 ! grid 00324 NF90(nf90_def_var(ncid, 'globalx', NCREAL, (/ xdimid, ydimid /), xvarid)) 00325 NF90(nf90_put_att(ncid, xvarid, 'units', 'm')) 00326 NF90(nf90_put_att(ncid, xvarid, 'long_name', 'local x coordinate')) 00327 NF90(nf90_put_att(ncid, xvarid, 'standard_name', 'projection_x_coordinate')) 00328 NF90(nf90_put_att(ncid, xvarid, 'axis', 'X')) 00329 ! For compatibility with CSIRO Dive software 00330 if (len(trim(par%projection)) .ne. 0) then 00331 NF90(nf90_put_att(ncid, xvarid, 'projection', par%projection)) 00332 NF90(nf90_put_att(ncid, xvarid, 'rotation',( s%alfa/atan(1.0d0)*45.d0))) 00333 end if 00334 00335 NF90(nf90_def_var(ncid, 'globaly', NCREAL, (/ xdimid, ydimid /), yvarid)) 00336 NF90(nf90_put_att(ncid, yvarid, 'units', 'm')) 00337 NF90(nf90_put_att(ncid, yvarid, 'long_name', 'local y coordinate')) 00338 NF90(nf90_put_att(ncid, yvarid, 'standard_name', 'projection_y_coordinate')) 00339 NF90(nf90_put_att(ncid, yvarid, 'axis', 'Y')) 00340 ! For compatibility with CSIRO Dive software 00341 if (len(trim(par%projection)) .ne. 0) then 00342 NF90(nf90_put_att(ncid, yvarid, 'projection', par%projection)) 00343 NF90(nf90_put_att(ncid, yvarid, 'rotation',( s%alfa/atan(1.0d0)*45.d0))) 00344 end if 00345 00346 ! Some metadata attributes 00347 NF90(nf90_put_att(ncid,nf90_global, "Conventions", "CF-1.4")) 00348 NF90(nf90_put_att(ncid,nf90_global, "Producer", "XBeach littoral zone wave model (http://www.xbeach.org)")) 00349 NF90(nf90_put_att(ncid,nf90_global, "Build-Revision", trim(Build_Revision))) 00350 NF90(nf90_put_att(ncid,nf90_global, "Build-Date", trim(Build_Date))) 00351 NF90(nf90_put_att(ncid,nf90_global, "URL", trim(Build_URL))) 00352 00353 ! Store all the parameters 00354 ! This part is awaiting comments from Robert McCall 00355 call getkeys(par, keys) 00356 do i=1,size(keys) 00357 rc = getkey(par, keys(i), val) 00358 if (val%type == 'i') then 00359 NF90(nf90_put_att(ncid, parvarid, keys(i), val%i0 )) 00360 elseif (val%type == 'c') then 00361 NF90(nf90_put_att(ncid, parvarid, keys(i), val%c0 )) 00362 elseif (val%type == 'r') then 00363 NF90(nf90_put_att(ncid, parvarid, keys(i), val%r0 )) 00364 end if 00365 end do 00366 00367 ! global 00368 if (outputg) then 00369 NF90(nf90_def_var(ncid, 'globaltime', NCREAL, (/ globaltimedimid /), globaltimevarid)) 00370 NF90(nf90_put_att(ncid, globaltimevarid, 'units', trim(par%tunits))) 00371 NF90(nf90_put_att(ncid, globaltimevarid, 'axis', 'T')) 00372 NF90(nf90_put_att(ncid, globaltimevarid, 'standard_name', 'time')) 00373 00374 00375 ! default global output variables 00376 do i=1,par%nglobalvar 00377 mnem = par%globalvars(i) 00378 coordinates = '' 00379 j = chartoindex(mnem) 00380 call indextos(s,j,t) 00381 call writelog('ls', '', 'Creating netcdf variable: ', trim(mnem) ) 00382 00383 ! Build the array with dimension ids 00384 allocate(dimids(t%rank+1)) 00385 select case(t%rank) 00386 case(0) 00387 dimids = (/ globaltimedimid /) 00388 coordinates = '' 00389 case(1) 00390 dimids = (/ dimensionid(t%dimensions(1)), globaltimedimid /) 00391 coordinates = '' 00392 if (dimids(1) .eq. xdimid) coordinates = 'globalx' 00393 if (dimids(1) .eq. ydimid) coordinates = 'globaly' 00394 case(2) 00395 dimids = (/ dimensionid(t%dimensions(1)), dimensionid(t%dimensions(2)), globaltimedimid /) 00396 coordinates = 'globalx globaly' 00397 case(3) 00398 dimids = (/ dimensionid(t%dimensions(1)), dimensionid(t%dimensions(2)), & 00399 &dimensionid(t%dimensions(3)), globaltimedimid /) 00400 coordinates = 'globalx globaly' 00401 ! Do we have a vertical level? 00402 if (dimids(3) .eq. bedlayersdimid) coordinates = trim(coordinates) // ' bed_layers' 00403 case(4) 00404 call writelog('ls', '', 'Variable ' // trim(mnem) // ' is of rank 4. This may not work due to an' // & 00405 &' unresolved issue. If so, remove the variable or use the fortran outputformat option.') 00406 dimids = (/ dimensionid(t%dimensions(1)), dimensionid(t%dimensions(2)), & 00407 &dimensionid(t%dimensions(3)), dimensionid(t%dimensions(4)), globaltimedimid /) 00408 coordinates = 'globalx globaly' 00409 ! Do we have a vertical level? 00410 if ((dimids(3) .eq. bedlayersdimid) .or. (dimids(4) .eq. bedlayersdimid)) then 00411 coordinates = trim(coordinates) // ' bed_layers' 00412 end if 00413 case default 00414 call writelog('lse', '', 'mnem: ' // mnem // ' not supported, rank:', t%rank) 00415 stop 1 00416 end select 00417 select case(t%type) 00418 case('i') 00419 NF90(nf90_def_var(ncid, trim(mnem), NF90_INT, dimids, globalvarids(i))) 00420 case('r') 00421 NF90(nf90_def_var(ncid, trim(mnem), NCREAL, dimids, globalvarids(i))) 00422 case default 00423 write(0,*) 'mnem', mnem, ' not supported, type:', t%type 00424 end select 00425 deallocate(dimids) 00426 NF90(nf90_put_att(ncid, globalvarids(i), 'coordinates', trim(coordinates))) 00427 NF90(nf90_put_att(ncid, globalvarids(i), 'units', trim(t%units))) 00428 if (.not.(trim(t%standardname) .eq. '')) then 00429 NF90(nf90_put_att(ncid, globalvarids(i), 'standard_name', trim(t%standardname))) 00430 endif 00431 NF90(nf90_put_att(ncid, globalvarids(i), 'long_name', trim(t%description))) 00432 select case(t%type) 00433 case('i') 00434 NF90(nf90_put_att(ncid, globalvarids(i), '_FillValue', iFill)) 00435 case('r') 00436 if(par%outputprecision == OUTPUTPRECISION_SINGLE) then 00437 NF90(nf90_put_att(ncid, globalvarids(i), '_FillValue', sFill)) 00438 else 00439 NF90(nf90_put_att(ncid, globalvarids(i), '_FillValue', dFill)) 00440 endif 00441 end select 00442 end do 00443 end if 00444 00445 ! ! points 00446 ! default global output variables 00447 if (outputp) then 00448 NF90(nf90_def_var(ncid, 'pointtime', NCREAL, (/ pointtimedimid /), pointtimevarid)) 00449 NF90(nf90_put_att(ncid, pointtimevarid, 'units', trim(par%tunits))) 00450 NF90(nf90_put_att(ncid, pointtimevarid, 'axis', 'T')) 00451 NF90(nf90_put_att(ncid, pointtimevarid, 'standard_name', 'time')) 00452 00453 ! points 00454 NF90(nf90_def_var(ncid, 'pointx', NCREAL, (/ pointsdimid /), xpointsvarid)) 00455 NF90(nf90_put_att(ncid, xpointsvarid, 'units', 'm')) 00456 NF90(nf90_put_att(ncid, xpointsvarid, 'long_name', 'local x coordinate')) 00457 NF90(nf90_put_att(ncid, xpointsvarid, 'standard_name', 'projection_x_coordinate')) 00458 NF90(nf90_put_att(ncid, xpointsvarid, 'axis', 'X')) 00459 00460 00461 NF90(nf90_def_var(ncid, 'pointy', NCREAL, (/ pointsdimid /), ypointsvarid)) 00462 NF90(nf90_put_att(ncid, ypointsvarid, 'units', 'm')) 00463 NF90(nf90_put_att(ncid, ypointsvarid, 'long_name', 'local y coordinate')) 00464 NF90(nf90_put_att(ncid, ypointsvarid, 'standard_name', 'projection_y_coordinate')) 00465 NF90(nf90_put_att(ncid, ypointsvarid, 'axis', 'Y')) 00466 00467 NF90(nf90_def_var(ncid, 'station_id', NF90_CHAR, (/ pointnamelengthdimid, pointsdimid /), stationidvarid)) 00468 NF90(nf90_put_att(ncid, stationidvarid, 'long_name', 'station identification code')) 00469 NF90(nf90_put_att(ncid, stationidvarid, 'standard_name', 'station_id')) 00470 00471 NF90(nf90_def_var(ncid, 'xpointindex', NF90_INT, (/ pointsdimid /), xpointindexvarid)) 00472 ! wwvv above was NF90_DOUBLE 00473 NF90(nf90_put_att(ncid, xpointindexvarid, 'long_name', 'nearest x grid cell')) 00474 NF90(nf90_def_var(ncid, 'ypointindex', NF90_INT, (/ pointsdimid /), ypointindexvarid)) 00475 ! wwvv above was NF90_DOUBLE 00476 NF90(nf90_put_att(ncid, ypointindexvarid, 'long_name', 'nearest y grid cell')) 00477 00478 NF90(nf90_def_var(ncid, 'pointtypes', NF90_INT, (/ pointsdimid /), pointtypesvarid)) 00479 ! wwvv above was NF90_DOUBLE 00480 NF90(nf90_put_att(ncid, pointtypesvarid, 'long_name', 'type of point (0=point, 1=rugauge)')) 00481 00482 do i=1,par%npointvar 00483 mnem = par%pointvars(i) 00484 j = chartoindex(mnem) 00485 coordinates = '' 00486 call indextos(s,j,t) 00487 select case(t%type) 00488 case('r') 00489 ! Build the array with dimension ids 00490 call writelog('ls', '', 'Creating netcdf variable: ', 'point_'// trim(mnem) ) 00491 allocate(dimids(t%rank)) 00492 ! Make sure the variable has x and y as the first 2 dimensions 00493 if ((dimensionid(t%dimensions(1)) .ne. xdimid) .or. (dimensionid(t%dimensions(2)) .ne. ydimid)) then 00494 call writelog('lse','', 'Tried to store variable ' // trim(mnem) // ', but it is not a function of x,y') 00495 endif 00496 select case(t%rank) 00497 case(2) 00498 dimids = (/ pointsdimid, pointtimedimid /) 00499 coordinates = 'pointx pointy' 00500 case(3) 00501 dimids = (/ pointsdimid, dimensionid(t%dimensions(3)), pointtimedimid /) 00502 coordinates = 'pointx pointy' 00503 ! Do we have a vertical level? 00504 if (dimids(3) .eq. bedlayersdimid) coordinates = trim(coordinates) // ' bed_layers' 00505 case(4) 00506 call writelog('ls', '', 'Variable ' // trim(mnem) // ' is of rank 4. This may not work due to an' // & 00507 &' unresolved issue. If so, remove the variable or use the fortran outputformat option.') 00508 dimids = (/ pointsdimid, dimensionid(t%dimensions(3)), dimensionid(t%dimensions(4)), pointtimedimid /) 00509 coordinates = 'pointx pointy' 00510 ! Do we have a vertical level? 00511 if ((dimids(3) .eq. bedlayersdimid) .or. (dimids(4) .eq. bedlayersdimid)) then 00512 coordinates = trim(coordinates) // ' bed_layers' 00513 end if 00514 case default 00515 call writelog('lse', '', 'mnem: ' // mnem // ' not supported, rank:', t%rank) 00516 stop 1 00517 end select 00518 NF90(nf90_def_var(ncid, 'point_' // trim(mnem), NCREAL, dimids, pointsvarids(i))) 00519 NF90(nf90_put_att(ncid, pointsvarids(i), 'coordinates', trim(coordinates))) 00520 deallocate(dimids) 00521 case default 00522 write(0,*) 'mnem', mnem, ' not supported, type:', t%type 00523 end select 00524 NF90(nf90_put_att(ncid, pointsvarids(i), 'units', trim(t%units))) 00525 if (.not.(trim(t%standardname) .eq. '')) then 00526 NF90(nf90_put_att(ncid, pointsvarids(i), 'standard_name', trim(t%standardname))) 00527 endif 00528 NF90(nf90_put_att(ncid, pointsvarids(i), 'long_name', trim(t%description))) 00529 end do 00530 endif ! outputp 00531 00532 if (par%ndrifter .gt. 0) then 00533 allocate(dimids(2)) 00534 00535 NF90(nf90_def_var(ncid, 'driftertime', NCREAL, (/ drifterstimedimid /), drifterstimevarid)) 00536 NF90(nf90_put_att(ncid, drifterstimevarid, 'units', trim(par%tunits))) 00537 NF90(nf90_put_att(ncid, drifterstimevarid, 'axis', 'T')) 00538 NF90(nf90_put_att(ncid, drifterstimevarid, 'standard_name', 'time')) 00539 00540 ! create netcdf variables for drifters: 00541 ! each variable is a 2-d array, containing the x-y values 00542 ! the netcdf name of e.g. the 3rd array is drifter_003 00543 00544 NF90(nf90_def_dim(ncid, 'drifterstime2', 2, driftersdimid2)) 00545 do i=1,par%ndrifter 00546 ! we need the following dimensions per drifter: 00547 ! driftersdimid2 : 2 00548 ! drifterstimedimid : number of output time steps 00549 dimids(1) = driftersdimid2 00550 dimids(2) = drifterstimedimid 00551 write(mnem,'("drifter_",I0.3)') i 00552 NF90(nf90_def_var(ncid, trim(mnem), NCREAL, dimids, driftersvarids(i))) 00553 NF90(nf90_put_att(ncid, driftersvarids(i), 'coordinates', 'pointx pointy')) 00554 NF90(nf90_put_att(ncid, driftersvarids(i), 'units', 'm')) 00555 enddo 00556 deallocate(dimids) 00557 00558 endif ! par%ndrifter 00559 00560 if(outputm) then 00561 NF90(nf90_def_var(ncid, 'meantime', NCREAL, (/ meantimedimid /), meantimevarid)) 00562 NF90(nf90_put_att(ncid, meantimevarid, 'units', trim(par%tunits))) 00563 NF90(nf90_put_att(ncid, meantimevarid, 'axis', 'T')) 00564 NF90(nf90_put_att(ncid, meantimevarid, 'standard_name', 'time')) 00565 ! default global output variables 00566 do i=1,par%nmeanvar 00567 ! Not sure if this is required here, but it is used in varoutput 00568 ! #ifdef USEMPI 00569 ! ! No need to collect here, we're just using the types 00570 ! ! call means_collect(sl,meansparsglobal(i),meansparslocal(i)) 00571 ! #else 00572 ! meansparsglobal(i)=meansparslocal(i) 00573 ! #endif 00574 coordinates = '' 00575 meanvar = meansparsglobal(i) 00576 t = meanvar%t 00577 select case(t%type) 00578 case('r') 00579 ! Build the array with dimension ids 00580 allocate(dimids(t%rank+1)) 00581 select case(t%rank) 00582 case(2) 00583 dimids = (/ dimensionid(t%dimensions(1)), dimensionid(t%dimensions(2)), meantimedimid /) 00584 coordinates = 'globalx globaly' 00585 case(3) 00586 dimids = (/ dimensionid(t%dimensions(1)), dimensionid(t%dimensions(2)), & 00587 &dimensionid(t%dimensions(3)), meantimedimid /) 00588 coordinates = 'globalx globaly' 00589 ! Do we have a vertical level? 00590 if (dimids(3) .eq. bedlayersdimid) coordinates = trim(coordinates) // ' bed_layers' 00591 case(4) 00592 call writelog('ls', '', 'Variable ' // trim(mnem) // ' is of rank 4. This may not work due to an' // & 00593 &' unresolved issue. If so, remove the variable or use the fortran outputformat option.') 00594 dimids = (/ dimensionid(t%dimensions(1)), dimensionid(t%dimensions(2)), & 00595 &dimensionid(t%dimensions(3)), dimensionid(t%dimensions(4)), meantimedimid /) 00596 coordinates = 'globalx globaly' 00597 ! Do we have a vertical level? 00598 if ((dimids(3) .eq. bedlayersdimid) .or. (dimids(4) .eq. bedlayersdimid)) then 00599 coordinates = trim(coordinates) // ' bed_layers' 00600 end if 00601 case default 00602 call writelog('lse', '', 'mnem: ' // mnem // ' not supported, rank:', t%rank) 00603 call halt_program 00604 stop 1 00605 end select 00606 00607 ! Create a variable for all types of meanvars (mean, var, min, max) 00608 do j = 1,nmeanvartypes 00609 cellmethod = meanvartypes(j) 00610 call writelog('ls', '', 'Creating netcdf variable: ', trim(t%name) // '_' // cellmethod) 00611 NF90(nf90_def_var(ncid, trim(t%name) // '_' // trim(cellmethod), NCREAL, dimids, meanvarids(i,j))) 00612 NF90(nf90_put_att(ncid, meanvarids(i,j), 'coordinates', trim(coordinates))) 00613 if (cellmethod .eq. 'var') then 00614 NF90(nf90_put_att(ncid, meanvarids(i,j), 'units', '(' // trim(t%units) // ')^2')) 00615 else 00616 NF90(nf90_put_att(ncid, meanvarids(i,j), 'units', trim(t%units))) 00617 endif 00618 if (.not.(trim(t%standardname) .eq. '')) then 00619 NF90(nf90_put_att(ncid, meanvarids(i,j), 'standard_name', trim(t%standardname))) 00620 endif 00621 NF90(nf90_put_att(ncid, meanvarids(i,j), 'long_name', trim(t%description))) 00622 ! For H and urms we don't compute the mean but the rms of the rms..... 00623 if (cellmethod .eq. 'mean' .and. ((t%name .eq. 'H') .or. (t%name .eq. 'urms'))) then 00624 cellmethod = 'rms' 00625 end if 00626 NF90(nf90_put_att(ncid, meanvarids(i,j), 'cell_methods', 'meantime: ' // trim(cellmethod))) 00627 end do 00628 00629 deallocate(dimids) 00630 case default 00631 write(0,*) 'mnem', mnem, ' not supported, type:', t%type 00632 end select 00633 end do 00634 endif ! outputm 00635 00636 ! done defining variables 00637 call writelog('ls', '', 'Writing file definition.') 00638 NF90(nf90_enddef(ncid)) 00639 00640 ! Fill meta variables 00641 ! Grid 00642 j = chartoindex('xz') 00643 call indextos(s,j,t) 00644 NF90(nf90_put_var(ncid, xvarid, CONVREAL(t%r2))) 00645 00646 j = chartoindex('yz') 00647 call indextos(s,j,t) 00648 NF90(nf90_put_var(ncid, yvarid, CONVREAL(t%r2))) 00649 00650 if (outputp) then 00651 call writelog('ls', '', 'Writing point vars.') 00652 NF90(nf90_put_var(ncid, xpointsvarid, CONVREAL(par%xpointsw))) 00653 NF90(nf90_put_var(ncid, ypointsvarid, CONVREAL(par%ypointsw))) 00654 NF90(nf90_put_var(ncid, pointtypesvarid, par%pointtypes)) 00655 NF90(nf90_put_var(ncid, xpointindexvarid, xpoints)) 00656 NF90(nf90_put_var(ncid, ypointindexvarid, ypoints)) 00657 NF90(nf90_put_var(ncid, pointtypesvarid, par%pointtypes)) 00658 NF90(nf90_put_var(ncid, stationidvarid, par%stationid(1:npointstotal))) 00659 end if 00660 00661 NF90(nf90_close(ncid)) 00662 ! wwvv sl is not used so it should be removed. 00663 ! wwvv for now, to avoid warning: 00664 00665 if (sl%nx .ne. -1) return 00666 end subroutine ncoutput_init 00667 #endif 00668 ! USENETCDF 00669 00670 00671 !___________________________________________________________________________________ 00672 00673 subroutine ncoutput(s,sl,par, tpar) 00674 use xmpi_module 00675 use params, only: parameters 00676 use spaceparamsdef, only: spacepars 00677 use spaceparams, only: indextos 00678 use timestep_module, only: timepars 00679 use mnemmodule, only: arraytype, mnem_Subg, mnem_Sutot, mnem_Svtot, maxnamelen, mnem_alfaz, mnem_Susg 00680 use mnemmodule, only: mnem_Svbg, mnem_Svsg, chartoindex, mnem_wetz 00681 use means_module, only: meansparsglobal, meansparslocal 00682 use postprocessmod, only: get_sister_mnem, gridrotate, postprocessvar_r2 00683 use spaceparamsdef, only: spacepars 00684 use xmpi_module 00685 use constants, only: pi 00686 use paramsconst 00687 #ifdef USEMPI 00688 use xmpi_module 00689 00690 use spaceparams, only: space_collect_mnem, space_local_to_global, space_collect_index 00691 use spaceparams, only: space_who_has, space_global_to_local 00692 use means_module, only: means_collect 00693 #endif 00694 00695 implicit none 00696 00697 type(spacepars), intent(inout) :: s ! s-> spaceparams, what is isl? 00698 ! s describes the global system, sl the local system 00699 ! in this MPI process 00700 type(spacepars), intent(inout) :: sl ! s-> spaceparams, what is isl? 00701 type(parameters), intent(inout) :: par 00702 type(timepars), intent(in) :: tpar 00703 00704 type(arraytype) :: t 00705 integer :: i,j,ii 00706 #ifdef USEMPI 00707 integer :: index 00708 #endif 00709 character(maxnamelen) :: mnem,sistermnemalloc 00710 real*8, dimension(:,:), allocatable :: points 00711 00712 ! some local variables to pass the data through the postprocessing function. 00713 integer :: i0 00714 integer, dimension(:,:), allocatable :: i2 00715 integer, dimension(:,:,:), allocatable :: i3 00716 CONVREALTYPE :: r0conv 00717 real*8, dimension(:), allocatable :: r1 00718 CONVREALTYPE, dimension(:), allocatable :: r1conv 00719 real*8, dimension(:,:), allocatable :: r2 00720 CONVREALTYPE, dimension(:,:), allocatable :: r2conv 00721 real*8, dimension(:,:,:), allocatable :: r3 00722 CONVREALTYPE, dimension(:,:,:), allocatable :: r3conv 00723 real*8, dimension(:,:,:,:), allocatable :: r4 00724 CONVREALTYPE, dimension(:,:,:,:), allocatable :: r4conv 00725 real*8, allocatable :: tempvectorr(:) 00726 real*8,dimension(size(tpar%tpg)+size(tpar%tpp)+size(tpar%tpm)) :: outputtimes 00727 00728 integer :: jtg,reclen,unit,idumhl,ird,iru,xmax,xmin 00729 integer :: iz,jz,idum 00730 real*8 :: di,dj,dx,dy 00731 integer :: ilocal, jlocal 00732 #ifdef USEMPI 00733 real*8, dimension(par%ndrifter) :: idriftlocal,jdriftlocal 00734 #endif 00735 integer :: pii 00736 00737 #ifdef USENETCDF 00738 integer :: status 00739 #endif 00740 00741 logical :: dofortran, donetcdf, dofortran_compat 00742 logical :: dooutput_global, dooutput_mean, dooutput_point, dooutput_drifter, dooutput_hotstart 00743 00744 type pointoutput 00745 integer :: rank ! rank of the data 00746 character(len=maxnamelen) :: name ! name of the variable 00747 real*8, dimension(:,:), allocatable :: r2 ! contains 2-d data (rank = 2) 00748 real*8, dimension(:,:,:), allocatable :: r3 ! contains 3-d data (rank = 3) 00749 real*8, dimension(:,:,:,:), allocatable :: r4 ! contains 4-d data (rank = 4) 00750 ! note that r2 will be allocated (1,1) 00751 ! r3 will be allocated (1,1,nnn) 00752 ! r4 will be allocated (1,1,mmm,nnn) 00753 ! why ? because gridrotate demands this 00754 end type pointoutput 00755 00756 integer :: xpii, ypii 00757 integer :: rugx, rugy 00758 00759 type (pointoutput), dimension(:,:), allocatable :: pointoutputs 00760 00761 real*8, dimension(:,:), allocatable :: runups,runups1 ! runups(:,i) will contain the 1+par%nrugdepth*3 runup values 00762 ! ! for runup number i(i=1 .. par%nrugauge) 00763 integer, dimension(:), allocatable:: xpoints1 ! for netcdf runup output 00764 00765 ! fortran output requested? 00766 dofortran = par%outputformat .eq. OUTPUTFORMAT_FORTRAN .or. & 00767 & par%outputformat .eq. OUTPUTFORMAT_DEBUG 00768 00769 ! for compatibility with older version, at some places 00770 ! fortran output is only done when tpar%output = .true. 00771 00772 dofortran_compat = dofortran .and. tpar%output 00773 00774 !netcdf output requested? 00775 donetcdf = par%outputformat .eq. OUTPUTFORMAT_NETCDF .or. & 00776 & par%outputformat .eq. OUTPUTFORMAT_DEBUG 00777 00778 ! time for global output? 00779 dooutput_global = tpar%outputg .and. par%nglobalvar .gt. 0 00780 00781 ! time for mean output? 00782 dooutput_mean = tpar%outputm .and. par%nmeanvar .gt. 0 .and. tpar%itm .gt. 1 00783 00784 ! time for point output? 00785 dooutput_point = tpar%outputp .and. par%npointvar .gt. 0 & 00786 .and. par%npoints+par%nrugauge .gt. 0 00787 00788 ! time for drifter output? 00789 dooutput_drifter = tpar%outputp .and. par%ndrifter .gt. 0 00790 00791 ! time for hotstart output? 00792 dooutput_hotstart = tpar%outputh .and. par%writehotstart == 1 00793 00794 #ifdef USEMPI 00795 ! clear collected items 00796 s%collected = s%precollected 00797 00798 ! If we're gonna write some global output 00799 if (dooutput_global) then 00800 ! we'll need to collect the information from all nodes. 00801 do i=1,par%nglobalvar 00802 mnem = par%globalvars(i) 00803 index = chartoindex(mnem) 00804 call space_collect_index(s,sl,par,index) 00805 if (par%rotate==1) then 00806 sistermnemalloc = get_sister_mnem(mnem) 00807 00808 select case (sistermnemalloc) 00809 case ('none') 00810 ! nothing 00811 case default 00812 call space_collect_mnem(s,sl,par,sistermnemalloc) 00813 end select 00814 00815 select case(mnem) 00816 case(mnem_Sutot,mnem_Svtot) 00817 call space_collect_mnem(s,sl,par,mnem_Subg) 00818 call space_collect_mnem(s,sl,par,mnem_Svbg) 00819 call space_collect_mnem(s,sl,par,mnem_Susg) 00820 call space_collect_mnem(s,sl,par,mnem_Svsg) 00821 end select 00822 00823 endif 00824 end do 00825 if (par%rotate==1) then 00826 call space_collect_mnem(s,sl,par,mnem_alfaz) 00827 endif 00828 if(par%remdryoutput==1) then 00829 call space_collect_mnem(s,sl,par,mnem_wetz) 00830 endif 00831 00832 endif 00833 00834 #endif 00835 ! USEMPI 00836 00837 00838 #ifdef USEMPI 00839 if(dooutput_point) then 00840 ! Collect all data for which we store the points. 00841 ! TODO: This will be a lot faster if nodes write their own point. Use the parallel netcdf for that. 00842 ! Let's wait till someone needs it... 00843 ! wwvv Collecting the data could als be done using sends from the compute processes 00844 ! wwvv and receives on the xomaster process. 00845 ! wwvv the info to do this properly could be computed in _init 00846 ! wwvv in stead of sends and receives, probably more simple is the use of mpi_alltoallw 00847 ! 00848 endif !dooutput_point 00849 #endif 00850 ! USEMPI 00851 00852 ! If we're gonna write some mean output 00853 if(dooutput_mean) then 00854 00855 do i=1,par%nmeanvar 00856 #ifdef USEMPI 00857 call means_collect(sl,meansparsglobal(i),meansparslocal(i)) 00858 #else 00859 meansparsglobal(i)=meansparslocal(i) 00860 #endif 00861 end do 00862 endif 00863 00864 if (dooutput_hotstart) then 00865 #ifdef USEMPI 00866 call hotstartfiles_collect(s,sl,par, tpar) 00867 #endif 00868 endif 00869 00870 00871 #ifdef USEMPI 00872 if (par%ndrifter .gt. 0) then 00873 if(xmaster) then 00874 idriftlocal = sl%idrift 00875 jdriftlocal = sl%jdrift 00876 endif 00877 call xmpi_send(xmpi_imaster,xmpi_omaster,idriftlocal) 00878 call xmpi_send(xmpi_imaster,xmpi_omaster,jdriftlocal) 00879 if (xomaster) then 00880 s%idrift = idriftlocal 00881 s%jdrift = jdriftlocal 00882 endif 00883 endif 00884 #endif 00885 00886 ! wwvv a more efficient implementation: 00887 ! step 1: determine which process contains the desired point 00888 ! step 2: gridrotate this point 00889 ! step 3: store this point 00890 ! step 4: send the value of this point to xomaster (xmpi_send) after 00891 ! if( .not. xomaster) return 00892 ! 00893 ! the points are assembled in pointoutputs, which is an 2-d array 00894 ! of type(pointoutput), which is described above. 00895 ! 00896 ! pointoutputs(i,j) contains the values for pointvariable i 00897 ! and coordinates j 00898 00899 if (dooutput_point) then 00900 if (xomaster) then 00901 allocate(pointoutputs(par%npointvar,par%npoints+par%nrugauge)) 00902 endif 00903 allocate(tempvectorr(1+par%nrugdepth*3)) 00904 allocate(runups (size(tempvectorr),par%nrugauge)) 00905 allocate(runups1(size(tempvectorr),par%nrugauge)) 00906 ! WD: new code 00907 ! first: compute run gauge 00908 if(xcompute) then 00909 allocate(xpoints1(par%nrugauge)) 00910 xpoints1 = s%nx+1 00911 ! the compute processes will determine their own runup values 00912 ! if no runp is found, a value of huge(0.0d0) is used 00913 ! The runup values are collected in xomaster, using xmpi_reduce, 00914 ! taking the minimum values 00915 do ii = par%npoints+1, par%npoints+par%nrugauge 00916 iru = ii-par%npoints ! iru will run from 1 to par%nrugauge 00917 tempvectorr = huge(0.d0) 00918 ! for netcdf output: 00919 do ird=1,par%nrugdepth 00920 #ifdef USEMPI 00921 call space_global_to_local(sl,1,rugrowindex(iru),rugx,rugy) 00922 #else 00923 rugx = 1 00924 rugy = rugrowindex(iru) 00925 #endif 00926 ! rugy now contains the local y-coordinate 00927 ! rugx is not used 00928 ! if rugy is within the computational domain of this process, 00929 ! we determine the local runup values 00930 ! the computational domain in this process: 00931 ! 1st dimension: sl%icls(xmpi_rank+1) .. sl%icle(xmpi_rank+1) 00932 ! 2nd dimension: sl%jcls(xmpi_rank+1) .. sl%jcle(xmpi_rank+1) 00933 #ifdef USEMPI 00934 if (rugy .ge. sl%jcls(xmpi_rank+1) .and. rugy .le. sl%jcle(xmpi_rank+1)) then 00935 xmin = sl%icls(xmpi_rank+1) 00936 xmax = sl%icle(xmpi_rank+1) 00937 #else 00938 if (.true.) then 00939 xmin = 2 00940 xmax = sl%nx+1 00941 #endif 00942 idumhl = -1 ! Set default 00943 if (rugrowindex(iru)>0) then ! master domain always contains this runup gauge 00944 ! local index of minimum location where hh<rugdepth 00945 do j=xmin+1,xmax 00946 if ((sl%hh(j, rugy)<=par%rugdepth(ird)) .and. & 00947 & (sl%hh(j-1,rugy)> par%rugdepth(ird)) ) then 00948 idumhl=j-1 00949 exit 00950 endif 00951 enddo 00952 endif 00953 00954 if (par%morfacopt==1) then 00955 tempvectorr(1) = par%t*max(par%morfac,1.d0) 00956 else 00957 tempvectorr(1)=par%t 00958 endif 00959 00960 ! for netcdf output: 00961 if(ird == 1) then 00962 ! convert from local coordinates to global, y-coordinate is not relevant: 00963 if (idumhl .gt. 0) then 00964 #ifdef USEMPI 00965 call space_local_to_global(sl,idumhl,10,xpoints1(iru),idum) 00966 #else 00967 xpoints1(iru) = idumhl 00968 idum = 10 00969 #endif 00970 else 00971 xpoints1(iru) = s%nx+1 00972 endif 00973 endif 00974 00975 if (idumhl .gt. 0) then 00976 tempvectorr((ird-1)*3+2) = sl%xz(idumhl,rugy) 00977 tempvectorr((ird-1)*3+3) = sl%yz(idumhl,rugy) 00978 tempvectorr((ird-1)*3+4) = sl%zs(idumhl,rugy) 00979 endif 00980 00981 endif ! rugy within computational domain 00982 enddo ! ird=1,par%nrugdepth 00983 runups(:,iru) = tempvectorr 00984 enddo ! ii = par%npoints+1, par%npoints+par%nrugauge 00985 ! reduce the runups to runup1 on xmaster: 00986 #ifdef USEMPI 00987 call xmpi_reduce(runups,runups1,MPI_MIN) 00988 #else 00989 runups1 = runups 00990 #endif 00991 00992 ! reduce xpoints1 to (par%npoints+1:) on xmaster: 00993 #ifdef USEMPI 00994 call xmpi_reduce(xpoints1,xpoints(par%npoints+1:),MPI_MIN) 00995 #else 00996 xpoints(par%npoints+1:) = xpoints1 00997 #endif 00998 endif ! xcompute 00999 ! send runups1 to xomaster who will receive it in runups 01000 #ifdef USEMPI 01001 if (xomaster) then 01002 call xmpi_send(xmpi_imaster,xmpi_omaster,runups) 01003 else 01004 call xmpi_send(xmpi_imaster,xmpi_omaster,runups1) 01005 endif 01006 ! same for xpoints: 01007 call xmpi_send(xmpi_imaster,xmpi_omaster,xpoints(par%npoints+1:)) 01008 #else 01009 runups = runups1 01010 #endif 01011 ! WD: /new code 01012 ! end runup gauge computations 01013 do i=1,par%npointvar 01014 mnem = par%pointvars(i) 01015 j = chartoindex(mnem) 01016 ! lookup the proper array 01017 call indextos(sl,j,t) 01018 ! get the proper output points .... 01019 select case(t%type) 01020 case('r') 01021 do ii = 1, par%npoints + par%nrugauge 01022 ! wwvv above line should probably be 01023 ! do ii = 1, par%npoints + par%nrugauge 01024 ! have to check this with trunk 01025 if(xomaster) then 01026 pointoutputs(i,ii)%name = mnem 01027 pointoutputs(i,ii)%rank = t%rank 01028 endif 01029 if (xomaster) then 01030 xpii = xpoints(ii) 01031 ypii = ypoints(ii) 01032 endif 01033 #ifdef USEMPI 01034 call xmpi_bcast(xpii,xmpi_omaster,xmpi_ocomm) 01035 call xmpi_bcast(ypii,xmpi_omaster,xmpi_ocomm) 01036 call space_who_has(sl,xpii,ypii,pii) 01037 if (xmpi_orank .eq. pii) then 01038 call space_global_to_local(sl, xpii, ypii, ilocal,jlocal) 01039 endif 01040 #else 01041 ilocal = xpii 01042 jlocal = ypii 01043 pii = xmpi_orank 01044 #endif 01045 select case(t%rank) 01046 case(2) 01047 #ifdef USEMPI 01048 if(xmpi_orank .eq. pii) then 01049 allocate(r2(1,1)) 01050 call gridrotate(par, sl, t, r2, ilocal, jlocal) 01051 call xmpi_send(pii,xmpi_omaster,r2) 01052 elseif(xomaster) then 01053 allocate(pointoutputs(i,ii)%r2(1,1)) 01054 call xmpi_send(pii,xmpi_omaster,pointoutputs(i,ii)%r2) 01055 endif 01056 if (xmpi_orank .eq. pii) deallocate(r2) 01057 #else 01058 allocate(pointoutputs(i,ii)%r2(1,1)) 01059 call gridrotate(par, sl, t, pointoutputs(i,ii)%r2, ilocal, jlocal) 01060 #endif 01061 case(3) 01062 #ifdef USEMPI 01063 if (xmpi_orank .eq. pii) then 01064 allocate(r3(1,1,size(t%r3,3))) 01065 call gridrotate(par, sl, t, r3, ilocal, jlocal) 01066 call xmpi_send(pii,xmpi_omaster,r3) 01067 elseif(xomaster) then 01068 allocate(pointoutputs(i,ii)%r3(1,1,size(t%r3,3))) 01069 call xmpi_send(pii,xmpi_omaster,pointoutputs(i,ii)%r3) 01070 endif 01071 if (xmpi_orank .eq. pii) deallocate(r3) 01072 #else 01073 allocate(pointoutputs(i,ii)%r3(1,1,size(t%r3,3))) 01074 call gridrotate(par, sl, t, pointoutputs(i,ii)%r3, ilocal, jlocal) 01075 #endif 01076 case(4) 01077 #ifdef USEMPI 01078 if (xmpi_orank .eq. pii) then 01079 allocate(r4(1,1,size(t%r4,3),size(t%r4,4))) 01080 call gridrotate(t, r4, ilocal, jlocal) 01081 call xmpi_send(pii,xmpi_omaster,r4) 01082 elseif(xomaster) then 01083 allocate(pointoutputs(i,ii)%r4(1,1,size(t%r4,3),size(t%r4,4))) 01084 call xmpi_send(pii,xmpi_omaster,pointoutputs(i,ii)%r4) 01085 endif 01086 if (xmpi_orank .eq. pii) deallocate(r4) 01087 #else 01088 allocate(pointoutputs(i,ii)%r4(1,1,size(t%r3,3),size(t%r4,4))) 01089 call gridrotate(t, pointoutputs(i,ii)%r4, ilocal, jlocal) 01090 #endif 01091 case default 01092 write(0,*) 'Can''t handle rank: ', t%rank, ' of mnemonic', mnem 01093 end select 01094 end do 01095 case default 01096 write(0,*) 'Can''t handle type: ', t%type, ' of mnemonic', mnem 01097 end select 01098 enddo ! i=1,par%npointvar 01099 ! Here follows go the code to determine the runup values and put them 01100 ! in runups(:,1:par%nrugauge) 01101 ! 01102 ! This code assumes that hh and zs are available on xomaster 01103 ! we will change this asap 01104 endif ! dooutput_point 01105 01106 ! writing is done by xomaster, the others processes go back to work 01107 01108 if( .not. xomaster) return 01109 01110 ! Open the output file 01111 01112 #ifdef USENETCDF 01113 if(donetcdf) then 01114 NF90(nf90_open(ncid=ncid, path=par%ncfilename, mode=nf90_write)) 01115 endif 01116 #endif 01117 01118 ! 01119 ! some variables can only be output if others are available, see the code 01120 ! for gridrotate. 01121 ! 01122 if (dooutput_global) then 01123 itg = itg+1 01124 ! Store the time (in morphological time) 01125 #ifdef USENETCDF 01126 if(donetcdf) then 01127 NF90(nf90_put_var(ncid, globaltimevarid, CONVREAL(par%t*max(par%morfac,1.d0)), (/tpar%itg/))) 01128 endif 01129 #endif 01130 ! write global output variables 01131 do i=1,par%nglobalvar 01132 mnem = par%globalvars(i) 01133 j = chartoindex(mnem) 01134 ! lookup the proper array (should have been collected already) 01135 call indextos(s,j,t) 01136 01137 select case(t%type) 01138 case('i') 01139 select case(t%rank) 01140 case(0) 01141 ! no need to allocate here 01142 call gridrotate(t,i0) 01143 #ifdef USENETCDF 01144 if(donetcdf) then 01145 NF90(nf90_put_var(ncid, globalvarids(i), i0, start=(/1,tpar%itg/) )) 01146 endif 01147 #endif 01148 if(dofortran) then 01149 inquire(iolength=reclen) i0 01150 call checkfile(i,unit,reclen,jtg) 01151 write(unit,rec=jtg) i0 01152 call flush(unit) 01153 endif 01154 case(2) 01155 allocate(i2(size(t%i2,1),size(t%i2,2))) 01156 call gridrotate(t, i2) 01157 #ifdef USENETCDF 01158 if(donetcdf) then 01159 NF90(nf90_put_var(ncid, globalvarids(i), i2, start=(/1,1,tpar%itg/) )) 01160 endif 01161 #endif 01162 if(dofortran) then 01163 inquire(iolength=reclen) i2 01164 call checkfile(i,unit,reclen,jtg) 01165 write(unit,rec=jtg) i2 01166 call flush(unit) 01167 endif 01168 deallocate(i2) 01169 case(3) 01170 allocate(i3(size(t%i3,1),size(t%i3,2),size(t%i3,3))) 01171 call gridrotate(t, i3) 01172 #ifdef USENETCDF 01173 if(donetcdf) then 01174 NF90(nf90_put_var(ncid, globalvarids(i), i3, start=(/1,1,1,tpar%itg/) )) 01175 endif 01176 #endif 01177 if(dofortran) then 01178 inquire(iolength=reclen) i3 01179 call checkfile(i,unit,reclen,jtg) 01180 write(unit,rec=jtg) i3 01181 call flush(unit) 01182 endif 01183 deallocate(i3) 01184 case default 01185 write(0,*) 'Can''t handle rank: ', t%rank, ' of mnemonic', mnem 01186 end select 01187 case('r') 01188 select case(t%rank) 01189 case(0) 01190 ! no need to allocate here 01191 r0conv = CONVREAL(t%r0) 01192 #ifdef USENETCDF 01193 if(donetcdf) then 01194 NF90(nf90_put_var(ncid, globalvarids(i), r0conv, start=(/1,tpar%itg/) )) 01195 endif 01196 #endif 01197 if(dofortran) then 01198 inquire(iolength=reclen) r0conv 01199 call checkfile(i,unit,reclen,jtg) 01200 write(unit,rec=jtg) r0conv 01201 call flush(unit) 01202 endif 01203 case(1) 01204 allocate(r1(size(t%r1,1))) 01205 allocate(r1conv(size(t%r1,1))) 01206 ! no need to rotate here 01207 r1conv = CONVREAL(t%r1) 01208 #ifdef USENETCDF 01209 if(donetcdf) then 01210 NF90(nf90_put_var(ncid, globalvarids(i), r1conv, start=(/1,tpar%itg/) )) 01211 endif 01212 #endif 01213 if(dofortran) then 01214 inquire(iolength=reclen) r1conv 01215 call checkfile(i,unit,reclen,jtg) 01216 write(unit,rec=jtg) r1conv 01217 call flush(unit) 01218 endif 01219 deallocate(r1) 01220 deallocate(r1conv) 01221 case(2) 01222 allocate(r2 (size(t%r2,1),size(t%r2,2))) 01223 allocate(r2conv(size(t%r2,1),size(t%r2,2))) 01224 call gridrotate(par, s, t, r2) 01225 if(par%remdryoutput==1) then 01226 call postprocessvar_r2(s%wetz,t%name,r2,dFill) 01227 endif 01228 r2conv = CONVREAL(r2) 01229 #ifdef USENETCDF 01230 if(donetcdf) then 01231 NF90(nf90_put_var(ncid, globalvarids(i), r2conv, start=(/1,1,tpar%itg/) )) 01232 endif 01233 #endif 01234 if(dofortran) then 01235 inquire(iolength=reclen) r2conv 01236 call checkfile(i,unit,reclen,jtg) 01237 write(unit,rec=jtg) r2conv 01238 call flush(unit) 01239 endif 01240 deallocate(r2) 01241 deallocate(r2conv) 01242 case(3) 01243 allocate(r3 (size(t%r3,1),size(t%r3,2),size(t%r3,3))) 01244 allocate(r3conv(size(t%r3,1),size(t%r3,2),size(t%r3,3))) 01245 call gridrotate(par, s, t, r3) 01246 r3conv = CONVREAL(r3) 01247 #ifdef USENETCDF 01248 if(donetcdf) then 01249 NF90(nf90_put_var(ncid, globalvarids(i), r3conv, start=(/1,1,1, tpar%itg/) )) 01250 endif 01251 #endif 01252 if(dofortran) then 01253 inquire(iolength=reclen) r3conv 01254 call checkfile(i,unit,reclen,jtg) 01255 write(unit,rec=jtg) r3conv 01256 call flush(unit) 01257 endif 01258 deallocate(r3) 01259 deallocate(r3conv) 01260 case(4) 01261 allocate(r4 (size(t%r4,1),size(t%r4,2),size(t%r4,3),size(t%r4,4))) 01262 allocate(r4conv(size(t%r4,1),size(t%r4,2),size(t%r4,3),size(t%r4,4))) 01263 call gridrotate(t, r4) 01264 r4conv = CONVREAL(r4) 01265 #ifdef USENETCDF 01266 if(donetcdf) then 01267 NF90(nf90_put_var(ncid, globalvarids(i), r4conv, start=(/1,1,1,1, tpar%itg/) )) 01268 endif 01269 #endif 01270 if(dofortran) then 01271 inquire(iolength=reclen) r4conv 01272 call checkfile(i,unit,reclen,jtg) 01273 write(unit,rec=jtg) r4conv 01274 call flush(unit) 01275 endif 01276 deallocate(r4) 01277 deallocate(r4conv) 01278 case default 01279 write(0,*) 'Can''t handle rank: ', t%rank, ' of mnemonic', mnem 01280 end select 01281 case default 01282 write(0,*) 'Can''t handle type: ', t%type, ' of mnemonic', mnem 01283 end select 01284 end do 01285 end if ! dooutput_global 01286 01287 01288 if(dooutput_point) then 01289 itp = itp+1 01290 01291 #ifdef USENETCDF 01292 if(donetcdf) then 01293 NF90(nf90_put_var(ncid, pointtimevarid, CONVREAL(par%t*max(par%morfac,1.d0)), (/tpar%itp/))) 01294 endif 01295 #endif 01296 01297 if(dofortran_compat) then 01298 if (par%npoints .gt. 0) then 01299 allocate(points(par%npoints,par%npointvar+1)) 01300 else 01301 allocate(points(0,0)) 01302 endif 01303 else 01304 allocate(points(0,0)) 01305 endif 01306 01307 do i=1,par%npointvar 01308 mnem = par%pointvars(i) 01309 j = chartoindex(mnem) 01310 ! lookup the proper array 01311 call indextos(s,j,t) 01312 ! get the proper output points .... 01313 ! I have no idea what is happening in varouput so I'll try it in a different way 01314 !TODO This is not very efficient because we are using the outer counters, reorder dimensions.... 01315 ! 01316 select case(t%type) 01317 case('r') 01318 do ii = 1, par%npoints + par%nrugauge 01319 select case(t%rank) 01320 case(2) 01321 ! This postprocessing creates an ugly dependency. 01322 ! it would be nice if we could call gridrotate as a function 01323 ! or if we could just have the postprocessing insert some reference processing routines to call 01324 ! or if we could split this out of the case statement (dry) 01325 ! or if we could defer this to a postprocessing routine (for example ncks) 01326 01327 #ifdef USENETCDF 01328 if (donetcdf) then 01329 NF90(nf90_put_var(ncid, pointsvarids(i), CONVREAL(pointoutputs(i,ii)%r2(1,1)), start=(/ii,tpar%itp/) )) 01330 endif 01331 #endif 01332 if(dofortran_compat) then 01333 if (ii .le. par%npoints) then 01334 !points(ii,i+1) = r2(xpoints(ii), ypoints(ii)) 01335 points(ii,i+1) = pointoutputs(i,ii)%r2(1,1) 01336 endif 01337 endif 01338 case(3) 01339 #ifdef USENETCDF 01340 if (donetcdf) then 01341 NF90(nf90_put_var(ncid, pointsvarids(i), CONVREAL(r3(xpoints(ii), ypoints(ii),:)), start=(/ii,1,tpar%itp/))) 01342 endif 01343 #endif 01344 if(dofortran_compat) then 01345 if (ii .le. par%npoints) then 01346 points(ii,i+1) = r3(xpoints(ii), ypoints(ii),1) ! wwvv todo 01347 endif 01348 endif 01349 case(4) 01350 #ifdef USENETCDF 01351 if (donetcdf) then 01352 status = nf90_put_var(ncid, pointsvarids(i), & 01353 & CONVREAL(r4(xpoints(ii), ypoints(ii),:,:)), & 01354 & start=(/ii,1,1,tpar%itp/) ) 01355 if (status /= nf90_noerr) call handle_err(status,__FILE__,__LINE__) 01356 endif 01357 #endif 01358 if(dofortran_compat) then 01359 if (ii .le. par%npoints) then 01360 points(ii,i+1) = r4(xpoints(ii), ypoints(ii),1,1) ! wwvv todo 01361 endif 01362 endif 01363 case default 01364 write(0,*) 'Can''t handle rank: ', t%rank, ' of mnemonic', mnem 01365 end select 01366 end do 01367 case default 01368 write(0,*) 'Can''t handle type: ', t%type, ' of mnemonic', mnem 01369 end select 01370 enddo ! i=1,par%npointvar 01371 01372 if(dofortran_compat) then 01373 do ii = 1,par%npoints 01374 if (par%morfacopt==1) then 01375 points(ii,1) = par%t*max(par%morfac,1.d0) 01376 else 01377 points(ii,1) = par%t 01378 endif 01379 enddo 01380 do ii = 1,par%npoints 01381 write(indextopointsunit(ii),rec=tpar%itp)CONVREAL(points(ii,:)) 01382 call flush(indextopointsunit(ii)) 01383 enddo 01384 deallocate(points) 01385 01386 ! WD: new code 01387 do ii=1,par%nrugauge 01388 write(indextopointsunit(ii+par%npoints),rec=tpar%itp)CONVREAL(runups(:,ii)) 01389 call flush(indextopointsunit(i+par%npoints)) 01390 enddo 01391 ! WD: /new code 01392 endif ! dofortran 01393 endif ! dooutput_point 01394 01395 01396 #ifdef USENETCDF 01397 if(dooutput_mean .and. donetcdf) then 01398 ! Store the time (in morphological time) 01399 NF90(nf90_put_var(ncid, meantimevarid, CONVREAL(par%t*max(par%morfac,1.d0)), (/tpar%itm-1/))) 01400 ! write global output variables 01401 do i=1,par%nmeanvar 01402 t = meansparsglobal(i)%t 01403 do j=1,nmeanvartypes 01404 01405 select case(t%type) 01406 case('r') 01407 select case(t%rank) 01408 case(2) 01409 select case(meanvartypes(j)) 01410 case('mean') 01411 if ((t%name .eq. 'H') .or. (t%name .eq. 'urms')) then 01412 status = nf90_put_var(ncid, meanvarids(i,j), & 01413 & CONVREAL(sqrt(meansparsglobal(i)%variancesquareterm2d)), & 01414 & start=(/1,1,tpar%itm-1/) ) 01415 if (status /= nf90_noerr) call handle_err(status,__FILE__,__LINE__) 01416 elseif (t%name .eq. 'thetamean') then 01417 status = nf90_put_var(ncid, meanvarids(i,j), & 01418 & CONVREAL( & 01419 & mod(2.d0*par%px & 01420 & + atan2(nint(meansparsglobal(i)%mean2d)/1d7, & 01421 & mod(meansparsglobal(i)%mean2d,1.d0)*1d1), 2.d0*par%px) / par%px * 180), & 01422 & start=(/1,1,tpar%itm-1/) ) 01423 if (status /= nf90_noerr) call handle_err(status,__FILE__,__LINE__) 01424 else 01425 NF90(nf90_put_var(ncid, meanvarids(i,j), CONVREAL(meansparsglobal(i)%mean2d),start=(/1,1,tpar%itm-1/))) 01426 end if 01427 case('var') 01428 NF90(nf90_put_var(ncid, meanvarids(i,j), CONVREAL(meansparsglobal(i)%variance2d),start=(/1,1,tpar%itm-1/))) 01429 case('min') 01430 NF90(nf90_put_var(ncid, meanvarids(i,j), CONVREAL(meansparsglobal(i)%min2d),start=(/1,1,tpar%itm-1/))) 01431 case('max') 01432 NF90(nf90_put_var(ncid, meanvarids(i,j), CONVREAL(meansparsglobal(i)%max2d),start=(/1,1,tpar%itm-1/))) 01433 case default 01434 write(0,*) 'Can''t handle cell method: ', trim(meanvartypes(j)), ' of mnemonic', trim(t%name) 01435 end select 01436 case(3) 01437 select case(meanvartypes(j)) 01438 case('mean') 01439 NF90(nf90_put_var(ncid, meanvarids(i,j), CONVREAL(meansparsglobal(i)%mean3d), start=(/1,1,1,tpar%itm-1/) )) 01440 case('var') 01441 NF90(nf90_put_var(ncid,meanvarids(i,j),CONVREAL(meansparsglobal(i)%variance3d),start=(/1,1,1,tpar%itm-1/))) 01442 case('min') 01443 NF90(nf90_put_var(ncid, meanvarids(i,j), CONVREAL(meansparsglobal(i)%min3d), start=(/1,1,1,tpar%itm-1/) )) 01444 case('max') 01445 NF90(nf90_put_var(ncid, meanvarids(i,j), CONVREAL(meansparsglobal(i)%max3d), start=(/1,1,1,tpar%itm-1/) )) 01446 case default 01447 write(0,*) 'Can''t handle cell method: ', trim(meanvartypes(j)), ' of mnemonic', trim(t%name) 01448 end select 01449 case(4) 01450 select case(meanvartypes(j)) 01451 case('mean') 01452 NF90(nf90_put_var(ncid, meanvarids(i,j),CONVREAL(meansparsglobal(i)%mean4d), start=(/1,1,1,1,tpar%itm-1/))) 01453 case('var') 01454 NF90(nf90_put_var(ncid,meanvarids(i,j),CONVREAL(meansparsglobal(i)%variance4d),start=(/1,1,1,1,tpar%itm-1/))) 01455 case('min') 01456 NF90(nf90_put_var(ncid, meanvarids(i,j), CONVREAL(meansparsglobal(i)%min4d), start=(/1,1,1,1,tpar%itm-1/))) 01457 case('max') 01458 NF90(nf90_put_var(ncid, meanvarids(i,j), CONVREAL(meansparsglobal(i)%max4d), start=(/1,1,1,1,tpar%itm-1/))) 01459 case default 01460 write(0,*) 'Can''t handle cell method: ', trim(meanvartypes(j)), ' of mnemonic', trim(t%name) 01461 end select 01462 case default 01463 write(0,*) 'Can''t handle rank: ', t%rank, ' of mnemonic', mnem 01464 end select 01465 case default 01466 write(0,*) 'Can''t handle type: ', t%type, ' of mnemonic', mnem 01467 end select 01468 end do 01469 end do 01470 endif ! dooutput_mean .and. donetcdf 01471 #endif 01472 01473 !if (par%nmeanvar>0 .and. dofortran .and. tpar%output) then 01474 ! Not at the first in tpm as this is the start of averaging. Only output after second in tpm 01475 if(dooutput_mean .and. dofortran_compat) then 01476 itm=itm+1 ! Note, this is a local counter, used to position in output file 01477 do i=1,par%nmeanvar 01478 select case (meansparsglobal(i)%rank) 01479 case (2) 01480 if (par%meanvars(i)=='H') then ! Hrms changed to H 01481 write(indextomeanunit(i),rec=itm)CONVREAL(sqrt(meansparsglobal(i)%variancesquareterm2d)) 01482 elseif (par%meanvars(i)=='urms') then ! urms 01483 write(indextomeanunit(i),rec=itm)CONVREAL(sqrt(meansparsglobal(i)%variancesquareterm2d)) 01484 elseif (par%meanvars(i)=='thetamean') then ! thetamean 01485 write(indextomeanunit(i),rec=itm) & 01486 & CONVREAL( & 01487 & mod(2.d0*par%px + atan2(nint(meansparsglobal(i)%mean2d)/1d7, & 01488 & mod(meansparsglobal(i)%mean2d,1.d0)*1d1), 2.d0*par%px) / par%px * 180 & 01489 & ) 01490 else ! non-rms variables 01491 write(indextomeanunit(i),rec=itm)CONVREAL(meansparsglobal(i)%mean2d) 01492 endif 01493 write(indextovarunit(i),rec=itm)CONVREAL(meansparsglobal(i)%variance2d) 01494 where(meansparsglobal(i)%min2d>0.99d0*huge(0.d0)) 01495 meansparsglobal(i)%min2d=-999.d0 01496 endwhere 01497 where(meansparsglobal(i)%max2d<-0.99d0*huge(0.d0)) 01498 meansparsglobal(i)%max2d=-999.d0 01499 endwhere 01500 write(indextominunit(i),rec=itm)CONVREAL(meansparsglobal(i)%min2d) 01501 write(indextomaxunit(i),rec=itm)CONVREAL(meansparsglobal(i)%max2d) 01502 case (3) 01503 write(indextomeanunit(i),rec=itm)CONVREAL(meansparsglobal(i)%mean3d) 01504 write(indextovarunit(i),rec=itm)CONVREAL(meansparsglobal(i)%variance3d) 01505 where(meansparsglobal(i)%min3d>0.99d0*huge(0.d0)) 01506 meansparsglobal(i)%min3d=-999.d0 01507 endwhere 01508 where(meansparsglobal(i)%max3d<-0.99d0*huge(0.d0)) 01509 meansparsglobal(i)%max3d=-999.d0 01510 endwhere 01511 write(indextominunit(i),rec=itm)CONVREAL(meansparsglobal(i)%min3d) 01512 write(indextomaxunit(i),rec=itm)CONVREAL(meansparsglobal(i)%max3d) 01513 case (4) 01514 write(indextomeanunit(i),rec=itm)CONVREAL(meansparsglobal(i)%mean4d) 01515 write(indextovarunit(i),rec=itm)CONVREAL(meansparsglobal(i)%variance4d) 01516 where(meansparsglobal(i)%min4d>0.99d0*huge(0.d0)) 01517 meansparsglobal(i)%min4d=-999.d0 01518 endwhere 01519 where(meansparsglobal(i)%max4d<-0.99d0*huge(0.d0)) 01520 meansparsglobal(i)%max4d=-999.d0 01521 endwhere 01522 write(indextominunit(i),rec=itm)CONVREAL(meansparsglobal(i)%min4d) 01523 write(indextomaxunit(i),rec=itm)CONVREAL(meansparsglobal(i)%max4d) 01524 end select 01525 call flush(indextomeanunit(i)) 01526 call flush(indextovarunit(i)) 01527 call flush(indextominunit(i)) 01528 call flush(indextomaxunit(i)) 01529 enddo 01530 par%tintm = tpar%tpm(min(itm+2,stpm))-tpar%tpm(itm+1) ! Next averaging period (min to stop array out of bounds) 01531 par%tintm = max(par%tintm,tiny(0.d0)) ! to prevent par%tintm=0 after last output 01532 endif ! dooutput_mean .and. dofortran_compat 01533 01534 if (dooutput_hotstart) then 01535 call hotstartfiles_write(s,par, tpar) 01536 endif 01537 01538 if (dooutput_drifter) then 01539 itd = itd+1 01540 #ifdef USENETCDF 01541 if (donetcdf) then 01542 ! output time: 01543 NF90(nf90_put_var(ncid, drifterstimevarid, CONVREAL(par%t), start=(/itd/))) 01544 endif 01545 #endif 01546 do i=1,par%ndrifter 01547 if ( par%t>=s%tdriftb(i) .and. par%t<=s%tdrifte(i) .and. & 01548 & s%idrift(i)>1 .and. s%idrift(i)<=s%nx .and. & 01549 & s%jdrift(i)>1 .and. s%jdrift(i)<=s%ny ) then 01550 01551 iz = int(s%idrift(i)) 01552 jz = int(s%jdrift(i)) 01553 01554 di = mod(s%idrift(i),1.d0) 01555 dj = mod(s%jdrift(i),1.d0) 01556 01557 dx = di*s%dsu(iz,jz)*cos(s%alfaz(iz,jz)) - & 01558 & dj*s%dnv(iz,jz)*sin(s%alfaz(iz,jz)) 01559 dy = di*s%dsu(iz,jz)*sin(s%alfaz(iz,jz)) + & 01560 & dj*s%dnv(iz,jz)*cos(s%alfaz(iz,jz)) 01561 01562 drift(1) = s%xz(iz,jz)+dx 01563 drift(2) = s%yz(iz,jz)+dy 01564 drift(3) = par%t 01565 else 01566 drift(1) = -999 01567 drift(2) = -999 01568 drift(3) = par%t 01569 endif 01570 01571 if (dofortran_compat) then 01572 write(indextodrifterunit(i),rec=itd) CONVREAL(drift) 01573 call flush(indextodrifterunit(i)) 01574 endif 01575 #ifdef USENETCDF 01576 if (donetcdf) then 01577 ! output x,y: 01578 NF90(nf90_put_var(ncid,driftersvarids(i),CONVREAL(drift(1:2)),start=(/1,itd/))) 01579 endif 01580 #endif 01581 01582 enddo 01583 endif 01584 #ifdef USENETCDF 01585 if(donetcdf) then 01586 NF90(nf90_close(ncid=ncid)) 01587 endif 01588 #endif 01589 01590 if(dofortran_compat) then 01591 outputtimes = -999.d0 01592 outputtimes(1:itg) = tpar%tpg(1:itg) 01593 outputtimes(itg+1:itg+itp) = tpar%tpp(1:itp) 01594 outputtimes(itg+itp+1:itg+itp+itm) = tpar%tpm(2:itm+1) ! mean output always shifted by 1 01595 01596 if (par%morfacopt==1) outputtimes=outputtimes*max(par%morfac,1.d0) 01597 01598 open(1998,file='dims.dat',form='unformatted',access='direct',recl=wordsize*(10+size(outputtimes))) 01599 write(1998,rec=1) CONVREAL(itg*1.d0),& 01600 & CONVREAL(s%nx*1.d0),& 01601 & CONVREAL(s%ny*1.d0),& 01602 & CONVREAL(s%ntheta*1.d0),& 01603 & CONVREAL(par%kmax*1.d0),& 01604 & CONVREAL(par%ngd*1.d0),& 01605 & CONVREAL(par%nd*1.d0), & 01606 & CONVREAL(itp*1.d0),& 01607 & CONVREAL(itc*1.d0),& 01608 & CONVREAL(itm*1.d0),& 01609 & CONVREAL(outputtimes) 01610 call flush(1998) 01611 endif 01612 ! wwvv avoid warning about unused sl: 01613 if (sl%nx .eq. -1) return 01614 end subroutine ncoutput 01615 01616 #ifdef USENETCDF 01617 character(slen) function dimensionnames(dimids) 01618 implicit none 01619 integer, dimension(:), intent(in) :: dimids ! store the dimids in a vector 01620 01621 integer :: i 01622 character(slen) :: dimensionname 01623 ! combine all the dimensionnames 01624 ! assumes all dimensions have an accompanying variable that should be used for coordinates. 01625 ! ",".join would have been nice here.... 01626 dimensionnames = '' 01627 ! Fortran array dimensions are in reverse order 01628 do i=size(dimids),2,-1 01629 NF90(nf90_inquire_dimension(ncid, dimids(i), name=dimensionname)) 01630 dimensionnames = trim(dimensionnames) // trim(dimensionname) // ',' 01631 end do 01632 NF90(nf90_inquire_dimension(ncid, dimids(1), name=dimensionname)) 01633 dimensionnames = trim(dimensionnames) // trim(dimensionname) 01634 end function dimensionnames 01635 01636 integer function dimensionid(expression) 01637 ! Function to transform the expression in spaceparams.tmpl to an id, we might want this in the 01638 ! makeincludes module 01639 use logging_module, only: writelog 01640 implicit none 01641 character(len=*),intent(in) :: expression 01642 select case(expression) 01643 case('s%nx+1') 01644 dimensionid = xdimid 01645 case('s%ny+1') 01646 dimensionid = ydimid 01647 case('s%ntheta') 01648 dimensionid = thetadimid 01649 case('s%tidelen') 01650 dimensionid = tidetimedimid 01651 case('par%tideloc') 01652 dimensionid = tidecornersdimid 01653 case('s%windlen') 01654 dimensionid = windtimedimid 01655 case('par%ngd') 01656 dimensionid = sedimentclassesdimid 01657 case('s%ntdisch') 01658 dimensionid = inoutdimid 01659 case('2') 01660 dimensionid = inoutdimid 01661 case('par%nd') 01662 dimensionid = bedlayersdimid 01663 case('par%ndrifter') 01664 dimensionid = driftersdimid 01665 case('par%nship') 01666 dimensionid = shipdimid 01667 case('par%nz') 01668 dimensionid = Q3Ddimid 01669 case('s%nsecvegmax') 01670 dimensionid = vegdimid 01671 case default 01672 call writelog('els','','Unknown dimension expression:' // expression) 01673 stop 1 01674 end select 01675 end function dimensionid 01676 #endif 01677 ! USENETCDF 01678 01679 subroutine points_output_init(s,par) 01680 ! this initialize things for point output 01681 ! has to be called before fourtoutput_init and ncoutput_init 01682 use spaceparams 01683 use params, only: parameters 01684 use postprocessmod, only: snappointstogrid 01685 use xmpi_module 01686 01687 type(spacepars), intent(in) :: s 01688 type(parameters),intent(in) :: par 01689 01690 allocate(rugrowindex(par%nrugauge)) 01691 01692 allocate(xpoints(par%npoints+par%nrugauge)) 01693 allocate(ypoints(par%npoints+par%nrugauge)) 01694 01695 ! Convert world coordinates of points to nearest (lsm) grid point 01696 if(xomaster) then 01697 call snappointstogrid(par, s, xpoints, ypoints) 01698 endif 01699 01700 if(xomaster) then 01701 if (par%nrugauge>0) then 01702 rugrowindex = ypoints(par%npoints+1:) 01703 !do i=1,par%nrugauge 01704 ! rugrowindex(i)=ypoints(par%npoints+i) 01705 !enddo 01706 endif 01707 endif 01708 #ifdef USEMPI 01709 call xmpi_bcast(rugrowindex,xmpi_omaster,xmpi_ocomm) 01710 call xmpi_bcast(xpoints, xmpi_omaster,xmpi_ocomm) 01711 call xmpi_bcast(ypoints, xmpi_omaster,xmpi_ocomm) 01712 #endif 01713 end subroutine points_output_init 01714 01715 subroutine fortoutput_init(s,par,tpar) 01716 use params, only: parameters 01717 use spaceparams 01718 use readkey_module 01719 use timestep_module, only: timepars 01720 use logging_module 01721 use postprocessmod 01722 use filefunctions, only: create_new_fid 01723 use mnemmodule, only: arraytype, chartoindex 01724 use xmpi_module 01725 01726 implicit none 01727 01728 type(spacepars),intent(in) :: s 01729 type(parameters),intent(in) :: par 01730 type(timepars),intent(in) :: tpar 01731 01732 integer :: i,fid 01733 integer :: reclen,reclenm,reclenp 01734 character(100) :: fname,fnamemean,fnamevar,fnamemin,fnamemax 01735 type(arraytype) :: t 01736 01737 if (.not. xomaster) return 01738 01739 reclenm = -123 01740 01741 ! Initialize places in output files 01742 itg = 0 01743 itm = 0 01744 itp = 0 01745 itc = 0 01746 itd = 0 01747 stpm = size(tpar%tpm) 01748 01749 ! Record size for global and mean output 01750 inquire(iolength=wordsize) CONVREAL(1.0d0) 01751 reclen=wordsize*(s%nx+1)*(s%ny+1) 01752 01753 open(100,file='xy.dat',form='unformatted',access='direct',recl=reclen,status='REPLACE') 01754 write(100,rec=1)CONVREAL(s%xz) 01755 write(100,rec=2)CONVREAL(s%yz) 01756 write(100,rec=3)CONVREAL(s%x) 01757 write(100,rec=4)CONVREAL(s%y) 01758 close(100) 01759 01760 ! GLOBAL VARS 01761 01762 noutnumbers = par%nglobalvar 01763 ! store all indices for the global variables 01764 do i= 1,noutnumbers 01765 outnumbers(i) = chartoindex(par%globalvars(i)) 01766 enddo 01767 01768 do i=1,par%npoints+par%nrugauge 01769 fname = trim(par%stationid(i)) // '.dat' 01770 if (par%pointtypes(i)==0) then 01771 reclenp=wordsize*(par%npointvar+1) 01772 else 01773 reclenp=wordsize*(1+par%nrugdepth*3) 01774 endif 01775 open(indextopointsunit(i),file=fname,& 01776 & form='unformatted',access='direct',recl=reclenp,status='REPLACE') 01777 enddo 01778 if (par%npoints>0) then 01779 ! write index file of point output variables 01780 fid = create_new_fid() 01781 open(fid,file='pointvars.idx',status='replace',action='write') 01782 do i=1,par%npointvar 01783 write(fid,*)trim(par%pointvars(i)) 01784 enddo 01785 close(fid) 01786 endif 01787 01788 ! TIME-AVERAGE, VARIANCE and MIN-MAX ARRAYS 01789 01790 if (par%nmeanvar>0) then 01791 !! First time file opening for time-average output 01792 do i=1,par%nmeanvar 01793 call makeaveragenames(chartoindex(par%meanvars(i)),fnamemean,fnamevar,fnamemin,fnamemax) 01794 call indextos(s,chartoindex(par%meanvars(i)),t) 01795 reclenm = wordsize 01796 select case(t%rank) 01797 case (2) 01798 reclenm = wordsize*size(t%r2) 01799 case (3) 01800 reclenm = wordsize*size(t%r3) 01801 case (4) 01802 reclenm = wordsize*size(t%r4) 01803 end select 01804 open(indextomeanunit(i),file=fnamemean, form='unformatted',access='direct',recl=reclenm,status='REPLACE') 01805 open(indextovarunit(i), file=fnamevar, form='unformatted',access='direct',recl=reclenm,status='REPLACE') 01806 open(indextominunit(i), file=fnamemin, form='unformatted',access='direct',recl=reclenm,status='REPLACE') 01807 open(indextomaxunit(i), file=fnamemax, form='unformatted',access='direct',recl=reclenm,status='REPLACE') 01808 enddo 01809 endif ! par%nmeanvar > 0 01810 01811 ! 01812 ! drifter output files 01813 ! 01814 if (par%ndrifter>0) then 01815 reclen=wordsize*3 01816 do i=1,par%ndrifter 01817 write(fname,'("drifter",i0.3,".dat")') i 01818 open(indextodrifterunit(i),file=fname,form='unformatted',access='direct',recl=reclen,status='REPLACE') 01819 enddo 01820 endif ! par%ndrifter >0 01821 01822 01823 end subroutine fortoutput_init 01824 01825 subroutine checkfile(index,unit,reclen,jtg) 01826 implicit none 01827 integer, intent(in) :: index,reclen 01828 integer, intent(out) :: unit,jtg 01829 logical :: lopen 01830 character(len=1000) :: filename 01831 01832 unit = indextoglobalunit(index) 01833 inquire(unit=unit, opened=lopen) 01834 if ( .not. lopen ) then 01835 filename = trim(mnemonics(outnumbers(index)))//'.dat' 01836 open(unit, file=filename,form='unformatted',& 01837 & access='direct',recl=reclen) 01838 endif 01839 inquire(unit=unit,nextrec=jtg) 01840 end subroutine checkfile 01841 01842 integer function outunit(ind,s) 01843 ! 01844 ! given the type of output file 's' and the index 'ind' 01845 ! returns the appopriate file unit number 01846 ! 01847 implicit none 01848 integer, intent(in) :: ind 01849 character(*), intent(in) :: s 01850 01851 integer, parameter :: offset = 10000 01852 select case(s) 01853 case('points') 01854 outunit = offset + ind 01855 case('global') 01856 outunit = offset + 10*numvars + ind 01857 case('mean') 01858 outunit = offset + 20*numvars + ind 01859 case('min') 01860 outunit = offset + 30*numvars + ind 01861 case('max') 01862 outunit = offset + 40*numvars + ind 01863 case('var') 01864 outunit = offset + 50*numvars + ind 01865 case('drifter') 01866 outunit = offset + 60*numvars + ind 01867 case('hotstart') 01868 outunit = offset + 70*numvars + ind 01869 case default 01870 print *,'internal error in outunit: no such type: ',trim(s) 01871 outunit = -1 01872 call halt_program 01873 end select 01874 end function outunit 01875 01876 integer function indextopointsunit(index) 01877 implicit none 01878 integer, intent(in) :: index 01879 indextopointsunit = outunit(index,'points') 01880 end function indextopointsunit 01881 01882 integer function indextoglobalunit(index) 01883 implicit none 01884 integer, intent(in) :: index 01885 indextoglobalunit = outunit(index,'global') 01886 end function indextoglobalunit 01887 01888 integer function indextomeanunit(index) 01889 implicit none 01890 integer, intent(in) :: index 01891 indextomeanunit = outunit(index,'mean') 01892 end function indextomeanunit 01893 01894 integer function indextominunit(index) 01895 implicit none 01896 integer, intent(in) :: index 01897 indextominunit = outunit(index,'min') 01898 end function indextominunit 01899 01900 integer function indextomaxunit(index) 01901 implicit none 01902 integer, intent(in) :: index 01903 indextomaxunit = outunit(index,'max') 01904 end function indextomaxunit 01905 01906 integer function indextovarunit(index) 01907 implicit none 01908 integer, intent(in) :: index 01909 indextovarunit = outunit(index,'var') 01910 end function indextovarunit 01911 01912 integer function indextodrifterunit(index) 01913 implicit none 01914 integer, intent(in) :: index 01915 indextodrifterunit = outunit(index,'drifter') 01916 end function indextodrifterunit 01917 01918 subroutine makeaveragenames(counter,fnamemean,fnamevar,fnamemin& 01919 & ,fnamemax) 01920 use mnemmodule 01921 01922 implicit none 01923 01924 character(*) :: fnamemean,fnamevar,fnamemin,fnamemax 01925 integer :: counter 01926 01927 fnamemean = trim(mnemonics(counter))//'_mean.dat' 01928 fnamevar = trim(mnemonics(counter))//'_var.dat' 01929 fnamemin = trim(mnemonics(counter))//'_min.dat' 01930 fnamemax = trim(mnemonics(counter))//'_max.dat' 01931 01932 end subroutine makeaveragenames 01933 01934 01935 subroutine hotstartfiles_collect(s,sl,par, tpar) 01936 use logging_module 01937 #ifdef USEMPI 01938 use xmpi_module 01939 #endif 01940 use params 01941 use spaceparams 01942 use timestep_module 01943 use mnemmodule 01944 use postprocessmod 01945 use paramsconst 01946 01947 type(spacepars), intent(inout) :: s ! s-> spaceparams, what is isl? 01948 ! s describes the global system, sl the local system 01949 ! in this MPI process 01950 type(spacepars), intent(inout) :: sl ! s-> spaceparams, what is isl? 01951 type(parameters), intent(inout) :: par 01952 type(timepars), intent(in) :: tpar 01953 01954 ! Water level, bed, velocities. Used in all models 01955 #ifdef USEMPI 01956 call space_collect_mnem(s,sl,par,mnem_zs) 01957 call space_collect_mnem(s,sl,par,mnem_zb) 01958 call space_collect_mnem(s,sl,par,mnem_uu) 01959 call space_collect_mnem(s,sl,par,mnem_vv) 01960 01961 01962 ! groundwater parameters 01963 if(par%gwflow==1) then 01964 call space_collect_mnem(s,sl,par,mnem_gwlevel) 01965 call space_collect_mnem(s,sl,par,mnem_gwhead) 01966 call space_collect_mnem(s,sl,par,mnem_gwcurv) 01967 endif 01968 01969 ! wave energy 01970 if (par%swave==1) then 01971 call space_collect_mnem(s,sl,par,mnem_ee) 01972 call space_collect_mnem(s,sl,par,mnem_rr) 01973 endif 01974 01975 ! single-dir 01976 if (par%single_dir==1) then 01977 call space_collect_mnem(s,sl,par,mnem_ee_s) 01978 endif 01979 01980 ! nonh parameters 01981 if (par%wavemodel == WAVEMODEL_NONH) then 01982 call space_collect_mnem(s,sl,par,mnem_breaking) 01983 call space_collect_mnem(s,sl,par,mnem_wb) 01984 call space_collect_mnem(s,sl,par,mnem_ws) 01985 endif 01986 01987 ! sediment transport 01988 if (par%sedtrans==1 .and. (par%sus==1 .or. par%bulk==1)) then 01989 call space_collect_mnem(s,sl,par,mnem_ccg) 01990 endif 01991 01992 ! turbulence 01993 if (par%turb==1) then 01994 call space_collect_mnem(s,sl,par,mnem_kturb) 01995 endif 01996 01997 ! structures 01998 if (par%struct==1) then 01999 call space_collect_mnem(s,sl,par,mnem_structdepth) 02000 endif 02001 02002 ! nhplus 02003 if (par%nonhq3d==1) then 02004 call space_collect_mnem(s,sl,par,mnem_dU) 02005 call space_collect_mnem(s,sl,par,mnem_dV) 02006 endif 02007 02008 ! boundary flow 02009 if ((par%flow==1).or.(par%wavemodel==WAVEMODEL_NONH)) then 02010 call space_collect_mnem(s,sl,par,mnem_umean) 02011 call space_collect_mnem(s,sl,par,mnem_vmean) 02012 endif 02013 #endif 02014 02015 end subroutine hotstartfiles_collect 02016 02017 subroutine hotstartfiles_write(s,par, tpar) 02018 use logging_module 02019 #ifdef USEMPI 02020 use xmpi_module 02021 #endif 02022 use params 02023 use spaceparams 02024 use timestep_module 02025 use mnemmodule 02026 use postprocessmod 02027 use paramsconst 02028 02029 type(spacepars), intent(inout) :: s ! s-> spaceparams, what is isl? 02030 ! s describes the global system, sl the local system 02031 ! in this MPI process 02032 type(parameters), intent(inout) :: par 02033 type(timepars), intent(in) :: tpar 02034 02035 if (xomaster) then 02036 ! Water level, bed, velocities. Used in all models 02037 if (xomaster) then 02038 call lowlevel_write_hotstart_2d(s%zs,s%nx,s%ny,'zs',tpar%ith) 02039 call lowlevel_write_hotstart_2d(s%zb,s%nx,s%ny,'zb',tpar%ith) 02040 call lowlevel_write_hotstart_2d(s%uu,s%nx,s%ny,'uu',tpar%ith) 02041 call lowlevel_write_hotstart_2d(s%vv,s%nx,s%ny,'vv',tpar%ith) 02042 endif 02043 02044 ! groundwater parameters 02045 if(par%gwflow==1) then 02046 call lowlevel_write_hotstart_2d(s%gwlevel,s%nx,s%ny,'gwlevel',tpar%ith) 02047 call lowlevel_write_hotstart_2d(s%gwhead,s%nx,s%ny,'gwhead',tpar%ith) 02048 call lowlevel_write_hotstart_2d(s%gwcurv,s%nx,s%ny,'gwcurv',tpar%ith) 02049 endif 02050 02051 ! wave energy 02052 if (par%swave==1) then 02053 call lowlevel_write_hotstart_3d(s%ee,s%nx,s%ny,s%ntheta,'ee',tpar%ith) 02054 call lowlevel_write_hotstart_3d(s%rr,s%nx,s%ny,s%ntheta,'rr',tpar%ith) 02055 endif 02056 02057 ! single-dir 02058 if (par%single_dir==1) then 02059 call lowlevel_write_hotstart_3d(s%ee_s,s%nx,s%ny,s%ntheta_s,'ee_s',tpar%ith) 02060 endif 02061 02062 ! nonh parameters 02063 if (par%wavemodel == WAVEMODEL_NONH) then 02064 call lowlevel_write_hotstart_2d_int(s%breaking,s%nx,s%ny,'breaking',tpar%ith) 02065 call lowlevel_write_hotstart_2d(s%wb,s%nx,s%ny,'wb',tpar%ith) 02066 call lowlevel_write_hotstart_2d(s%ws,s%nx,s%ny,'ws',tpar%ith) 02067 endif 02068 02069 ! sediment transport 02070 if (par%sedtrans==1 .and. (par%sus==1 .or. par%bulk==1)) then 02071 call lowlevel_write_hotstart_3d(s%ccg,s%nx,s%ny,par%ngd,'ccg',tpar%ith) 02072 endif 02073 02074 ! turbulence 02075 if (par%turb==1) then 02076 call lowlevel_write_hotstart_2d(s%kturb,s%nx,s%ny,'kturb',tpar%ith) 02077 endif 02078 02079 ! structures 02080 if (par%struct==1) then 02081 call lowlevel_write_hotstart_2d(s%structdepth,s%nx,s%ny,'structdepth',tpar%ith) 02082 endif 02083 02084 ! nhplus 02085 if (par%nonhq3d==1) then 02086 call lowlevel_write_hotstart_2d(s%dU,s%nx,s%ny,'dU',tpar%ith) 02087 call lowlevel_write_hotstart_2d(s%dV,s%nx,s%ny,'dV',tpar%ith) 02088 endif 02089 02090 ! boundary flow 02091 if ((par%flow==1).or.(par%wavemodel==WAVEMODEL_NONH)) then 02092 call lowlevel_write_hotstart_2d(s%umean,s%nx,s%ny,'umean',tpar%ith) 02093 call lowlevel_write_hotstart_2d(s%vmean,s%nx,s%ny,'vmean',tpar%ith) 02094 endif 02095 endif ! xomaster 02096 end subroutine hotstartfiles_write 02097 02098 subroutine lowlevel_write_hotstart_2d(var,nx,ny,varname,ith) 02099 02100 integer,intent(in) :: nx,ny,ith 02101 real*8,dimension(nx+1,ny+1),intent(in) :: var 02102 character(*), intent(in) :: varname 02103 02104 integer :: unit,i,j,ierr 02105 character(128) :: fname,rowfmt 02106 02107 unit = outunit(1,'hotstart') 02108 write(fname,'(a,a,I6.6,a)') 'hotstart_',trim(varname),ith,'.dat' 02109 write(rowfmt,'(a,I12,a)') '(',ny+1,'(1X,D))' 02110 open(unit,file=fname,form='FORMATTED', iostat=ierr) 02111 do j=1,ny+1 02112 write(unit,FMT=rowfmt) (var(i,j), i=1,nx+1) 02113 enddo 02114 close(unit) 02115 02116 end subroutine lowlevel_write_hotstart_2d 02117 02118 subroutine lowlevel_write_hotstart_2d_int(var,nx,ny,varname,ith) 02119 02120 integer,intent(in) :: nx,ny,ith 02121 integer,dimension(nx+1,ny+1),intent(in) :: var 02122 character(*), intent(in) :: varname 02123 02124 integer :: unit,i,j,ierr 02125 character(128) :: fname,rowfmt 02126 02127 unit = outunit(2,'hotstart') 02128 write(fname,'(a,a,I6.6,a)') 'hotstart_',trim(varname),ith,'.dat' 02129 write(rowfmt,'(a,I12,a)') '(',ny+1,'(1X,I4))' 02130 open(unit,file=fname,form='FORMATTED', iostat=ierr) 02131 do j=1,ny+1 02132 write(unit,FMT=rowfmt) (var(i,j), i=1,nx+1) 02133 enddo 02134 close(unit) 02135 02136 end subroutine lowlevel_write_hotstart_2d_int 02137 02138 subroutine lowlevel_write_hotstart_3d(var,nx,ny,nd,varname,ith) 02139 02140 integer,intent(in) :: nx,ny,nd,ith 02141 real*8,dimension(nx+1,ny+1,nd),intent(in) :: var 02142 character(*), intent(in) :: varname 02143 02144 integer :: unit,i,j,k,ierr 02145 character(128) :: fname,rowfmt 02146 02147 unit = outunit(3,'hotstart') 02148 write(fname,'(a,a,I6.6,a)') 'hotstart_',trim(varname),ith,'.dat' 02149 write(rowfmt,'(a,I12,a)') '(',ny+1,'(1X,D))' 02150 open(unit,file=fname,form='FORMATTED', iostat=ierr) 02151 do k=1,nd 02152 do j=1,ny+1 02153 write(unit,FMT=rowfmt) (var(i,j,k), i=1,nx+1) 02154 enddo 02155 enddo 02156 close(unit) 02157 02158 end subroutine lowlevel_write_hotstart_3d 02159 02160 02161 end module ncoutput_module
1.7.6.1
