algae_transp.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\algae_transp.f
00002 !
00075       MODULE ALGAE_TRANSP
00076         USE BIEF_DEF, ONLY : BIEF_OBJ
00077       IMPLICIT NONE
00078 !***********************************************************************
00079 ! TELEMAC2D   V6P3                                   14/06/2013
00080 ! EDF R&D                                           antoine.joly@edf.fr
00081 !***********************************************************************
00082 !
00083 !
00084 !
00085 ! SUBROUTINES MADE AVAILABLE
00086       PUBLIC :: ALLOC_ALGAE,INTERP_ALGAE,DISP_ALGAE
00087 ! MODEL VARIABLES
00088       INTEGER :: ALGAE_START
00089       DATA ALGAE_START / 1 / ! THIS VALUE NEEDS TO BE UPDATED IN FLOT
00090 ! MEAN FLUID VARIABLES AT THE POSITION OF EACH ALGAE PARTICLES
00091       TYPE(BIEF_OBJ):: U_X_AV_0
00092       TYPE(BIEF_OBJ):: U_Y_AV_0
00093       TYPE(BIEF_OBJ):: U_Z_AV_0
00094       TYPE(BIEF_OBJ):: U_X_AV
00095       TYPE(BIEF_OBJ):: U_Y_AV
00096       TYPE(BIEF_OBJ):: U_Z_AV
00097       TYPE(BIEF_OBJ):: K_AV_0
00098       TYPE(BIEF_OBJ):: EPS_AV_0
00099       TYPE(BIEF_OBJ):: K_AV
00100       TYPE(BIEF_OBJ):: EPS_AV
00101       TYPE(BIEF_OBJ):: H_FLU
00102 ! VARIABLES OF THE BODIES
00103       TYPE(BIEF_OBJ):: U_X_0
00104       TYPE(BIEF_OBJ):: U_Y_0
00105       TYPE(BIEF_OBJ):: U_Z_0
00106       TYPE(BIEF_OBJ):: U_X
00107       TYPE(BIEF_OBJ):: U_Y
00108       TYPE(BIEF_OBJ):: U_Z
00109       TYPE(BIEF_OBJ):: V_X_0
00110       TYPE(BIEF_OBJ):: V_Y_0
00111       TYPE(BIEF_OBJ):: V_Z_0
00112       TYPE(BIEF_OBJ):: V_X
00113       TYPE(BIEF_OBJ):: V_Y
00114       TYPE(BIEF_OBJ):: V_Z
00115       TYPE(BIEF_OBJ):: DX_A
00116       TYPE(BIEF_OBJ):: DY_A
00117       TYPE(BIEF_OBJ):: DZ_A
00118       TYPE(BIEF_OBJ):: I_A_GL
00119 ! VARIABLES USED TO CALCULATE THE BASSET HISTORY FORCE
00120       INTEGER:: NB
00121       INTEGER:: IB
00122       INTEGER:: NP
00123       INTEGER:: IP
00124       DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE:: T_TIL_P
00125       DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE:: A_P
00126       INTEGER:: NWIN
00127       INTEGER:: IWIN
00128       DOUBLE PRECISION:: TWIN
00129       DOUBLE PRECISION,DIMENSION(:,:,:),ALLOCATABLE:: PSI
00130       DOUBLE PRECISION:: PHI_PLUS
00131       DOUBLE PRECISION:: PHI_MOINS
00132       DOUBLE PRECISION:: CB
00133       DOUBLE PRECISION,DIMENSION(:,:,:),ALLOCATABLE:: FI_P
00134       DOUBLE PRECISION:: F_TAIL
00135       DOUBLE PRECISION:: CI_BAS
00136 !
00137       CONTAINS
00138 !
00139 !                    **********************
00140                      SUBROUTINE ALLOC_ALGAE
00141 !                    **********************
00142 !
00143      &(NP_TOT,MESH,DT)
00144 !
00145 !***********************************************************************
00146 ! TELEMAC 2D VERSION 6.3    MAI 2013                       ANTOINE JOLY
00147 ! EDF R&D                                           antoine.joly@edf.fr
00148 !***********************************************************************
00149 !
00150 !
00151 !-----------------------------------------------------------------------
00152 !                             ARGUMENTS
00153 ! .________________.____.______________________________________________.
00154 ! |  NAME          |MODE|                  ROLE                        |
00155 ! |________________|____|______________________________________________|
00156 ! | NP_TOT         | -->| TOTAL NUMBER OF ALGAE PARTICLES              |
00157 ! | MESH           | -->| MESH STRUCTURE WITH ALL THE INFORMATIONS     |
00158 ! |                |    | OF THE MESH TREATED                          |
00159 ! | DT             | -->| NUMERICAL TIME STEP OF THE SIMULATIONS       |
00160 ! |________________|____|______________________________________________|
00161 ! MODE : -->(NON MODIFIED DATA), <--(RESULT), <-->(MODIFIED DATA)
00162 !
00163 !-----------------------------------------------------------------------
00164 !
00165 ! CALLED BY : DERIVE
00166 !
00167 ! SUBROUTINE CALLED : INIT_BASSET
00168 !
00169 !***********************************************************************
00170 !
00171       USE BIEF
00172       USE STREAMLINE
00173 !
00174       IMPLICIT NONE
00175       INTEGER LNG,LU
00176       COMMON/INFO/LNG,LU
00177 !
00178 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00179 !
00180       INTEGER         , INTENT(IN) :: NP_TOT
00181       TYPE(BIEF_MESH) , INTENT(IN) :: MESH
00182       DOUBLE PRECISION, INTENT(IN) :: DT
00183 !
00184 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00185 !
00186 !=======================================================================
00187 ! ALLOCATE THE VARIABLES
00188 !=======================================================================
00189 ! MEAN FLUID VARIABLES
00190       CALL BIEF_ALLVEC(1,U_X_AV_0,'UX_AV0',NP_TOT,1,0,MESH)
00191       CALL BIEF_ALLVEC(1,U_Y_AV_0,'UY_AV0',NP_TOT,1,0,MESH)
00192       CALL BIEF_ALLVEC(1,U_Z_AV_0,'UZ_AV0',NP_TOT,1,0,MESH)
00193       CALL BIEF_ALLVEC(1,U_X_AV,'U_X_AV',NP_TOT,1,0,MESH)
00194       CALL BIEF_ALLVEC(1,U_Y_AV,'U_Y_AV',NP_TOT,1,0,MESH)
00195       CALL BIEF_ALLVEC(1,U_Z_AV,'U_Z_AV',NP_TOT,1,0,MESH)
00196       CALL BIEF_ALLVEC(1,K_AV_0,'K_AV_0',NP_TOT,1,0,MESH)
00197       CALL BIEF_ALLVEC(1,EPS_AV_0,'E_AV_0',NP_TOT,1,0,MESH)
00198       CALL BIEF_ALLVEC(1,K_AV,'K_AVER',NP_TOT,1,0,MESH)
00199       CALL BIEF_ALLVEC(1,EPS_AV,'E_AVER',NP_TOT,1,0,MESH)
00200       CALL BIEF_ALLVEC(1,H_FLU,'H_FL_A',NP_TOT,1,0,MESH)
00201 ! VARIABLES OF THE BODIES
00202       CALL BIEF_ALLVEC(1,U_X_0,'UF_X_0',NP_TOT,1,0,MESH)
00203       CALL BIEF_ALLVEC(1,U_Y_0,'UF_Y_0',NP_TOT,1,0,MESH)
00204       CALL BIEF_ALLVEC(1,U_Z_0,'UF_Z_0',NP_TOT,1,0,MESH)
00205       CALL BIEF_ALLVEC(1,U_X,'UFLU_X',NP_TOT,1,0,MESH)
00206       CALL BIEF_ALLVEC(1,U_Y,'UFLU_Y',NP_TOT,1,0,MESH)
00207       CALL BIEF_ALLVEC(1,U_Z,'UFLU_Z',NP_TOT,1,0,MESH)
00208       CALL BIEF_ALLVEC(1,V_X_0,'VC_X_0',NP_TOT,1,0,MESH)
00209       CALL BIEF_ALLVEC(1,V_Y_0,'VC_Y_0',NP_TOT,1,0,MESH)
00210       CALL BIEF_ALLVEC(1,V_Z_0,'VC_Z_0',NP_TOT,1,0,MESH)
00211       CALL BIEF_ALLVEC(1,V_X,'VCOR_X',NP_TOT,1,0,MESH)
00212       CALL BIEF_ALLVEC(1,V_Y,'VCOR_Y',NP_TOT,1,0,MESH)
00213       CALL BIEF_ALLVEC(1,V_Z,'VCOR_Y',NP_TOT,1,0,MESH)
00214       CALL BIEF_ALLVEC(1,DX_A,'DX_ALG',NP_TOT,1,0,MESH)
00215       CALL BIEF_ALLVEC(1,DY_A,'DY_ALG',NP_TOT,1,0,MESH)
00216       CALL BIEF_ALLVEC(1,DZ_A,'DZ_ALG',NP_TOT,1,0,MESH)
00217       CALL BIEF_ALLVEC(2,I_A_GL,'I_A_GL',NP_TOT,1,0,MESH)
00218 ! VARIABLES USED TO CALCULATE THE BASSET HISTORY FORCE
00219       CALL INIT_BASSET(NP_TOT,MESH%DIM,DT)
00220 !
00221 !=======================================================================
00222 ! INITIALISE THE VARIABLES
00223 !=======================================================================
00224 ! MEAN FLUID VARIABLES
00225       CALL OS('X=C     ',X=U_X_AV,C=0.D0)
00226       CALL OS('X=C     ',X=U_Y_AV,C=0.D0)
00227       CALL OS('X=C     ',X=U_Z_AV,C=0.D0)
00228       CALL OS('X=C     ',X=K_AV,C=0.D0)
00229       CALL OS('X=C     ',X=EPS_AV,C=1.D0)
00230       CALL OS('X=C     ',X=H_FLU,C=1.D0)
00231 ! VARIABLES OF THE BODIES
00232       CALL OS('X=C     ',X=U_X,C=0.D0)
00233       CALL OS('X=C     ',X=U_Y,C=0.D0)
00234       CALL OS('X=C     ',X=U_Z,C=0.D0)
00235       CALL OS('X=C     ',X=V_X,C=0.D0)
00236       CALL OS('X=C     ',X=V_Y,C=0.D0)
00237       CALL OS('X=C     ',X=V_Z,C=0.D0)
00238       CALL OS('X=C     ',X=DX_A,C=0.D0)
00239       CALL OS('X=C     ',X=DY_A,C=0.D0)
00240       CALL OS('X=C     ',X=DZ_A,C=0.D0)
00241 !=======================================================================
00242 ! ALLOCATE VARIABLES USED IN PARALLEL TRANSPORT
00243 !=======================================================================
00244       IF(NCSIZE.GT.1) CALL ORGANISE_ALGS(NP_TOT,NWIN*MESH%DIM)
00245 !
00246       RETURN
00247       END SUBROUTINE ALLOC_ALGAE
00248 !
00249 !                    **********************
00250                      SUBROUTINE INIT_BASSET
00251 !                    **********************
00252 !
00253      &(N_A,NDIM,DT)
00254 !
00255 !***********************************************************************
00256 ! TELEMAC 2D VERSION 6.3    MAI 2013                       ANTOINE JOLY
00257 ! EDF R&D                                           antoine.joly@edf.fr
00258 !***********************************************************************
00259 !
00260 !
00261 !-----------------------------------------------------------------------
00262 !                             ARGUMENTS
00263 ! .________________.____.______________________________________________.
00264 ! |  NAME          |MODE|                  ROLE                        |
00265 ! |________________|____|______________________________________________|
00266 ! | NA             | -->| TOTAL NUMBER OF ALGAE PARTICLES              |
00267 ! | NDIM           | -->| NUMBER OF DIMENSIONS TREATED IN THE          |
00268 ! |                |    | SIMULATION                                   |
00269 ! | DT             | -->| NUMERICAL TIME STEP OF THE SIMULATIONS       |
00270 ! |________________|____|______________________________________________|
00271 ! MODE : -->(NON MODIFIED DATA), <--(RESULT), <-->(MODIFIED DATA)
00272 !
00273 !-----------------------------------------------------------------------
00274 !
00275 ! CALLED BY : ALLOC_ALGAE
00276 !
00277 ! SUBROUTINE CALLED : NONE
00278 !
00279 !***********************************************************************
00280 !
00281 !
00282       IMPLICIT NONE
00283 !
00284 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00285 !
00286       INTEGER         ,INTENT(IN)    :: N_A
00287       INTEGER         ,INTENT(IN)    :: NDIM
00288       DOUBLE PRECISION,INTENT(IN)    :: DT
00289       INTEGER                        :: I_A
00290       INTEGER                        :: I_DIM
00291 !
00292 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00293 !
00294 !
00295       NP=10
00296 !
00297       IF(ALLOCATED(T_TIL_P))DEALLOCATE(T_TIL_P)
00298       ALLOCATE(T_TIL_P(NP))
00299       T_TIL_P(1)=0.1D0
00300       T_TIL_P(2)=0.3D0
00301       T_TIL_P(3)=1.D0
00302       T_TIL_P(4)=3.D0
00303       T_TIL_P(5)=10.D0
00304       T_TIL_P(6)=40.D0
00305       T_TIL_P(7)=190.D0
00306       T_TIL_P(8)=1000.D0
00307       T_TIL_P(9)=6500.D0
00308       T_TIL_P(10)=50000.D0
00309 !
00310       IF(ALLOCATED(A_P))DEALLOCATE(A_P)
00311       ALLOCATE(A_P(NP))
00312       A_P(1)=0.23477481312586D0
00313       A_P(2)=0.28549576238194D0
00314       A_P(3)=0.28479416718255D0
00315       A_P(4)=0.26149775537574D0
00316       A_P(5)=0.32056200511938D0
00317       A_P(6)=0.35354490689146D0
00318       A_P(7)=0.39635904496921D0
00319       A_P(8)=0.42253908596514D0
00320       A_P(9)=0.48317384225265D0
00321       A_P(10)=0.63661146557001D0
00322 !
00323       IF(ALLOCATED(FI_P))DEALLOCATE(FI_P)
00324       ALLOCATE(FI_P(N_A,2,NP))
00325       DO I_A=1,N_A
00326         DO I_DIM=1,2
00327           DO IP=1,NP
00328             FI_P(I_A,I_DIM,IP)=0.D0
00329           END DO
00330         END DO
00331       END DO
00332 !
00333       NWIN=100
00334       TWIN=REAL(NWIN)*DT
00335 !
00336       IF(ALLOCATED(PSI))DEALLOCATE(PSI)
00337       ALLOCATE(PSI(N_A,NDIM,NWIN+1))
00338       DO I_A=1,N_A
00339         DO I_DIM=1,NDIM
00340           DO IWIN=1,NWIN+1
00341             PSI(I_A,I_DIM,IWIN)=0.D0
00342           END DO
00343         END DO
00344       END DO
00345 !
00346       RETURN
00347       END SUBROUTINE INIT_BASSET
00348 !
00349 !                    ************************
00350                      SUBROUTINE INTERP_ALGAE
00351 !                    ************************
00352 !
00353      &(NP,NP_TOT,SHP_P,SHZ_P,ELT_P,U_X_AV,U_Y_AV,U_Z_AV,K_AV,EPS_AV,
00354      & H_FLU,NPOIN,IELM,NDP,NPLAN,NELMAX,IKLE,W1,
00355      & IELMU,NPOINU,UCONV,VCONV,WCONV,
00356      & AK,EP,H)
00357 !
00358 !***********************************************************************
00359 ! TELEMAC 2D VERSION 6.3    MAI 2013                       ANTOINE JOLY
00360 ! EDF R&D                                           antoine.joly@edf.fr
00361 !***********************************************************************
00362 !
00363 !
00364 !
00365 !-----------------------------------------------------------------------
00366 !                             ARGUMENTS
00367 ! .________________.____.______________________________________________.
00368 ! |  NAME          |MODE|                  ROLE                        |
00369 ! |________________|____|______________________________________________|
00370 ! | NP_TOT         | -->| TOTAL NUMBER OF ALGAE PARTICLES              |
00371 ! | NP             | -->| NUMBER OF ALGAE PARTICLES IN THE CURRENT     |
00372 ! |                |    | PROCESSOR                                    |
00373 ! | SHP_P          | -->| 2D BARYCENTRIC COORDINATES FOR EACH ALGAE    |
00374 ! |                |    | PARTICLE                                     |
00375 ! | ELT_P          | -->| NUMBER OF THE ELEMENT CONTAINING THE ALGAE   |
00376 ! |                |    | PARTICLE                                     |
00377 ! | U_X_AV,U_Y_AV  |<-- | MEAN FLUID VELOCITIES AT THE POSITION OF     |
00378 ! |                |    | EACH ALGAE PARTICLE                          |
00379 ! | K_AV,EPS_AV    |<-- | TURBULENT PROPERTIES OF THE FLOW AT THE      |
00380 ! |                |    | POSITION OF EACH ALGAE PARTICLE              |
00381 ! | H_FLU          |<-- | WATER DEPTH AT THE POSITION OF EACH ALGAE    |
00382 ! |                |    | PARTICLE                                     |
00383 ! | NPOIN          | -->| NUMBER OF POINTS IN THE MESH                 |
00384 ! | IELM           | -->| ELEMENT TYPE OF THE MESH (ASSUMED THE SAME   |
00385 ! |                |    | FOR K, EPS AND H)                            |
00386 ! | NDP            | -->| NUMBER OF NODES PER ELEMENTS                 |
00387 ! | NPLAN          | -->| NUMBER OF PLANES IN THE 3D MESH OF PRISMS    |
00388 ! | NELMAX         | -->| MAXIMUM NUMBER OF ELEMENTS IN 2D             |
00389 ! | IKLE           | -->| CONNECTIVITY TABLE                           |
00390 ! | W1             | -->| TEMPORARY WORK ARRAY USED AS A BUFFER FOR    |
00391 ! |                |    | FREQUENCY                                    |
00392 ! | IELMU          | -->| ELEMENT TYPE FOR U, V AND W                  |
00393 ! | NPOINU         | -->| NUMBER OF POINTS IN THE MESH FOR VARIABLES   |
00394 ! |                |    | U, V AND W                                   |
00395 ! | UCONV,VCONV,   | -->| THE X, Y AND Z COMPONENTS OF THE FLUID       |
00396 ! |   WCONV        |    | VELOCITIES AT EACH NODE OF THE MESH          |
00397 ! | AK,EPS         | -->| TURBULENT PROPERTIES OF THE FLOW AT EACH     |
00398 ! |                |    | NODE OF THE MESH                             |
00399 ! | H              | -->| WATER DEPTH AT EACH NODE OF THE MESH         |
00400 ! |________________|____|______________________________________________|
00401 ! MODE : -->(NON MODIFIED DATA), <--(RESULT), <-->(MODIFIED DATA)
00402 !
00403 !-----------------------------------------------------------------------
00404 !
00405 ! CALLED BY : DERIVE
00406 !
00407 ! SUBROUTINE CALLED : INIT_BASSET
00408 !
00409 !***********************************************************************
00410 !
00411       USE BIEF
00412       USE STREAMLINE, ONLY : BIEF_INTERP
00413 !
00414       IMPLICIT NONE
00415       INTEGER LNG,LU
00416       COMMON/INFO/LNG,LU
00417 
00418 !
00419 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00420 !
00421 ! FLOW VARIABLES THAT ARE USED IN THE INTERPOLATION
00422       INTEGER         , INTENT(IN)     :: NDP,NPLAN,NELMAX
00423       INTEGER         , INTENT(IN)     :: NPOIN,IELM,IELMU,NPOINU
00424       DOUBLE PRECISION, INTENT(IN)     :: UCONV(NPOIN),VCONV(NPOIN)
00425       DOUBLE PRECISION, INTENT(IN)     :: WCONV(NPOIN)
00426       DOUBLE PRECISION, INTENT(IN)     :: AK(NPOIN),EP(NPOIN)
00427       DOUBLE PRECISION, INTENT(IN)     :: H(NPOIN)
00428 ! ALGAE VARIABLES THAT NEED TO BE INTERPOLATED
00429       INTEGER         , INTENT(IN)     :: NP,NP_TOT
00430       DOUBLE PRECISION, INTENT(IN)     :: SHP_P(NDP,NP_TOT)
00431       INTEGER         , INTENT(IN)     :: IKLE(NELMAX,NDP)
00432       DOUBLE PRECISION, INTENT(IN)     :: SHZ_P(NP_TOT)
00433       INTEGER         , INTENT(IN)     :: ELT_P(NP_TOT)
00434       DOUBLE PRECISION, INTENT(OUT)    :: U_X_AV(NP_TOT)
00435       DOUBLE PRECISION, INTENT(OUT)    :: U_Y_AV(NP_TOT)
00436       DOUBLE PRECISION, INTENT(OUT)    :: U_Z_AV(NP_TOT)
00437       DOUBLE PRECISION, INTENT(OUT)    :: K_AV(NP_TOT)
00438       DOUBLE PRECISION, INTENT(OUT)    :: EPS_AV(NP_TOT)
00439       DOUBLE PRECISION, INTENT(OUT)    :: H_FLU(NP_TOT)
00440 ! BUFFERS (USED ON FREQUENCY)
00441       DOUBLE PRECISION, INTENT(IN)    :: W1(NP_TOT)
00442       INTEGER FRE(1)
00443 !
00444 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00445 !
00446 !=======================================================================
00447 ! INITIALISE THE VARIABLES
00448 !=======================================================================
00449 ! U,V,W
00450       CALL BIEF_INTERP(UCONV,U_X_AV,SHP_P,NDP,SHZ_P,ELT_P,W1,FRE,
00451      &           ELT_P,NP,NPOINU,NPLAN,IELMU,IKLE,NELMAX,
00452      &           .FALSE.,.FALSE.)
00453       CALL BIEF_INTERP(VCONV,U_Y_AV,SHP_P,NDP,SHZ_P,ELT_P,W1,FRE,
00454      &           ELT_P,NP,NPOINU,NPLAN,IELMU,IKLE,NELMAX,
00455      &           .FALSE.,.FALSE.)
00456       CALL BIEF_INTERP(WCONV,U_Z_AV,SHP_P,NDP,SHZ_P,ELT_P,W1,FRE,
00457      &           ELT_P,NP,NPOINU,NPLAN,IELMU,IKLE,NELMAX,
00458      &           .FALSE.,.FALSE.)
00459 ! K, EPS
00460       CALL BIEF_INTERP(AK,K_AV,SHP_P,NDP,SHZ_P,ELT_P,W1,FRE,
00461      &           ELT_P,NP,NPOIN,1,IELM,IKLE,NELMAX,
00462      &           .FALSE.,.FALSE.)
00463       CALL BIEF_INTERP(EP,EPS_AV,SHP_P,NDP,SHZ_P,ELT_P,W1,FRE,
00464      &           ELT_P,NP,NPOIN,1,IELM,IKLE,NELMAX,
00465      &           .FALSE.,.FALSE.)
00466 ! H
00467       CALL BIEF_INTERP(H,H_FLU,SHP_P,NDP,SHZ_P,ELT_P,W1,FRE,
00468      &           ELT_P,NP,NPOIN,1,IELM,IKLE,NELMAX,
00469      &           .FALSE.,.FALSE.)
00470 !
00471       RETURN
00472       END SUBROUTINE INTERP_ALGAE
00473 !
00474 !                    ********************
00475                      SUBROUTINE DISP_ALGAE
00476 !                    ********************
00477 !
00478      & (NA_TOT,NA,NDIM,DT,ALGAE_START,U_X_AV_0,U_Y_AV_0,U_Z_AV_0,K_AV_0,
00479      &  EPS_AV_0,H_FLU,U_X_AV,U_Y_AV,U_Z_AV,U_X_0,U_Y_0,U_Z_0,V_X_0,
00480      &  V_Y_0,V_Z_0,DX_A,DY_A,DZ_A,ELEM_ALG,U_X,U_Y,U_Z,V_X,V_Y,V_Z,
00481      &  X_A,Y_A,Z_A,LT,DALGAE,RALGAE,EALGAE,ALGTYP)
00482 !
00483 !***********************************************************************
00484 ! TELEMAC 2D VERSION 6.3    MAI 2013                       ANTOINE JOLY
00485 ! EDF R&D                                           antoine.joly@edf.fr
00486 !***********************************************************************
00487 !
00488 !
00489 !
00490 !
00491 !-----------------------------------------------------------------------
00492 !                             ARGUMENTS
00493 ! .________________.____.______________________________________________.
00494 ! |  NAME          |MODE|                  ROLE                        |
00495 ! |________________|____|______________________________________________|
00496 ! | NA_TOT         | -->| TOTAL NUMBER OF ALGAE PARTICLES              |
00497 ! | NA             | -->| NUMBER OF ALGAE PARTICLES IN THE CURRENT     |
00498 ! |                |    | PROCESSOR                                    |
00499 ! | NDIM           | -->| NUMBER OF DIMENSIONS TREATED                 |
00500 ! | DT             | -->| TIME STEP OF THE SIMULATION                  |
00501 ! | ALGAE_START    | -->| TIME STEP AT WHICH ALGAE PARTICLES ARE       |
00502 ! |                |    | RELEASED                                     |
00503 ! | U_X_AV_0,      | -->| MEAN FLUID VELOCITIES AT THE POSITION OF     |
00504 ! |   U_Y_AV_0,    |    | EACH ALGAE PARTICLE ONE TIME STEP BEFORE     |
00505 ! |     U_Z_AV_0   |    | THE CALCULATIONS                             |
00506 ! | U_X_AV,U_Y_AV, | -->| MEAN FLUID VELOCITIES AT THE POSITION OF     |
00507 ! |   U_Y_AV       |    | EACH ALGAE PARTICLE                          |
00508 ! | K_AV_0,        | -->| TURBULENT PROPERTIES OF THE FLOW AT THE      |
00509 ! |   EPS_AV_0     |    | POSITION OF EACH ALGAE PARTICLE ONE TIME     |
00510 ! |                |    | STEP BEFORE THE CALCULATIONS                 |
00511 ! | H_FLU          | -->| WATER DEPTH AT THE POSITION OF EACH ALGAE    |
00512 ! |                |    | PARTICLE                                     |
00513 ! | U_X_0,U_Y_0,   | -->| TURBULENT FLUID VELOCITIES AT THE POSITION   |
00514 ! |   U_Z_0        |    | OF EACH ALGAE PARTICLE ONE TIME STEP BEFORE  |
00515 ! |                |    | THE CALCULATIONS                             |
00516 ! | U_X,U_Y,U_Z    |<-- | TURBULENT FLUID VELOCITIES AT THE POSITION   |
00517 ! |                |    | OF EACH ALGAE PARTICLE                       |
00518 ! | V_X_0,V_Y_0,   | -->| BODY VELOCITIES AT THE POSITION OF EACH      |
00519 ! |   V_Z_0        |    | ALGAE PARTICLE ONE TIME STEP BEFORE          |
00520 ! |                |    | THE CALCULATIONS                             |
00521 ! | V_X,V_Y,V_Z    |<-- | BODY VELOCITIES AT THE POSITION OF EACH      |
00522 ! |                |    | ALGAE PARTICLE                               |
00523 ! | X_A,Y_A,Z_A    |<-->| POSITION OF THE ALGAE PARTICLES BEFORE AND   |
00524 ! |                |    | AFTER DISPLACEMENT                           |
00525 ! | DX_A,DY_A,DZ_A |<-- | DISPLACEMENT OF EACH ALGAE PARTICLE          |
00526 ! | ELEM_ALG       |<-->| NUMBER OF THE ELEMENT CONTAINING THE ALGAE   |
00527 ! |                |    | PARTICLE                                     |
00528 ! | LT             | -->| TIME ITERATION OF THE SIMULATION             |
00529 ! | DALGAE         | -->| DIAMETER OF THE ALGAE PARTICLES              |
00530 ! | RALGAE         | -->| DENSITY OF THE ALGAE PARTICLES               |
00531 ! | EALGAE         | -->| THICKNESS OF THE ALGAE PARTICLES             |
00532 ! | ALGTYP         | -->| ALGAE TYPE OF THE PARTICLES                  |
00533 ! |________________|____|______________________________________________|
00534 ! MODE : -->(NON MODIFIED DATA), <--(RESULT), <-->(MODIFIED DATA)
00535 !
00536 !-----------------------------------------------------------------------
00537 !
00538 ! CALLED BY : DERIVE
00539 !
00540 ! SUBROUTINE CALLED : NONE
00541 !
00542 !***********************************************************************
00543 !
00544       INTEGER LNG,LU
00545       COMMON/INFO/LNG,LU
00546 !
00547 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00548 !
00549 ! PARAMETRES
00550       DOUBLE PRECISION,PARAMETER     :: PI=3.141592654
00551       DOUBLE PRECISION,PARAMETER     :: E=EXP(1.D0)
00552       DOUBLE PRECISION,PARAMETER     :: C0=2.1D0
00553       DOUBLE PRECISION,PARAMETER     :: RHO_F=1000.D0
00554       DOUBLE PRECISION,PARAMETER     :: NU=0.000001D0
00555 ! VARIABLES DEPENDENT ON THE ALGAE SIMULATIONS
00556       INTEGER         ,INTENT(IN)    :: NA_TOT
00557       INTEGER         ,INTENT(IN)    :: NA
00558       INTEGER         ,INTENT(IN)    :: ALGAE_START
00559       DOUBLE PRECISION,INTENT(IN)    :: DT
00560       INTEGER                        :: I_A
00561       INTEGER         ,INTENT(IN)    :: LT
00562       DOUBLE PRECISION,INTENT(IN)    :: DALGAE
00563       DOUBLE PRECISION,INTENT(IN)    :: RALGAE
00564       DOUBLE PRECISION,INTENT(IN)    :: EALGAE
00565       INTEGER         ,INTENT(IN)    :: ALGTYP
00566 ! CONSTANTS OF THE BODIES
00567       DOUBLE PRECISION               :: S
00568       DOUBLE PRECISION               :: OMEGA
00569       DOUBLE PRECISION               :: M
00570       DOUBLE PRECISION               :: MASS
00571 ! TEMPORARY VARIABLES TO CLACULATE VARIABLES FOR EACH DIRECTION
00572       INTEGER         ,INTENT(IN)    :: NDIM
00573       INTEGER                        :: I_DIM
00574       DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE :: U_I_0
00575       DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE :: U_I
00576       DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE :: V_I_0
00577       DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE :: V_I
00578       DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE :: X_I_0
00579       DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE :: X_I
00580 ! PROPERTIES OF THE BODIES
00581       DOUBLE PRECISION               :: F_A
00582       DOUBLE PRECISION               :: F_B
00583       DOUBLE PRECISION               :: TAU_PART
00584       DOUBLE PRECISION               :: FI_C
00585 ! PROPERTIES OF THE FLOW
00586       DOUBLE PRECISION               :: T_I
00587       DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE :: C_I
00588       DOUBLE PRECISION               :: B_I
00589       DOUBLE PRECISION               :: NORME_U0_V0
00590 ! MEAN FLUID VARIABLES
00591       DOUBLE PRECISION,INTENT(IN)    :: U_X_AV_0(NA_TOT)
00592       DOUBLE PRECISION,INTENT(IN)    :: U_Y_AV_0(NA_TOT)
00593       DOUBLE PRECISION,INTENT(IN)    :: U_Z_AV_0(NA_TOT)
00594       DOUBLE PRECISION,INTENT(IN)    :: U_X_AV(NA_TOT)
00595       DOUBLE PRECISION,INTENT(IN)    :: U_Y_AV(NA_TOT)
00596       DOUBLE PRECISION,INTENT(IN)    :: U_Z_AV(NA_TOT)
00597       DOUBLE PRECISION,INTENT(IN)    :: K_AV_0(NA_TOT)
00598       DOUBLE PRECISION,INTENT(IN)    :: EPS_AV_0(NA_TOT)
00599       DOUBLE PRECISION,INTENT(IN)    :: H_FLU(NA_TOT)
00600 ! VARIABLES OF THE BODIES
00601       DOUBLE PRECISION,INTENT(IN)    :: U_X_0(NA_TOT)
00602       DOUBLE PRECISION,INTENT(IN)    :: U_Y_0(NA_TOT)
00603       DOUBLE PRECISION,INTENT(IN)    :: U_Z_0(NA_TOT)
00604       DOUBLE PRECISION,INTENT(OUT)   :: U_X(NA_TOT)
00605       DOUBLE PRECISION,INTENT(OUT)   :: U_Y(NA_TOT)
00606       DOUBLE PRECISION,INTENT(OUT)   :: U_Z(NA_TOT)
00607       DOUBLE PRECISION,INTENT(IN)    :: V_X_0(NA_TOT)
00608       DOUBLE PRECISION,INTENT(IN)    :: V_Y_0(NA_TOT)
00609       DOUBLE PRECISION,INTENT(IN)    :: V_Z_0(NA_TOT)
00610       DOUBLE PRECISION,INTENT(OUT)   :: V_X(NA_TOT)
00611       DOUBLE PRECISION,INTENT(OUT)   :: V_Y(NA_TOT)
00612       DOUBLE PRECISION,INTENT(OUT)   :: V_Z(NA_TOT)
00613       DOUBLE PRECISION,INTENT(INOUT) :: X_A(NA_TOT)
00614       DOUBLE PRECISION,INTENT(INOUT) :: Y_A(NA_TOT)
00615       DOUBLE PRECISION,INTENT(INOUT) :: Z_A(NA_TOT)
00616       DOUBLE PRECISION,INTENT(OUT)   :: DX_A(NA_TOT)
00617       DOUBLE PRECISION,INTENT(OUT)   :: DY_A(NA_TOT)
00618       DOUBLE PRECISION,INTENT(OUT)   :: DZ_A(NA_TOT)
00619       INTEGER         ,INTENT(IN)    :: ELEM_ALG(NA_TOT)
00620 ! RANDOM NUMBERS
00621       DOUBLE PRECISION               :: RAND1
00622       DOUBLE PRECISION               :: XI_G_I
00623       DOUBLE PRECISION               :: RAND2
00624       DOUBLE PRECISION               :: XI_CAPG_I
00625       DOUBLE PRECISION               :: RAND3
00626       DOUBLE PRECISION               :: XI_CAPP_I
00627 ! VARIABLES USED TO FIND THE DRAG COEFFICIENT
00628       DOUBLE PRECISION               :: RE
00629       DOUBLE PRECISION               :: CD
00630       DOUBLE PRECISION               :: PHI_1
00631       DOUBLE PRECISION               :: PHI_2
00632       DOUBLE PRECISION               :: PHI_3
00633       DOUBLE PRECISION               :: PHI_4
00634 ! VARIABLES USED TO CALCULATE THE STOCHASTIC INTEGRALS
00635       DOUBLE PRECISION               :: GAMMA_I
00636       DOUBLE PRECISION               :: CAPGAMMA_I
00637       DOUBLE PRECISION               :: CAPPHI_I
00638       DOUBLE PRECISION               :: COV_G_I2
00639       DOUBLE PRECISION               :: COV_CAPG_I2
00640       DOUBLE PRECISION               :: COV_CAPP_I2
00641       DOUBLE PRECISION               :: COV_G_ICAPG_I
00642       DOUBLE PRECISION               :: COV_G_ICAPP_I
00643       DOUBLE PRECISION               :: COV_CAPG_ICAPP_I
00644       DOUBLE PRECISION               :: L11
00645       DOUBLE PRECISION               :: L21
00646       DOUBLE PRECISION               :: L22
00647       DOUBLE PRECISION               :: L31
00648       DOUBLE PRECISION               :: L32
00649       DOUBLE PRECISION               :: L33
00650 ! VARIABLES USED FOR THE EXACT INTEGRATOR
00651       DOUBLE PRECISION               :: ALPHA
00652       DOUBLE PRECISION               :: BETA
00653       DOUBLE PRECISION               :: C_CHECK
00654       DOUBLE PRECISION               :: K_CHECK
00655       DOUBLE PRECISION               :: Q_CHECK
00656       DOUBLE PRECISION               :: G_CHECK
00657 ! SAVE ALLOCATABLE VARIABLES
00658       SAVE :: U_I_0,U_I,V_I_0,V_I,X_I_0,X_I,C_I
00659 ! DAJ
00660 !     CHARACTER(LEN=100)              :: LOG_NAME
00661 ! FAJ
00662 !
00663 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00664 !
00665 ! C
00666 ! C FICHIER TXT AVEC LE JOURNAL DES ERREURS
00667 ! !       WRITE(LOG_NAME,'(A,I1,A)')'../erreur_log',RANG,'.txt'
00668 ! !       OPEN(9999,FILE=LOG_NAME)
00669 ! !       OPEN(9999,FILE='../erreur_log.txt')
00670 ! !
00671 ! !             WRITE(9999,*)LT
00672 ! !             WRITE(9999,*)PSI(1,1,:)
00673 
00674 !=======================================================================
00675 ! PREAMBLE OF THE CALCULATIONS
00676 !=======================================================================
00677 ! ALLOCATE VARIABLES
00678       IF(.NOT.ALLOCATED(U_I_0))ALLOCATE(U_I_0(NDIM))
00679       IF(.NOT.ALLOCATED(U_I))ALLOCATE(U_I(NDIM))
00680       IF(.NOT.ALLOCATED(V_I_0))ALLOCATE(V_I_0(NDIM))
00681       IF(.NOT.ALLOCATED(V_I))ALLOCATE(V_I(NDIM))
00682       IF(.NOT.ALLOCATED(X_I_0))ALLOCATE(X_I_0(NDIM))
00683       IF(.NOT.ALLOCATED(X_I))ALLOCATE(X_I(NDIM))
00684       IF(.NOT.ALLOCATED(C_I))ALLOCATE(C_I(NDIM))
00685 ! CONSTANTS
00686       IF(ALGTYP.EQ.1)THEN !SPHERE
00687         S=PI*DALGAE**2/4.D0
00688         OMEGA=PI*DALGAE**3/6.D0
00689         M=0.5D0*RHO_F*OMEGA
00690       ELSEIF(ALGTYP.EQ.2)THEN !IRIDAEA FLACCIDA (CLOSE TO ULVA)
00691         S=PI*DALGAE**2/4.D0
00692         OMEGA=S*EALGAE
00693         M=3.57D0*RHO_F*OMEGA
00694       ELSEIF(ALGTYP.EQ.3)THEN !PELVETIOPSIS LIMITATA
00695         S=PI*DALGAE**2/4.D0
00696         OMEGA=S*EALGAE
00697         M=4.64D0*RHO_F*OMEGA
00698       ELSEIF(ALGTYP.EQ.4)THEN !GIGARTINA LEPTORHYNCHOS
00699         S=PI*DALGAE**2/4.D0
00700         OMEGA=S*EALGAE
00701         M=3.26D0*RHO_F*OMEGA
00702       END IF
00703       MASS=RALGAE*OMEGA
00704       CB=6.D0*DALGAE**2*RHO_F*SQRT(PI*NU)
00705 !
00706 !=======================================================================
00707 ! START THE CALCULATIONS FOR EACH PARTICLE
00708 !=======================================================================
00709       DO I_A=1,NA
00710 !=======================================================================
00711 ! CHECK TO SEE IF THE TRANSPORT NEEDS TO BE CALCULATED
00712 !=======================================================================
00713         IF(ELEM_ALG(I_A).LE.0)GOTO 12
00714         IF(H_FLU(I_A).LT.DALGAE)GOTO 12
00715 !
00716 !=======================================================================
00717 ! REDEFINE THE PREVIOUS VARIABLES
00718 !=======================================================================
00719         U_I_0(1)=U_X_0(I_A)
00720         U_I_0(2)=U_Y_0(I_A)
00721         V_I_0(1)=V_X_0(I_A)
00722         V_I_0(2)=V_Y_0(I_A)
00723         X_I_0(1)=X_A(I_A)
00724         X_I_0(2)=Y_A(I_A)
00725         IF(NDIM.EQ.3)THEN
00726           U_I_0(3)=U_Z_0(I_A)
00727           V_I_0(3)=V_Z_0(I_A)
00728           X_I_0(3)=Z_A(I_A)
00729         END IF
00730 !=======================================================================
00731 ! REDEFINE THE FLUID PROPERTIES
00732 !=======================================================================
00733         T_I=1.D0/(0.5D0+0.75D0*C0)*K_AV_0(I_A)/EPS_AV_0(I_A)
00734 ! IN C_I IT IS ASSUMED THAT 1/RHO_F*dP/dX_i = (U_X_AV(I_A)-U_X_AV_0(I_A))/DT
00735 !    => maybe use variation in surface elevation ALONG X (I.E. HYDROSTATIC PRESSURE)
00736         C_I(1)=(U_X_AV(I_A)-U_X_AV_0(I_A))/DT+1.D0/T_I*U_X_AV_0(I_A)
00737         C_I(2)=(U_Y_AV(I_A)-U_Y_AV_0(I_A))/DT+1.D0/T_I*U_Y_AV_0(I_A)
00738         B_I=(C0*EPS_AV_0(I_A))**0.5
00739         NORME_U0_V0=SQRT((U_X_0(I_A)-V_X_0(I_A))**2+(U_Y_0(I_A)
00740      &   -V_Y_0(I_A))**2)
00741         IF(NDIM.EQ.3)THEN
00742           C_I(2)=(U_Z_AV(I_A)-U_Z_AV_0(I_A))/DT+1.D0/T_I*U_Z_AV_0(I_A)
00743           NORME_U0_V0=SQRT((U_X_0(I_A)-V_X_0(I_A))**2+(U_Y_0(I_A)
00744      &     -V_Y_0(I_A))**2+(U_Z_0(I_A)-V_Z_0(I_A))**2)
00745         END IF
00746 !=======================================================================
00747         DO I_DIM=1,NDIM
00748 !=======================================================================
00749 !=======================================================================
00750 ! REDEFINE THE BODY PROPERTIES
00751 !=======================================================================
00752           CI_BAS=0.D0
00753           F_TAIL=0.D0
00754 !
00755           IB=LT-ALGAE_START+1
00756 !
00757           IF(IB.LE.NWIN)THEN
00758             NB=IB
00759 !
00760             CI_BAS=CI_BAS+2.D0/3.D0*CB*PSI(I_A,I_DIM,NB)*SQRT(DT)
00761      &       *(3.D0*SQRT(REAL(NB))+2.D0*(REAL(NB)-1.D0)**(1.5D0)
00762      &       -2.D0*REAL(NB)**(1.5D0))
00763             DO IWIN=1,NB-1
00764               CI_BAS=CI_BAS+4.D0/3.D0*CB*SQRT(DT)*PSI(I_A,I_DIM,IWIN)
00765      &         *((REAL(IWIN)-1.D0)**(1.5D0)-2.D0*REAL(IWIN)**(1.5D0)
00766      &         +(REAL(IWIN)+1.D0)**(1.5D0))
00767             END DO
00768 !
00769             CI_BAS=CI_BAS+F_TAIL
00770 !
00771             F_A=(RHO_F*OMEGA+M+4.D0/3.D0*CB*SQRT(DT))/(MASS+M+4.D0/3.D0
00772      &        *CB*SQRT(DT))
00773 !
00774             RE=NORME_U0_V0*DALGAE/NU
00775             IF(ALGTYP.EQ.1)THEN !SPHERE
00776               IF(RE.EQ.0.D0)THEN
00777                 CD=0.D0
00778                 F_B=(RHO_F*S*24.D0*NU)/(2.D0*(MASS+M+4.D0/3.D0*CB
00779      &           *SQRT(DT))*DALGAE)
00780               ELSEIF(RE.LT.0.4)THEN
00781                 CD=24.D0/RE
00782                 F_B=(RHO_F*S*24.D0*NU)/(2.D0*(MASS+M+4.D0/3.D0*CB
00783      &           *SQRT(DT))*DALGAE)
00784               ELSEIF(RE.GT.1000000)THEN
00785                 CD=0.2
00786                 F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0*CB
00787      &           *SQRT(DT)))
00788               ELSE
00789                 PHI_1=(24.D0*RE**(-1))**10+(21.D0*RE**(-0.67))**10
00790      &           +(4.D0*RE**(-0.33))**10+0.4D0**10
00791                 PHI_2=1.D0/((0.148D0*RE**0.11)**(-10)+0.5D0**(-10))
00792                 PHI_3=(1.57D0*10.D0**8*RE**(-1.625))**10
00793                 PHI_4=1.D0/((6.D0*10.D0**(-17)*RE**2.63)**(-10)
00794      &           +0.2**(-10))
00795                 CD=(1.D0/((PHI_1+PHI_2)**(-1)+PHI_3**(-1))+PHI_4)**0.1
00796                 F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0*CB
00797      &           *SQRT(DT)))
00798               END IF
00799             ELSEIF(ALGTYP.EQ.2)THEN !IRIDAEA FLACCIDA (CLOSE TO ULVA)
00800               IF(RE.GE.14073.D0)THEN
00801                   CD=EXP(6.822121D0-0.800627D0*LOG(RE))
00802                   F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0
00803      &             *CB*SQRT(DT)))
00804               ELSE
00805                 IF(RE.EQ.0.D0)THEN
00806                   CD=0.D0
00807                   F_B=(RHO_F*S*24.D0*NU)/(2.D0*(MASS+M+4.D0/3.D0*CB
00808      &             *SQRT(DT))*DALGAE)
00809                 ELSEIF(RE.LT.0.4)THEN
00810                   CD=24.D0/RE
00811                   F_B=(RHO_F*S*24.D0*NU)/(2.D0*(MASS+M+4.D0/3.D0*CB
00812      &             *SQRT(DT))*DALGAE)
00813                 ELSEIF(RE.GT.1000000)THEN
00814                   CD=0.2
00815                   F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0
00816      &             *CB*SQRT(DT)))
00817                 ELSE
00818                   PHI_1=(24.D0*RE**(-1))**10+(21.D0*RE**(-0.67))**10
00819      &             +(4.D0*RE**(-0.33))**10+0.4D0**10
00820                   PHI_2=1.D0/((0.148D0*RE**0.11)**(-10)+0.5D0**(-10))
00821                   PHI_3=(1.57D0*10.D0**8*RE**(-1.625))**10
00822                   PHI_4=1.D0/((6.D0*10.D0**(-17)*RE**2.63)**(-10)
00823      &             +0.2**(-10))
00824                   CD=(1.D0/((PHI_1+PHI_2)**(-1)+PHI_3**(-1))+PHI_4)
00825      &             **0.1
00826                   F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0
00827      &             *CB*SQRT(DT)))
00828                 END IF
00829               END IF
00830             ELSEIF(ALGTYP.EQ.3)THEN !PELVETIOPSIS LIMITATA
00831               IF(RE.GE.28611.D0)THEN
00832                   CD=EXP(8.214783D0-0.877036D0*LOG(RE))
00833                   F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0
00834      &             *CB*SQRT(DT)))
00835               ELSE
00836                 IF(RE.EQ.0.D0)THEN
00837                   CD=0.D0
00838                   F_B=(RHO_F*S*24.D0*NU)/(2.D0*(MASS+M+4.D0/3.D0*CB
00839      &             *SQRT(DT))*DALGAE)
00840                 ELSEIF(RE.LT.0.4)THEN
00841                   CD=24.D0/RE
00842                   F_B=(RHO_F*S*24.D0*NU)/(2.D0*(MASS+M+4.D0/3.D0*CB
00843      &             *SQRT(DT))*DALGAE)
00844                 ELSEIF(RE.GT.1000000)THEN
00845                   CD=0.2
00846                   F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0
00847      &             *CB*SQRT(DT)))
00848                 ELSE
00849                   PHI_1=(24.D0*RE**(-1))**10+(21.D0*RE**(-0.67))**10
00850      &             +(4.D0*RE**(-0.33))**10+0.4D0**10
00851                   PHI_2=1.D0/((0.148D0*RE**0.11)**(-10)+0.5D0**(-10))
00852                   PHI_3=(1.57D0*10.D0**8*RE**(-1.625))**10
00853                   PHI_4=1.D0/((6.D0*10.D0**(-17)*RE**2.63)**(-10)
00854      &             +0.2**(-10))
00855                   CD=(1.D0/((PHI_1+PHI_2)**(-1)+PHI_3**(-1))+PHI_4)
00856      &             **0.1
00857                   F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0
00858      &             *CB*SQRT(DT)))
00859                 END IF
00860               END IF
00861             ELSEIF(ALGTYP.EQ.4)THEN !GIGARTINA LEPTORHYNCHOS
00862               IF(RE.GE.17981.D0)THEN
00863                   CD=EXP(6.773712D0-0.774252D0*LOG(RE))
00864                   F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0
00865      &             *CB*SQRT(DT)))
00866               ELSE
00867                 IF(RE.EQ.0.D0)THEN
00868                   CD=0.D0
00869                   F_B=(RHO_F*S*24.D0*NU)/(2.D0*(MASS+M+4.D0/3.D0*CB
00870      &             *SQRT(DT))*DALGAE)
00871                 ELSEIF(RE.LT.0.4)THEN
00872                   CD=24.D0/RE
00873                   F_B=(RHO_F*S*24.D0*NU)/(2.D0*(MASS+M+4.D0/3.D0*CB
00874      &             *SQRT(DT))*DALGAE)
00875                 ELSEIF(RE.GT.1000000)THEN
00876                   CD=0.2
00877                   F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0
00878      &             *CB*SQRT(DT)))
00879                 ELSE
00880                   PHI_1=(24.D0*RE**(-1))**10+(21.D0*RE**(-0.67))**10
00881      &             +(4.D0*RE**(-0.33))**10+0.4D0**10
00882                   PHI_2=1.D0/((0.148D0*RE**0.11)**(-10)+0.5D0**(-10))
00883                   PHI_3=(1.57D0*10.D0**8*RE**(-1.625))**10
00884                   PHI_4=1.D0/((6.D0*10.D0**(-17)*RE**2.63)**(-10)
00885      &             +0.2**(-10))
00886                   CD=(1.D0/((PHI_1+PHI_2)**(-1)+PHI_3**(-1))+PHI_4)
00887      &             **0.1
00888                   F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0
00889      &             *CB*SQRT(DT)))
00890                 END IF
00891               END IF
00892             END IF
00893 !
00894             TAU_PART=1.D0/F_B
00895 !
00896             FI_C=CI_BAS/(MASS+M+4.D0/3.D0*CB*SQRT(DT))
00897           ELSE
00898             NB=NWIN
00899 !
00900             DO IP=1,NP
00901               PHI_PLUS=(EXP(DT/(2.D0*T_TIL_P(IP)*TWIN))-1.D0)/(DT
00902      &         /(2.D0*T_TIL_P(IP)*TWIN))
00903               PHI_MOINS=(EXP(-DT/(2.D0*T_TIL_P(IP)*TWIN))-1.D0)/(-DT
00904      &         /(2.D0*T_TIL_P(IP)*TWIN))
00905 !
00906               FI_P(I_A,I_DIM,IP)=2.D0*CB*SQRT(E*T_TIL_P(IP)*TWIN)
00907      &         *EXP(-1.D0/(2.D0*T_TIL_P(IP)))*(PSI(I_A,I_DIM,NB)*(1.D0
00908      &         -PHI_MOINS)+PSI(I_A,I_DIM,NB+1)*EXP(-DT/(2.D0*T_TIL_P(IP)
00909      &         *TWIN))*(PHI_PLUS-1))+EXP(-DT/(2.D0*T_TIL_P(IP)*TWIN))
00910      &         *FI_P(I_A,I_DIM,IP)
00911 !
00912               F_TAIL=F_TAIL+A_P(IP)*FI_P(I_A,I_DIM,IP)
00913             END DO
00914 !
00915             CI_BAS=CI_BAS+2.D0/3.D0*CB*PSI(I_A,I_DIM,NB)*SQRT(DT)
00916      &       *(3.D0*SQRT(REAL(NB))+2.D0*(REAL(NB)-1.D0)**(1.5D0)
00917      &       -2.D0*REAL(NB)**(1.5D0))
00918             DO IWIN=1,NB-1
00919               CI_BAS=CI_BAS+4.D0/3.D0*CB*SQRT(DT)*PSI(I_A,I_DIM,IWIN)
00920      &         *((REAL(IWIN)-1.D0)**(1.5D0)-2.D0*REAL(IWIN)**(1.5D0)
00921      &         +(REAL(IWIN)+1.D0)**(1.5D0))
00922             END DO
00923             CI_BAS=CI_BAS+F_TAIL
00924 !
00925             F_A=(RHO_F*OMEGA+M+4.D0/3.D0*CB*SQRT(DT))/(MASS+M+4.D0/3.D0
00926      &        *CB*SQRT(DT))
00927 !
00928             RE=NORME_U0_V0*DALGAE/NU
00929             IF(ALGTYP.EQ.1)THEN !SPHERE
00930               IF(RE.EQ.0.D0)THEN
00931                 CD=0.D0
00932                 F_B=(RHO_F*S*24.D0*NU)/(2.D0*(MASS+M+4.D0/3.D0*CB
00933      &           *SQRT(DT))*DALGAE)
00934               ELSEIF(RE.LT.0.4)THEN
00935                 CD=24.D0/RE
00936                 F_B=(RHO_F*S*24.D0*NU)/(2.D0*(MASS+M+4.D0/3.D0*CB
00937      &           *SQRT(DT))*DALGAE)
00938               ELSEIF(RE.GT.1000000)THEN
00939                 CD=0.2
00940                 F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0*CB
00941      &           *SQRT(DT)))
00942               ELSE
00943                 PHI_1=(24.D0*RE**(-1))**10+(21.D0*RE**(-0.67))**10
00944      &           +(4.D0*RE**(-0.33))**10+0.4D0**10
00945                 PHI_2=1.D0/((0.148D0*RE**0.11)**(-10)+0.5D0**(-10))
00946                 PHI_3=(1.57D0*10.D0**8*RE**(-1.625))**10
00947                 PHI_4=1.D0/((6.D0*10.D0**(-17)*RE**2.63)**(-10)
00948      &           +0.2**(-10))
00949                 CD=(1.D0/((PHI_1+PHI_2)**(-1)+PHI_3**(-1))+PHI_4)**0.1
00950                 F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0*CB
00951      &           *SQRT(DT)))
00952               END IF
00953             ELSEIF(ALGTYP.EQ.2)THEN !IRIDAEA FLACCIDA (CLOSE TO ULVA)
00954               IF(RE.GE.14073.D0)THEN
00955                   CD=EXP(6.822121D0-0.800627D0*LOG(RE))
00956                   F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0
00957      &             *CB*SQRT(DT)))
00958               ELSE
00959                 IF(RE.EQ.0.D0)THEN
00960                   CD=0.D0
00961                   F_B=(RHO_F*S*24.D0*NU)/(2.D0*(MASS+M+4.D0/3.D0*CB
00962      &             *SQRT(DT))*DALGAE)
00963                 ELSEIF(RE.LT.0.4)THEN
00964                   CD=24.D0/RE
00965                   F_B=(RHO_F*S*24.D0*NU)/(2.D0*(MASS+M+4.D0/3.D0*CB
00966      &             *SQRT(DT))*DALGAE)
00967                 ELSEIF(RE.GT.1000000)THEN
00968                   CD=0.2
00969                   F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0
00970      &             *CB*SQRT(DT)))
00971                 ELSE
00972                   PHI_1=(24.D0*RE**(-1))**10+(21.D0*RE**(-0.67))**10
00973      &             +(4.D0*RE**(-0.33))**10+0.4D0**10
00974                   PHI_2=1.D0/((0.148D0*RE**0.11)**(-10)+0.5D0**(-10))
00975                   PHI_3=(1.57D0*10.D0**8*RE**(-1.625))**10
00976                   PHI_4=1.D0/((6.D0*10.D0**(-17)*RE**2.63)**(-10)
00977      &             +0.2**(-10))
00978                   CD=(1.D0/((PHI_1+PHI_2)**(-1)+PHI_3**(-1))+PHI_4)
00979      &             **0.1
00980                   F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0
00981      &             *CB*SQRT(DT)))
00982                 END IF
00983               END IF
00984             ELSEIF(ALGTYP.EQ.3)THEN !PELVETIOPSIS LIMITATA
00985               IF(RE.GE.28611.D0)THEN
00986                   CD=EXP(8.214783D0-0.877036D0*LOG(RE))
00987                   F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0
00988      &             *CB*SQRT(DT)))
00989               ELSE
00990                 IF(RE.EQ.0.D0)THEN
00991                   CD=0.D0
00992                   F_B=(RHO_F*S*24.D0*NU)/(2.D0*(MASS+M+4.D0/3.D0*CB
00993      &             *SQRT(DT))*DALGAE)
00994                 ELSEIF(RE.LT.0.4)THEN
00995                   CD=24.D0/RE
00996                   F_B=(RHO_F*S*24.D0*NU)/(2.D0*(MASS+M+4.D0/3.D0*CB
00997      &             *SQRT(DT))*DALGAE)
00998                 ELSEIF(RE.GT.1000000)THEN
00999                   CD=0.2
01000                   F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0
01001      &             *CB*SQRT(DT)))
01002                 ELSE
01003                   PHI_1=(24.D0*RE**(-1))**10+(21.D0*RE**(-0.67))**10
01004      &             +(4.D0*RE**(-0.33))**10+0.4D0**10
01005                   PHI_2=1.D0/((0.148D0*RE**0.11)**(-10)+0.5D0**(-10))
01006                   PHI_3=(1.57D0*10.D0**8*RE**(-1.625))**10
01007                   PHI_4=1.D0/((6.D0*10.D0**(-17)*RE**2.63)**(-10)
01008      &             +0.2**(-10))
01009                   CD=(1.D0/((PHI_1+PHI_2)**(-1)+PHI_3**(-1))+PHI_4)
01010      &             **0.1
01011                   F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0
01012      &             *CB*SQRT(DT)))
01013                 END IF
01014               END IF
01015             ELSEIF(ALGTYP.EQ.4)THEN !GIGARTINA LEPTORHYNCHOS
01016               IF(RE.GE.17981.D0)THEN
01017                   CD=EXP(6.773712D0-0.774252D0*LOG(RE))
01018                   F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0
01019      &             *CB*SQRT(DT)))
01020               ELSE
01021                 IF(RE.EQ.0.D0)THEN
01022                   CD=0.D0
01023                   F_B=(RHO_F*S*24.D0*NU)/(2.D0*(MASS+M+4.D0/3.D0*CB
01024      &             *SQRT(DT))*DALGAE)
01025                 ELSEIF(RE.LT.0.4)THEN
01026                   CD=24.D0/RE
01027                   F_B=(RHO_F*S*24.D0*NU)/(2.D0*(MASS+M+4.D0/3.D0*CB
01028      &             *SQRT(DT))*DALGAE)
01029                 ELSEIF(RE.GT.1000000)THEN
01030                   CD=0.2
01031                   F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0
01032      &             *CB*SQRT(DT)))
01033                 ELSE
01034                   PHI_1=(24.D0*RE**(-1))**10+(21.D0*RE**(-0.67))**10
01035      &             +(4.D0*RE**(-0.33))**10+0.4D0**10
01036                   PHI_2=1.D0/((0.148D0*RE**0.11)**(-10)+0.5D0**(-10))
01037                   PHI_3=(1.57D0*10.D0**8*RE**(-1.625))**10
01038                   PHI_4=1.D0/((6.D0*10.D0**(-17)*RE**2.63)**(-10)
01039      &             +0.2**(-10))
01040                   CD=(1.D0/((PHI_1+PHI_2)**(-1)+PHI_3**(-1))+PHI_4)
01041      &             **0.1
01042                   F_B=(RHO_F*S*CD*NORME_U0_V0)/(2.D0*(MASS+M+4.D0/3.D0
01043      &             *CB*SQRT(DT)))
01044                 END IF
01045               END IF
01046             END IF
01047 !
01048             TAU_PART=1.D0/F_B
01049 !
01050             FI_C=CI_BAS/(MASS+M+4.D0/3.D0*CB*SQRT(DT))
01051 !
01052           END IF
01053 ! T_I SHOULD NEVER EQUAL TAU_PART, PEFORM A CHECK JUST IN CASE
01054           IF(T_I-TAU_PART.EQ.0)THEN
01055             WRITE(LU,*)''
01056             WRITE(LU,*)'                           **************'
01057             WRITE(LU,*)'                           *T_I=TAU_PART*'
01058             WRITE(LU,*)'                           **************'
01059             WRITE(LU,*)''
01060             CALL PLANTE(1)
01061             STOP
01062           END IF
01063 !=======================================================================
01064 ! RANDOM NUMBERS
01065 !=======================================================================
01066           CALL RANDOM_NUMBER(RAND1)
01067           XI_G_I=(RAND1*2.D0-1.D0)*(12.D0**0.5)/2.D0
01068           CALL RANDOM_NUMBER(RAND2)
01069           XI_CAPG_I=(RAND2*2.D0-1.D0)*(12.D0**0.5)/2.D0
01070           CALL RANDOM_NUMBER(RAND3)
01071           XI_CAPP_I=(RAND3*2.D0-1.D0)*(12.D0**0.5)/2.D0
01072 !=======================================================================
01073 ! FLUID VELOCITIES
01074 !=======================================================================
01075 ! VARIABLES OF THE EXACT INTEGRATOR
01076           ALPHA=EXP(-DT/T_I)
01077 !
01078           COV_G_I2=(1.D0-ALPHA**2)*B_I**2*T_I/2.D0
01079 !
01080           L11=SQRT(COV_G_I2)
01081           GAMMA_I=L11*XI_G_I
01082 !FLUID VELOCITY
01083           U_I(I_DIM)=ALPHA*U_I_0(I_DIM)+(1.D0-ALPHA)*C_I(I_DIM)*T_I
01084      &               +GAMMA_I
01085 !
01086 !=======================================================================
01087 ! BODY VELOCITIES
01088 !=======================================================================
01089 ! VARIABLES OF THE EXACT INTEGRATOR
01090           BETA=EXP(-DT/TAU_PART)
01091           C_CHECK=(T_I-F_A*TAU_PART)/(T_I-TAU_PART)
01092           K_CHECK=F_A/C_CHECK-1.D0
01093           Q_CHECK=K_CHECK*T_I*TAU_PART/(T_I+TAU_PART)
01094 !
01095           COV_CAPG_I2=(B_I*C_CHECK)**2*((1.D0-ALPHA**2)*T_I/2.D0+(1.D0
01096      &     -BETA**2)*K_CHECK**2*TAU_PART/2.D0+2.D0*(1.D0-ALPHA*BETA)
01097      &     *Q_CHECK)
01098           COV_G_ICAPG_I=B_I**2*C_CHECK*((1.D0-ALPHA**2)*T_I/2.D0
01099      &     +(1.D0-ALPHA*BETA)*Q_CHECK)
01100 !
01101           L21=COV_G_ICAPG_I/L11
01102 ! TO STOP ROUND OFF ERRORS CREATING NaN
01103           IF(COV_CAPG_I2.GE.L21**2)THEN
01104             L22=SQRT(COV_CAPG_I2-L21**2)
01105           ELSE
01106             L22=SQRT(L21**2-COV_CAPG_I2)
01107           END IF
01108           IF(L22.EQ.0.D0)THEN
01109              L22=1.D-12
01110           END IF
01111 !
01112           CAPGAMMA_I=L21*XI_G_I+L22*XI_CAPG_I
01113 ! BODY VELOCITY
01114           V_I(I_DIM)=BETA*V_I_0(I_DIM)+(1.D0-BETA)*(C_I(I_DIM)*T_I+FI_C)
01115      &               +(ALPHA-BETA)*C_CHECK*(U_I_0(I_DIM)-C_I(I_DIM)*T_I)
01116      &               +CAPGAMMA_I
01117 !
01118 !=======================================================================
01119 ! ALGAE POSITIONS
01120 !=======================================================================
01121 ! VARIABLES OF THE EXACT INTEGRATOR
01122           G_CHECK=T_I+K_CHECK*TAU_PART
01123 !
01124           COV_CAPP_I2=(B_I*C_CHECK)**2*(G_CHECK**2*DT+(1.D0-ALPHA**2)
01125      &     *T_I**3/2.D0+(1.D0-BETA**2)*K_CHECK**2*TAU_PART**3/2.D0
01126      &     -2.D0*G_CHECK*((1.D0-ALPHA)*T_I**2+(1.D0-BETA)*K_CHECK
01127      &     *TAU_PART**2)+2.D0*(1.D0-ALPHA*BETA)*Q_CHECK*T_I*TAU_PART)
01128 !
01129           COV_G_ICAPP_I=B_I**2*C_CHECK*((1.D0-ALPHA)*G_CHECK*T_I-(1.D0
01130      &     -ALPHA**2)*T_I**2/2.D0-(1.D0-ALPHA*BETA)*Q_CHECK*TAU_PART)
01131 !
01132           COV_CAPG_ICAPP_I=(B_I*C_CHECK)**2*(((1.D0-ALPHA)*T_I+(1.D0
01133      &     -BETA)*K_CHECK*TAU_PART)*G_CHECK-(1.D0-ALPHA**2)*T_I**2
01134      &     /2.D0-(1.D0-BETA**2)*K_CHECK**2*TAU_PART**2/2.D0-(1.D0
01135      &     -ALPHA*BETA)*K_CHECK*T_I*TAU_PART)
01136 !
01137           L31=COV_G_ICAPP_I/L11
01138           L32=(COV_CAPG_ICAPP_I-L21*L31)/L22
01139 ! TO STOP ROUND OFF ERRORS CREATING NaN
01140           IF(COV_CAPP_I2.GE.L31**2+L32**2)THEN
01141             L33=SQRT(COV_CAPP_I2-L31**2-L32**2)
01142           ELSE
01143             L33=SQRT(ABS(COV_CAPP_I2-L31**2-L32**2))
01144           END IF
01145 !
01146           CAPPHI_I=L31*XI_G_I+L32*XI_CAPG_I+L33*XI_CAPP_I
01147 !ALGAE POSITION
01148           X_I(I_DIM)=X_I_0(I_DIM)+(1.D0-BETA)*TAU_PART*V_I_0(I_DIM)+(DT
01149      &     -(1.D0-BETA)*TAU_PART)*(C_I(I_DIM)*T_I+FI_C*TAU_PART)+C_CHECK
01150      &     *(U_I_0(I_DIM)-C_I(I_DIM)*T_I)*((1.D0-ALPHA)*T_I-(1.D0-BETA)
01151      &     *TAU_PART)+CAPPHI_I
01152 !
01153 !=======================================================================
01154 ! CHECK FOR POSSIBLE ERRORS
01155 !=======================================================================
01156           IF(ABS(X_I(I_DIM)).GT.10.D0**10)THEN
01157             WRITE(LU,*)''
01158             WRITE(LU,*)'                           **************'
01159             WRITE(LU,*)'                           *POSITION INF*'
01160             WRITE(LU,*)'                           **************'
01161             WRITE(LU,*)''
01162             CALL PLANTE(1)
01163             STOP
01164           END IF
01165 !
01166 !=======================================================================
01167 ! REDEFINE PSI
01168 !=======================================================================
01169           DO IWIN=1,NWIN
01170              PSI(I_A,I_DIM,NWIN+1-IWIN+1)=PSI(I_A,I_DIM,NWIN+1-IWIN)
01171           END DO
01172           PSI(I_A,I_DIM,1)=((U_I(I_DIM)-V_I(I_DIM))-(U_I_0(I_DIM)
01173      &     -V_I_0(I_DIM)))/DT
01174 !=======================================================================
01175         END DO
01176 !=======================================================================
01177 !=======================================================================
01178 ! REDEFINIR THE CURRENT VARIABLES VARIABLES ACTUELLES
01179 !=======================================================================
01180         U_X(I_A)=U_I(1)
01181         U_Y(I_A)=U_I(2)
01182         V_X(I_A)=V_I(1)
01183         V_Y(I_A)=V_I(2)
01184         X_A(I_A)=X_I(1)
01185         Y_A(I_A)=X_I(2)
01186 !
01187         DX_A(I_A)=X_I(1)-X_I_0(1)
01188         DY_A(I_A)=X_I(2)-X_I_0(2)
01189         IF(NDIM.EQ.3)THEN
01190           Z_A(I_A)=X_I(3)
01191           DZ_A(I_A)=X_I(3)-X_I_0(3)
01192         END IF
01193 
01194 !
01195         GOTO 13
01196 12      CONTINUE
01197 ! IF THE PARTICLES ARE NOT TRANSPORTED
01198         DX_A(I_A)=0.D0
01199         DY_A(I_A)=0.D0
01200         DZ_A(I_A)=0.D0
01201         V_X(I_A)=0.D0
01202         V_Y(I_A)=0.D0
01203         V_Z(I_A)=0.D0
01204         U_X(I_A)=0.D0
01205         U_Y(I_A)=0.D0
01206         U_Z(I_A)=0.D0
01207 !
01208         DO I_DIM=1,NDIM
01209           DO IWIN=1,NWIN+1
01210             PSI(I_A,I_DIM,IWIN)=0.D0
01211           END DO
01212         END DO
01213 !
01214 13      CONTINUE
01215       END DO
01216 !
01217       RETURN
01218       END SUBROUTINE DISP_ALGAE
01219 !
01220       END MODULE ALGAE_TRANSP

Generated on Fri Aug 31 2013 18:12:58 by S.E.Bourban (HRW) using doxygen 1.7.0