!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! - Ariane - (May - 2007) !! !! bruno.blanke@univ-brest.fr and nicolas.grima@univ-brest.fr !! !! This software is a computer program whose purpose is !! the computation of 3D streamlines in a given velocity field !! (as the output of an Ocean General Circulation Model) and !! subsequent water masses analyses. !! !! This software is governed by the CeCILL license under French law and !! abiding by the rules of distribution of free software. You can use, !! modify and/ or redistribute the software under the terms of the CeCILL !! license as circulated by CEA, CNRS and INRIA at the following URL !! "http://www.cecill.info". !! !! As a counterpart to the access to the source code and rights to copy, !! modify and redistribute granted by the license, users are provided only !! with a limited warranty and the software's author, the holder of the !! economic rights, and the successive licensors have only limited !! liability. !! !! In this respect, the user's attention is drawn to the risks associated !! with loading, using, modifying and/or developing or reproducing the !! software by the user in light of its specific status of free software, !! that may mean that it is complicated to manipulate, and that also !! therefore means that it is reserved for developers and experienced !! professionals having in-depth computer knowledge. Users are therefore !! encouraged to load and test the software's suitability as regards their !! requirements in conditions enabling the security of their systems and/or !! data to be ensured and, more generally, to use and operate it in the !! same conditions as regards security. !! !! The fact that you are presently reading this means that you have had !! knowledge of the CeCILL license and that you accept its terms. !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Module mod_init_particules Use mod_precision USE mod_memory Use mod_namelist Use mod_lun Use mod_input_data Use mod_posin Use mod_seq Implicit None Integer( kind = iprec) :: & ! maxsect = 0 , & ! quantitative mode = max of sections maxsegm = 0 , & ! quantitative mode = max of segments ntraj = 0 , & ! number of trajectories nsect = 0 , & ! number of sections nfnt = 0 , & ! icrit1 = 0 Real(kind = rprec) :: & ! trtot = 0._rprec Contains !========================================================================= !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !======================================================================== ! four options to initialize the positions: ! - qualitative: on ASCII file "initial_positions.txt" ! - qualitative: on BINARY file "initial.bin" ! - quantitative: auto positionning (from file "segments") ! - quantitative: from file "initial.bin" !========================================================================= !========================================================================= Subroutine sub_init_particules_positions() Integer(kind = iprec) :: nb_dims Integer(kind = iprec), Dimension(1) :: dims !================! !- Declarations -! !================! Integer(kind=iprec) :: & ! n , & ! ind , & ! nb_subset , & ! alloc_size , & ! ncid , & ! varid Integer(kind=iprec) :: is_loop Real(kind=rprec) :: & ! trmin = 1._rprec ! A particle is set if the transport is ! higher thant this value. ! Reduce this value if your transport is ! very small. Integer(kind=iprec), Dimension(:), Allocatable :: ind_subset Write(lun_standard,*)'' Write(lun_standard,*)'===============================' Write(lun_standard,*)'= INITIAL PARTICULE POSITIONS =' Write(lun_standard,*)'===============================' Write(lun_standard,*)' (read or computed)' !=====================! !- Quantitative MODE -! !=====================! If (Trim(mode) == 'quantitative') Then !-------------------------------------------------! !- Compute in sections.txt : maxsect and maxsegm -! !-------------------------------------------------! Call sub_init_part_maxsect_maxsegm() !---------------------! !- Auto positioning - ! !---------------------! If (key_sequential) Then id_comments = .True. !- Write input data reading information -! Call posini_seq(trmin,ntraj,nsect,trtot,nfnt) Else Call posini(trmin, ntraj, nsect, trtot, nfnt) Endif !---------------------------------------- ! initialization of transparent sections !---------------------------------------- !!NG IF (nfnt > 0) THEN !!NG DO lu = 91, 90 + nfnt !!NG WRITE(fname,'(a12,a15)')'transparent.',fntname(lu-90) !!NG OPEN(lu,form='FORMATTED',file=fname) !!NG END DO !!NG ENDIF Endif !=====================================================! !- indices read on ASCII file (QUALITATIVE analysis) -! !=====================================================! If ((Trim(mode) == 'qualitative').And.(Trim(bin) == 'nobin')) Then Write(lun_output,*) Write(lun_output,*)'======================================================' Write(lun_output,*)'ASCII positions are read on file initial_positions.txt' Write(lun_output,*)'======================================================' Open(unit=lun_dummy, file='initial_positions.txt', action="read") Do n = 1, nmax + 1 Read(lun_dummy,*,End=7054,ERR=7054)tfi(n),tfj(n),tfk(n),tfl(n),ttr(n) End Do Write(lun_error,*) Write(lun_error,*)'# of positions in initial.positions.txt exceeds nmax' Write(lun_error,*)'--- STOP ARIANE ---' Stop 7054 ntraj = n - 1 Write(lun_output,*) Write(lun_output,*)'# of particles that could be read: ',ntraj If ( ntraj == 0) Then Write(lun_error,*) Write(lun_error,*)')-: The number of particule is 0 :-( ' Write(lun_error,*)'Please verify that you have correctly set the file' Write(lun_error,*)'initial.positions.txt (param0 is obsoleted).' Write(lun_error,*)'--- STOP ARIANE ---' Endif If (key_sequential) Then IF (key_alltracers) THEN IF ((key_roms).OR.(key_symphonie))THEN alloc_size = 7 ELSE IF (TRIM(w_surf_option) == 'E-P-R') THEN alloc_size = 7 ELSE alloc_size = 6 ENDIF ENDIF ELSE IF ((key_roms).OR.(key_symphonie))THEN alloc_size = 3 ELSE IF (TRIM(w_surf_option) == 'E-P-R') THEN alloc_size = 4 ELSE alloc_size = 3 ENDIF ENDIF ENDIF !!$ If (key_alltracers) Then !!$ IF ((key_roms).OR.(key_symphonie))THEN !!$ alloc_size = 7 !!$ ELSE !!$ alloc_size = 6 !!$ ENDIF !!$ Else !!$ alloc_size = 3 !!$ Endif Call sub_seq_alloc(alloc_size) Endif Endif !===========================================! !- initial positions in file "initial.bin" -! !===========================================! If (Trim(bin) /= 'nobin') Then !----------------------------- ! indices read on BINARY file !----------------------------- !!NG WRITE(lun_output,*) !!NG WRITE(lun_output,*)'===============================================' !!NG WRITE(lun_output,*)'BINARY positions are read from file initial.bin' !!NG WRITE(lun_output,*)'===============================================' !!NG !!NG OPEN(UNIT=lun_initial,file='initial.bin', form='UNFORMATTED', ACTION='READ') !!NG !!NG DO n = 1, nmax + 1 !!NG !!NG READ(UNIT=lun_initial,END=7053,ERR=7053)tfi(n),tfj(n),tfk(n),tfl(n),ttr(n) !!NG !!NG trtot = trtot + ttr(n) !!NG !!NG END DO !!NG !!NG WRITE(lun_error,*) !!NG WRITE(lun_error,*)'# of available positions in initial.bin exceeds nmax' !!NG WRITE(lun_error,*)'--- STOP ARIANE ---' !!NG STOP !!NG !!NG7053 ntraj = n - 1 Write(lun_output,*) Write(lun_output,*)'==================================================' Write(lun_output,*)'= Positions are read from ariane_initial.nc file =' Write(lun_output,*)'==================================================' !------------------------! !- Open the netcdf file -! !------------------------! Call sub_open_netcdf_file('.', 'ariane_initial.nc', ncid) Write(lun_standard, *)' --- ariane_initial.nc is opened ---' !-------! !- TFI -! (rien à voir avec la chaine TF1 ;-) !-------! Call sub_read_netcdf_varid_ndims( & ncid = ncid , & nc_var = Trim(init_final)//'_x', & varid = varid , & ndims = nb_dims ) Call sub_read_netcdf_var_dims( & ncid = ncid , & varid = varid , & ndims = nb_dims , & dims = dims(:) ) ntraj = dims(1) Write(lun_output,*) Write(lun_output,*)'# of particles that could be read: ',ntraj Call sub_read_netcdf_var1d( ncid, varid, tfi(1:ntraj)) Write(lun_standard, *)' - X positions are read from ', & &Trim(init_final)//'_x', ' -' !-------! !- TFJ -! !-------! Call sub_read_netcdf_varid_ndims(ncid, Trim(init_final)//'_y', varid) Call sub_read_netcdf_var1d( ncid, varid, tfj(1:ntraj)) Write(lun_standard, *)' - Y positions are read from ', & &Trim(init_final)//'_y', ' -' !-------! !- TFK -! !-------! Call sub_read_netcdf_varid_ndims(ncid, Trim(init_final)//'_z', varid) Call sub_read_netcdf_var1d( ncid, varid, tfk(1:ntraj)) Write(lun_standard, *)' - Z positions are read from ', & &Trim(init_final)//'_z', ' -' !-------! !- TFL -! !-------! Call sub_read_netcdf_varid_ndims(ncid, Trim(init_final)//'_t', varid) Call sub_read_netcdf_var1d( ncid, varid, tfl(1:ntraj)) Write(lun_standard, *)' - Time positions read from ', & &Trim(init_final)//'_t', ' -' !-------! !- AGE -! !-------! If (key_read_age) Then Call sub_read_netcdf_varid_ndims(ncid, Trim(init_final)//'_age', varid) Call sub_read_netcdf_var1d( ncid, varid, tage(1:ntraj)) trtot = Sum(tage(1:ntraj)) / Real(ntraj, kind = rprec) Write(lun_standard, *)' - Age are read from ', & &Trim(init_final)//'_age', ' - (mean =', trtot, ')' Endif !-------! !- TTR -! (transport) !-------! Call sub_read_netcdf_varid_ndims(ncid, Trim(init_final)//'_transp', varid) Call sub_read_netcdf_var1d( ncid, varid, ttr(1:ntraj)) trtot = Sum(ttr(1:ntraj)) Write(lun_standard, *)' - Transports are read from ', & &Trim(init_final)//'_transp', ' - (total = ', trtot, ')' If (key_alltracers.And.(Trim(mode) == 'quantitative')) Then !---------------! !- TEMPERATURE -! !---------------! Call sub_read_netcdf_varid_ndims(ncid, Trim(init_final)//'_temp', varid) Call sub_read_netcdf_var1d( ncid, varid, init_temp(1:ntraj)) Write(lun_standard, *)' - Temperature values are read from ', & &Trim(init_final)//'_temp', ' -' !------------! !- SALYNITY -! !------------! Call sub_read_netcdf_varid_ndims(ncid, Trim(init_final)//'_salt', varid) Call sub_read_netcdf_var1d( ncid, varid, init_salt(1:ntraj)) Write(lun_standard, *)' - Salinity values are read from ', & &Trim(init_final)//'_salt', ' -' !-----------! !- DENSITY -! !-----------! Call sub_read_netcdf_varid_ndims(ncid, Trim(init_final)//'_dens', varid) Call sub_read_netcdf_var1d( ncid, varid, init_dens(1:ntraj)) Write(lun_standard, *)' - Density values are read from ', & &Trim(init_final)//'_dens', ' -' Endif !-------------------------! !- Close the netcdf file -! !-------------------------! Call sub_close_netcdf_file(ncid) Write(lun_standard, *)' --- ariane_initial.nc is closed ---' Endif !=========================! !- Read a subset of data -! !=========================! If (Trim(bin) == 'subbin') Then Allocate(ind_subset(nmax)) call sub_memory(size(ind_subset),'i','ind_subset','sub_init_particules_positions') ind_subset(:)=0 Open(lun_subset, file='subset.txt') Write(lun_output,*)'A subset of indices is reading on file subset.txt.' Do nb_subset = 1, ntraj+1 Read(lun_subset,*,End=1010,ERR=1010)ind_subset(nb_subset) End Do Write(lun_error,*) Write(lun_error,*)'# of indices in subset exceeds ntraj' Write(lun_error,*)'--- STOP ARIANE ---' Stop 1010 Continue Close (lun_subset) nb_subset = nb_subset - 1 If (Maxval(ind_subset(:)) > ntraj) Then Write(lun_error,*)'Error: mod_init_particules.f90' Write(lun_error,*)' An indice in the subset.txt file is' Write(lun_error,*)' higher than the number of particules !!!' Write(lun_error,*)'--- STOP ARIANE ---' Stop Endif trtot = 0._rprec ! Because the indices of the subset particules increase ! we can not use an intermediate array. ! This means that ind-subset(ind) is always >= at ind. ! A little bit dangerous... Do ind = 1, nb_subset tfi(ind) = tfi(ind_subset(ind)) tfj(ind) = tfj(ind_subset(ind)) tfk(ind) = tfk(ind_subset(ind)) tfl(ind) = tfl(ind_subset(ind)) ttr(ind) = ttr(ind_subset(ind)) tage(ind) = tage(ind_subset(ind)) trtot = trtot + ttr(ind) if (key_alltracers.And.(Trim(mode) == 'quantitative')) THEN init_temp(ind) = init_temp(ind_subset(ind)) init_salt(ind) = init_salt(ind_subset(ind)) init_dens(ind) = init_dens(ind_subset(ind)) endif End Do ntraj = nb_subset Write(lun_output,*) Write(lun_output,*)'# of particles that were kept as a subset: ', ntraj call sub_memory(-size(ind_subset),'i','ind_subset','sub_init_particules_positions') Deallocate(ind_subset) Endif !========================================! !- Write initial positions in ASCI file -! !========================================! If ((Trim(mode) == 'quantitative').And.key_ascii_outputs ) Then Write(lun_standard,*)' - Writing initial positions...' Do is_loop = 1, ntraj Write(lun_init_pos,7050) & tfi(is_loop),tfj(is_loop),tfk(is_loop),tfl(is_loop), & ttr(is_loop) Enddo Close(lun_init_pos) Endif 7050 Format(5(1x,f0.3)) !---------------------------------------------! !- Compute lmin and lmax in qualitative mode -! !---------------------------------------------! If (Trim(mode) == 'qualitative') Then If(Trim(forback) == 'forward') Then Write(lun_standard,*) & ' - Computing lmin in qualitative mode...' lmin = Nint(Minval(tfl, mask = tfl > 0 )) Write(lun_standard,*) ' - lmin =', lmin IF (lmin == 0) THEN Write(lun_error,*) & ' lmin = 0, there is a problem in the initial_positions.txt file' STOP ENDIF Elseif(Trim(forback) == 'backward') Then Write(lun_standard,*) & ' - Computing lmax in qualitative mode...' lmax = Nint(Maxval(tfl, mask = tfl < lmt+1)) Write(lun_standard,*) ' - lmax =', lmax Endif Endif !======================================! !- Allocate arrays in sequential mode -! !======================================! If ((Trim(bin) /= 'nobin')) Then If (key_sequential) Then IF (key_alltracers) THEN IF ((key_roms).OR.(key_symphonie))THEN alloc_size = 7 ELSE IF (TRIM(w_surf_option) == 'E-P-R') THEN alloc_size = 7 ELSE alloc_size = 6 ENDIF ENDIF ELSE IF ((key_roms).OR.(key_symphonie))THEN alloc_size = 3 ELSE IF (TRIM(w_surf_option) == 'E-P-R') THEN alloc_size = 4 ELSE alloc_size = 3 ENDIF ENDIF ENDIF !!$ If (key_alltracers) Then !!$ IF ((key_roms).OR.(key_symphonie))THEN !!$ alloc_size = 7 !!$ ELSE !!$ alloc_size = 6 !!$ ENDIF !!$ Else !!$ alloc_size = 3 !!$ Endif Call sub_seq_alloc(alloc_size) Endif Endif !======================================================================== ! I/O - I/O - I/O - I/O - I/O - I/O - I/O - I/O - I/O - I/O - I/O - ! [QUANT] !======================================================================== If (Trim(mode) == 'quantitative') Then Write(lun_output,*)' ' Write(lun_output,*)'QUANTITATIVE EXPERIMENT' Write(lun_output,*)'=======================' Write(lun_output,*)' minimum gridcell transport: ', trmin Write(lun_output,*)' maximum individual transport: ', max_transport Write(lun_output,*)' exact number of segments: ', maxsegm Write(lun_output,*)' exact number of particles: ', ntraj Write(lun_output,*)' exact number of sections: ', nsect Write(lun_output,*)' total documented transport: ', trtot Write(lun_output,*)'number of transparent sections: ', nfnt Write(lun_output,*)' ' Endif End Subroutine sub_init_particules_positions !========================================================================= !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !========================================================================= Subroutine sub_init_part_maxsect_maxsegm() Integer(kind=iprec) :: idum !-----------------------------------------------------------------! !- Open the file sections.txt and compute the number of segments -! !- and the number of sections. -! !- (remember that a section could be composed by more than -! !- one segment) -! !-----------------------------------------------------------------! Open(unit=lun_dummy, file='sections.txt', action="read") Do Read(lun_dummy,*, End=7777, ERR=7777) idum maxsegm = maxsegm + 1 If (idum > maxsect) Then maxsect = idum Endif Enddo 7777 Continue Write(lun_output,*)'Number of segments in sections.txt file = ', maxsegm Close(lun_dummy) End Subroutine sub_init_part_maxsect_maxsegm End Module mod_init_particules