The TELEMAC-MASCARET system  trunk
gaia_init.F
Go to the documentation of this file.
1 ! ********************
2  SUBROUTINE gaia_init
3 ! ********************
4 !
5  &(grafcount,listcount,telnit,u_tel,v_tel,h_tel,zf_tel,uetcar,
6  & deltar,cf_tel,ks_tel,code,u3d,v3d,t_tel,dt_tel,charr_tel,
7  & susp_tel,thetaw_tel,hw_tel,tw_tel,uw_tel,yagout,xmve_tel,
8  & grav_tel)
9 !
10 !***********************************************************************
11 ! GAIA
12 !***********************************************************************
13 !
15 !
16 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
42 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
43 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
44 !
45  USE bief
48  USE interface_gaia, ex_gaia_init => gaia_init
51 !
52  USE interface_parallel, ONLY : p_max,p_sum
53  IMPLICIT NONE
54 !
55 !!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
56 !
57  INTEGER, INTENT(IN) :: GRAFCOUNT,LISTCOUNT
58  INTEGER, INTENT(IN) :: TELNIT
59  CHARACTER(LEN=24), INTENT(IN) :: CODE
60  TYPE(bief_obj), INTENT(IN) :: U_TEL,V_TEL,H_TEL
61  TYPE(bief_obj), INTENT(INOUT) :: ZF_TEL,UETCAR,KS_TEL
62  TYPE(bief_obj), INTENT(IN) :: DELTAR
63  TYPE(bief_obj), INTENT(IN) :: U3D,V3D
64  TYPE(bief_obj), INTENT(INOUT) :: CF_TEL
65  DOUBLE PRECISION, INTENT(IN) :: T_TEL,XMVE_TEL,GRAV_TEL
66  LOGICAL, INTENT(INOUT) :: CHARR_TEL,SUSP_TEL
67  DOUBLE PRECISION, INTENT(IN) :: DT_TEL
68  TYPE(bief_obj), INTENT(IN) :: THETAW_TEL,HW_TEL,TW_TEL
69  TYPE(bief_obj), INTENT(IN) :: UW_TEL
70  LOGICAL, INTENT(IN) :: YAGOUT
71 !
72 !!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
73  INTEGER :: I,K,IELEB,KP1,IVAR
74  INTEGER :: NUMENX,IERR
75  INTEGER :: TROUVE(maxvar+10)
76  INTEGER :: ISAND,IPOIN,IMUD,ICLA
77  DOUBLE PRECISION :: AT !only used for wave
78 !!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
79  INTEGER, ALLOCATABLE :: NULLT(:)
80  DOUBLE PRECISION, ALLOCATABLE :: NULLD(:)
81 !======================================================================!
82 !
83 #if defined COMPAD
84  CALL ad_gaia_begin
86 #endif
87  ! Set to 1 for Telemac3d
88  hind_type = 1
89 ! INITIALISES THE CONSTANTS
90 !
92 !
93 ! READS THE WAVE DATA IN THE HYDRODYNAMIC FILE
94 !
95  IF(houle.AND.gai_files(gaicou)%NAME(1:1).EQ.' ') THEN
96  alire(12)=1
97  alire(13)=1
98  alire(14)=1
99  alire(22)=1
100  ENDIF
101 !
102 ! READS THE SEDIMENTOLOGICAL DATA IN THE CONTINUATION FILE
103 ! IN CASE OF COMPUTATION CONTINUED
104 !
105  IF(debu) THEN
106  alire(5)=1 ! bottom
107  alire(9)=1 ! rigid bed
108  alire(15)=1
109  alire(16)=1
110  alire(17)=1
111  alire(18)=1
112  alire(19)=1
113  alire(20)=1
114  alire(21)=1
115 ! READS REFERENCE LEVEL FOR NESTOR
116  IF(nestor) alire(23)=1 ! Nestor
117 ! READS RATIO_SAND FROM THE PREVIOUS COMPUTATION FILE
118  DO i=1,nsand*nomblay
119  alire(nvar_ratios+i)=1
120  ENDDO
121 ! READS RATIO_MUD FROM THE PREVIOUS COMPUTATION FILE
122  DO i=1,nmud*nomblay
123  alire(nvar_ratiom+i)=1
124  ENDDO
125 ! READS THE LAYER THICKNESSES
126  DO i=1,nomblay
127  alire(nvar_laythi+i)=1
128  ENDDO
129 ! READS MASS_S FROM THE PREVIOUS COMPUTATION FILE
130  DO i=1,nsand*nomblay
131  alire(nvar_mass_s+i)=1
132  ENDDO
133 ! READS MASS_M FROM THE PREVIOUS COMPUTATION FILE
134  DO i=1,nmud*nomblay
135  alire(nvar_mass_m+i)=1
136  ENDDO
137 ! READS LAYER CONCENTRATION, CRITICAL SHEAR STRESS,
138 ! PARTHENIADES CONSTANT EROSION RATE, MASS TRANSFER RATE FROM
139 ! THE PREVIOUS COMPUTATION FILE
140  DO i=1,nomblay
141  alire(nvar_layconc+i)=1
142  alire(nvar_tocemud+i)=1
143  alire(nvar_parthe+i)=1
144  alire(nvar_mtrans+i)=1
145  ENDDO
146 ! DIFFERENTIATED VARIABLES
147 ! FOR READING GRADIENTS IN SELAFIN FILES
148 !
149  IF(nadvar.GT.0) THEN
150  DO ivar=1,nadvar
151 ! SEE NOMVAR_GAIA
152  alire(nvar_advar+ivar) = 1
153  ENDDO
154  ENDIF
155 !
156  ENDIF
157 !
158 ! INITIALISES (SETS TO 0) THE ARRAYS
159 !
160  CALL init_zero_gaia
161 
162 ! -- INITIALISES THE REFERENCE LEVEL ZRL FOR NESTOR
163 ! IT MUST BE EXACLY 123456789.0D0
164  zrl%R(:) = 123456789.0d0
165 !
166 !
167 ! DISCRETISATION : LINEAR FOR THE TIME BEING
168 !
169 ! IELMT HARD-CODED IN LECDON
170 !
172  nelmax = nelem
173 !
174 !=======================================================================
175 !
176 ! READS, PREPARES AND CONTROLS THE DATA
177 !
178 !=======================================================================
179 !
180 ! READS THE BOUNDARY CONDITIONS AND INDICES FOR THE BOUNDARY NODES
181 ! EBOR IS READ HERE FOR THE FIRST CLASS ONLY
182 ! SEE CONLIT_GAIA FOR OTHER CLASSES
183 !
184  IF(debug.GT.0) WRITE(lu,*) 'APPEL DE LECLIM'
185  ALLOCATE(nullt(nptfr),stat=ierr)
186  CALL check_allocate(ierr,'NULLT')
187  ALLOCATE(nulld(nptfr*2),stat=ierr)
188  CALL check_allocate(ierr,'NULLD')
189  CALL leclim(lihbor%I,liqbor%I,nullt,liebor%I,q2bor%R,nulld,
190  & nulld,ebor%ADR(1)%P%R,nulld,afbor%R,bfbor%R,
191  & mesh%NPTFR,'GAI',.true.,
192  & gai_files(gaigeo)%FMT,gai_files(gaigeo)%LU,
193  & kent,kent,-1,-1,-1,-1,
195  DEALLOCATE(nullt)
196  DEALLOCATE(nulld)
197  IF(debug.GT.0) WRITE(lu,*) 'END_LECLIM'
198 !
199 !-----------------------------------------------------------------------
200 !
201 ! COMPLEMENTS THE DATA STRUCTURE FOR BIEF
202 !
203  IF(debug.GT.0) WRITE(lu,*) 'INBIEF'
204  CALL inbief(it1%I,klog,it2,it3,it4,lvmac,ielmx,
205  & 0.d0,spheri,mesh,t1,t2,optass,produc,equa)
206 !
207  IF(debug.GT.0) WRITE(lu,*) 'END_INBIEF'
208 !
209 !-----------------------------------------------------------------------
210 !
211 ! LOCATES THE BOUNDARIES
212 !
213  IF(ncsize.GT.1) THEN
214  nfrliq=0
215  DO i=1,nptfr
216  nfrliq=max(nfrliq,numliq%I(i))
217  ENDDO
219  WRITE(lu,*) ' '
220  WRITE(lu,*) 'LIQUID BOUNDARIES:',nfrliq
221 ! NFRLIQ CHECKED (NFRSOL NOT USED IN PARALLEL)
222  IF(nfrliq.GT.maxfro) THEN
223  WRITE(lu,*) 'INCREASE THE MAXIMUM NUMBER OF BOUNDARIES'
224  WRITE(lu,*) 'CURRENTLY AT ',maxfro
225  WRITE(lu,*) 'TO THE VALUE ',nfrliq
226  CALL plante(1)
227  stop
228  ENDIF
229  ELSE
230  IF(debug.GT.0) WRITE(lu,*) 'APPEL DE FRONT2'
231  CALL front2(nfrliq,
232  & liebor%I,liebor%I,
233  & mesh%X%R,mesh%Y%R,mesh%NBOR%I,mesh%KP1BOR%I,
234  & it1%I,npoin,nptfr,klog,.true.,numliq%I,maxfro)
235  IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE FRONT2'
236  ENDIF
237  IF(nfrliq.GT.maxfro) THEN
238  WRITE(lu,*) 'FRONT2: SIZE OF ARRAYS EXCEEDED'
239  WRITE(lu,*) ' INCREASE THE KEYWORD'
240  WRITE(lu,*) ' MAXIMUM NUMBER OF BOUNDARIES'
241  WRITE(lu,*) ' IN THE CALLING PROGRAM'
242  WRITE(lu,*) ' THE CURRENT VALUE IS ',maxfro
243  WRITE(lu,*) ' THE VALUE SHOULD BE ',nfrliq
244  CALL plante(1)
245  stop
246  ENDIF
247 !
248 !-----------------------------------------------------------------------
249 ! LOOKS FOR RADSEC FOR SECONDARY CURRENTS IN THE GEOMETRY FILE :
250 !-----------------------------------------------------------------------
251 !
252 !
253 ! RADSEC IS LANGUAGE INDEPENDENT
254 ! CALL FIND_VARIABLE(RADSEC,'RADSEC ',GAI_FILES(GAIGEO)%LU,
255 ! & GAI_FILES(GAIGEO)%FMT,W,OKRADSEC,TIME=BID)
256 ! Initialising
257  radsec%R=0.d0
259  & 'RADSEC ', radsec%R,npoin,ierr)
260 !
261 !-----------------------------------------------------------------------
262 ! LOOKS FOR RFERENCE LEVEL FOR NESTOR IN THE GEOMETRY FILE :
263 !-----------------------------------------------------------------------
264 ! REFERENCE LEVEL IS LANGUAGE INDEPENDENT
265 ! Initialising
267  & 'REFERENCE LEVEL ', zrl%R,npoin,ierr)
268  IF( ierr == 0 ) THEN
269  WRITE(6,*)'-------------------------------------------'
270  WRITE(6,*)
271  & 'Found REFERENCE LEVEL (ZRL) for Nestor in Sis Geometry file'
272  WRITE(6,*)' max value of ZRL = ',maxval(zrl%R)
273  WRITE(6,*)' min value of ZRL = ',minval(zrl%R)
274  WRITE(6,*)'-------------------------------------------'
275  ENDIF
276 !
277 !-----------------------------------------------------------------------
278 !-----------------------------------------------------------------------
279 !
280 ! PREPARES THE RESULTS FILE (OPTIONAL)
281 !
282 ! STANDARD SELAFIN FORMAT
283 !
284  IF(debug.GT.0) WRITE(lu,*) 'WRITE_HEADER'
285 ! CREATES DATA FILE USING A GIVEN FILE FORMAT : FORMAT_RES
286 ! THE DATA ARE CREATED IN THE FILE: GAIRES, AND ARE
287 ! CHARACTERISED BY A TITLE AND NAME OF OUTPUT VARIABLES
288 ! CONTAINED IN THE FILE.
289  CALL write_header(gai_files(gaires)%FMT, ! RESULTS FILE FORMAT
290  & gai_files(gaires)%LU, ! LU FOR RESULTS FILE
291  & titca, ! TITLE
292  & maxvar, ! MAX NUMBER OF OUTPUT VARIABLES
293  & texte, ! NAMES OF OUTPUT VARIABLES
294  & sorleo) ! PRINT TO FILE OR NOT
295  IF(debug.GT.0) WRITE(lu,*) 'END_WRITE_HEADER'
296  IF(debug.GT.0) WRITE(lu,*) 'WRITE_MESH'
297 ! WRITES THE MESH IN THE OUTPUT FILE :
298 ! IN PARALLEL, REQUIRES NCSIZE AND NPTIR.
299 ! THE REST OF THE INFORMATION IS IN MESH.
300 ! ALSO WRITES : START DATE/TIME AND COORDINATES OF THE
301 ! ORIGIN.
302  CALL write_mesh(gai_files(gaires)%FMT, ! RESULTS FILE FORMAT
303  & gai_files(gaires)%LU, ! LU FOR RESULTS FILE
304  & mesh,
305  & 1, ! NUMBER OF PLANES /NA/
306  & mardat, ! START DATE
307  & martim, ! START TIME
308  & t1,t2,
309  & ncsize.GT.1,nptir,
310  & ngeo=gai_files(gaigeo)%LU,
311  & geoformat=gai_files(gaigeo)%FMT)
312  IF(debug.GT.0) WRITE(lu,*) 'END_WRITE_MESH'
313 !
314 ! FILLS IN MASKEL BY DEFAULT
315 !
316  IF(msk) CALL os ('X=C ', x=maskel, c=1.d0)
317 !
318 ! BUILDING THE MASK FOR LIQUID BOUNDARIES
319 ! A SEGMENT IS LIQUID IF BOTH ENDS ARE NOT SOLID
320 !
321  DO ieleb = 1, mesh%NELEB
322  k=mesh%IKLBOR%I(ieleb)
323  kp1=mesh%IKLBOR%I(ieleb+mesh%NELEBX)
324  IF(liebor%I(k).NE.2.AND.liebor%I(kp1).NE.2) THEN
325  mask%R(ieleb) = 1.d0
326  ELSE
327  mask%R(ieleb) = 0.d0
328  ENDIF
329  ENDDO
330 !
331 !=======================================================================
332 !
333 ! INITIALISES
334 !
335 !=======================================================================
336 !
337  pass = .true.
338 !
339  second_susp_step = .false.
340  massnestor = 0.d0
341  IF(bilma) THEN
342  summcumucl= 0.d0
343  sumbedload_b = 0.d0
344  sumrmascl = 0.d0
345  sum_deposition = 0.d0
346  sum_erosion = 0.d0
347  summassnestor = 0.d0
348  ENDIF
349  xmve = xmve_tel
350  grav = grav_tel
351 !
352 ! DETERMINES THE TOTAL NUMBER OF TIMESTEPS : NIT
353 ! INITIAL TIME : AT0
354 ! TIMESTEP : DT
355  at0 = t_tel
356  dt = dt_tel
357  nit = telnit
358 ! VALIDATES AGAINST THE LAST TIMESTEP
359 !PTRK -- VALNIT has to be checked if the computation is efficient! -> CHECK THIS PART!!
360  valnit = nit
361 ! VALNIT WILL BE USED FOR CALLING VALIDA
362 ! JMH ON 21/06/2016: IT SEEMS THAT THE FORMULA ON NEXT LINE IS OUTDATED...
363 ! VALNIT = (TELNIT/PERICOU)*PERICOU-PERICOU+1
364  valnit = telnit
365 ! MODIFICATION JMH + CV: TO AVOID 2 SUCCESSIVE CALLS TO VALIDA
366 ! WHEN BEDLOAD AND SUSPENSION
367  IF(grafcount.GT.telnit) valnit=nit+1
368  leopr = grafcount
369  lispr = listcount
370 ! CHARR, SUSP AND TIME STEP MONITORED BY CALLING PROGRAM
371 ! AT0=T_TEL
372  dt=mofac*dt_tel
373 !
374 ! RESUMES GAIA COMPUTATION
375 !
376  IF(gai_files(gaipre)%NAME(1:1).NE.' ') THEN
377 !
378 ! READS THE HYDRO AND SEDIMENTOLOGICAL VARIABLES
379 !
380  IF(debug.GT.0) WRITE(lu,*) 'READ_DATASET'
382  & varsor,npoin,numenx,at0,textpr,trouve,alire,
383  & .true.,.true.,maxvar)
384  IF(debug.GT.0) WRITE(lu,*) 'END_READ_DATASET'
385 !
386  IF(debug.GT.0) WRITE(lu,*) 'RESCUE_GAIA'
387  CALL rescue_gaia(hn%R,z%R,zf%R,zr%R,es,hw%R,tw%R,
388  & thetaw%R,npoin,nomblay,trouve,alire,pass,icf,.true.,maxvar)
389  IF(debug.GT.0) WRITE(lu,*) 'SORTIE DE RESCUE_GAIA'
390 !
391  ENDIF
392 !
393 ! READS THE LAST RECORD : WAVE FILE
394 !
395 ! NOTE : GAI_FILES(GAICOU)%NAME SET TO ' ' IF HOULE=NO
396 !
397  IF (gai_files(gaicou)%NAME(1:1).NE.' ') THEN
398 !
399  IF(debug.GT.0) WRITE(lu,*) 'SUITE_HOULE'
400  WRITE(lu,*) ' LECTURE HOULE :',gai_files(gaicou)%NAME
402  & varsor,npoin,numenx,at,textpr,trouve,alirh,
403  & .true.,.true.,maxvar)
404  IF(debug.GT.0) WRITE(lu,*) 'END_SUITE_HOULE'
405 ! TRACES IF WAVE DATA HAVE BEEN FOUND
406  IF(trouve(11).EQ.1) hw%TYPR='Q'
407  IF(trouve(12).EQ.1) tw%TYPR='Q'
408  IF(trouve(13).EQ.1) thetaw%TYPR='Q'
409 !
410  ENDIF
411 !
412  at0=t_tel
413 !
414 ! INFORMATION ON SUSPENSION SENT BACK
415  charr_tel = charr
416  susp_tel = susp
417 !
418 ! OV INSTEAD OF OS IN ORDER TO AVOID PROBLEMS WITH QUASI-BUBBLE ELEMENTS
419 ! OPERATES ONLY ON THE (1:NPOIN) RANGE OF THE TELEMAC FIELDS
420 ! IT IS A HIDDEN DISCRETISATION CHANGE
421 !
422  CALL ov('X=Y ', x=u2d%R, y=u_tel%R, dim1=npoin)
423  CALL ov('X=Y ', x=v2d%R, y=v_tel%R, dim1=npoin)
424  CALL ov('X=Y ', x=hn%R, y=h_tel%R, dim1=npoin)
425 ! BOTTOM GIVEN BY CALLING PROGRAMME
426  CALL os('X=Y ', x=zf, y=zf_tel)
427 !
428 ! CLIPS NEGATIVE DEPTHS
429 !
430  IF(optban.GT.0) THEN
431  DO i = 1,npoin
432  IF(hn%R(i).LT.hmin) THEN
433  u2d%R(i)=0.d0
434  v2d%R(i)=0.d0
435  hn%R(i) = hmin
436  ENDIF
437  ENDDO
438  ENDIF
439 !
440 ! CASE OF TRIPLE COUPLING
441 !
442  IF(inclus(coupling,'TOMAWAC')) THEN
443 ! incident wave direction
444  CALL os( 'X=Y ',thetaw,thetaw_tel)
445 ! Wave period
446  CALL os( 'X=Y ', tw, tw_tel)
447 ! significant wave height
448  CALL os( 'X=Y ', hw , hw_tel)
449 ! Bottom orbital velocity
450  CALL os( 'X=Y ', uw , uw_tel)
451  hw%TYPR='Q'
452  tw%TYPR='Q'
453  thetaw%TYPR='Q'
454  uw%TYPR='Q'
455  ENDIF
456 !
457 ! END COUPLING
458 !
459  IF(debug.GT.0) WRITE(lu,*) 'USER_FORCING_GAIA'
460  IF(.NOT.debu) THEN
461  CALL user_forcing_gaia
462  ENDIF
463  IF(debug.GT.0) WRITE(lu,*) 'END_USER_FORCING_GAIA'
464 !
465 ! AT THIS LEVEL U2D,V2D,H AND ZF MUST HAVE BEEN DEFINED
466 ! EITHER BY READ_DATASET, USER_FORCING_GAIA OR CALLING PROGRAM
467 !
468 ! NOW COMPUTES FUNCTIONS OF U2D,V2D,HN AND ZF
469 !
470 ! FREE SURFACE
471  IF(debug.GT.0) WRITE(lu,*) 'FREE SURFACE'
472  CALL os('X=Y+Z ', x=z, y=zf, z=hn)
473  IF(debug.GT.0) WRITE(lu,*) 'END FREE SURFACE'
474 !
475 ! CHECKS THE WAVE DATA
476 !
477  IF(houle) THEN
478  IF(uw%TYPR.EQ.'Q') THEN
479 ! IF(HW%TYPR .NE.'Q'.OR.
480  ELSEIF(hw%TYPR .NE.'Q'.OR.
481  & tw%TYPR .NE.'Q'.OR.
482  & thetaw%TYPR.NE.'Q') THEN
483  WRITE(lu,*) ' '
484  WRITE(lu,*) ' '
485  WRITE(lu,*) 'MISSING WAVE DATA'
486  IF(hw%TYPR.NE.'Q') WRITE(lu,*) 'WAVE HEIGHT HM0'
487  IF(tw%TYPR.NE.'Q') WRITE(lu,*) 'PEAK PERIOD TPR5'
488  IF(thetaw%TYPR.NE.'Q') WRITE(lu,*) 'MEAN DIRECTION'
489  CALL plante(1)
490  stop
491  ENDIF
492  ENDIF
493 !
494 ! END OF HYDRODYNAMIC INITIALISATION
495 !
496 ! COMPUTES AREAS (WITHOUT MASKING)
497 !
498  IF(debug.GT.0) WRITE(lu,*) 'VECTOR FOR VOLU2D'
499  CALL vector(volu2d,'=','MASBAS ',
500  & ielmh_gai,1.d0,
501  & t1,t1,t1,t1,t1,t1,mesh,.false.,maskel)
502  IF(debug.GT.0) WRITE(lu,*) 'END VECTOR FOR VOLU2D'
503 ! V2DPAR : LIKE VOLU2D BUT IN PARALLEL VALUES COMPLETED AT
504 ! INTERFACES BETWEEN SUBDOMAINS
505  CALL os('X=Y ',x=v2dpar,y=volu2d)
506  IF(ncsize.GT.1) CALL parcom(v2dpar,2,mesh)
507 ! INVERSE OF VOLUMES (DONE WITHOUT MASKING)
508  CALL os('X=1/Y ',x=unsv2d,y=v2dpar,
509  & iopt=2,infini=0.d0,zero=1.d-12)
510 !
511 ! INITIALIZATION FOR SETTLING VELOCITIES
512 !
513  CALL settling_vel
514 !
515 !
516 ! INITIALIZATION FOR SHIELDS PARAMETER
517 !
518  CALL shields
519 !
520 ! INITIALISATION FOR SEDIMENT IN BED MODEL
521 !
522  IF (bed_model.EQ.1.OR.bed_model.EQ.2)THEN
523 ! model multilayer or multilayer with consolidation
524  IF(debug.GT.0) WRITE(lu,*) 'BED1_INIT_SEDIMENT_GAIA'
527  & maxvar)
528  IF(debug.GT.0) WRITE(lu,*) 'END BED1_INIT_SEDIMENT_GAIA'
529 
530 
531 ! CONVERSION MASS TO THICKNESS
532 
533  IF(nestor.OR.vsmtype==1) THEN
534 
535 ! | faktor = 1 / ( density * area * ( 1 - porosity ) ) |
536 ! | => mass * faktor = thickness |
537 ! | but VARIOUS MASSES ARE STILL IN [kg/m**2] ==>
538 !
539  DO icla = 1,nsicla
540  DO i = 1,npoin
541  m2t(icla,i) = 1.d0/(xmvs0(icla)*volu2d%R(i)*(1.d0-xkv0(1)))
542  ENDDO
543  mpa2t(icla) = 1.d0/( xmvs0(icla) * (1.d0-xkv0(1)) )
544  ENDDO
545  ENDIF
546 !
547  IF(hirano) THEN
548  IF(debug.GT.0) WRITE(lu,*)'BED1_UPDATE 1'
549  CALL bed1_update(zr,zf,volu2d)
550  IF(debug.GT.0) WRITE(lu,*)'END BED1_UPDATE 1'
551  IF(debug.GT.0) WRITE(lu,*)'BED1_UPDATE_ACTIVELAYER_HIRANO'
553  IF(debug.GT.0) THEN
554  WRITE(lu,*)'END BED1_UPDATE_ACTIVELAYER_HIRANO'
555  ENDIF
556  IF(debug.GT.0) WRITE(lu,*)'BED1_UPDATE 2'
557  CALL bed1_update(zr,zf,volu2d)
558  IF(debug.GT.0) WRITE(lu,*)'END BED1_UPDATE 2'
559 !
560 ! INITIALIZATION OF CVSM MODEL
561  IF(vsmtype.EQ.1) THEN
562  IF(.NOT.hirano) THEN
563  WRITE(lu,*)'WITH CVSM MODEL HIRANO MUST BE CHOSEN!'
564  CALL plante(1)
565  ENDIF
566  IF(nmud.GT.0) THEN
567  WRITE(lu,*)'CVSM MODEL ONLY WITH SAND/GRAVEL CLASSES'
568  CALL plante(1)
569  ENDIF
570 
571 ! ---- PRINTOUTS TO C-VSP
572 !
573  IF(cvsmpperiod.EQ.0) cvsmpperiod = leopr
574  cvsm_out = .false.
575  IF(cvsmpperiod*(lt/cvsmpperiod).EQ.lt) cvsm_out = .true.
576  WRITE(lu,*) ' '
577  WRITE(lu,*) '------------------------------------------------'
578  WRITE(lu,*) 'C-VSM MODEL'
579  WRITE(lu,*) 'CONTINUOUS VERTICAL GRAIN SORTING STRATIGRAPHY'
580  WRITE(lu,*) ' '
581  WRITE(lu,*) 'ACTIVE LAYER THICKNESS MODEL:', alt_model
582 !
583  IF (pro_max_max .GT. 250) THEN
584  WRITE(lu,*) 'HIGH NUMBER OF SECTIONS IS EXPENSIVE',
585  &pro_max_max
586  WRITE(lu,*) 'BETTER < 250 and > 4 + 4 x NUMBER OF FRACTIONS'
587  ENDIF
588 !
589  IF ((cvsmpperiod / nit) .GT. 5) THEN
590  WRITE(lu,*) 'HIGH NUMBER OF FULL CVSM PRINTOUTS'
591  WRITE(lu,*) 'ATTENTION: DISK SPACE AND SIMULATION TIME'
592  WRITE(lu,*) 'ADAPT C-VSM FULL PRINTOUT PERIOD'
593  ENDIF
594 !
595  WRITE(lu,*) ' '
596  CALL cvsp_init_gaia()
597  WRITE(lu,*)'CVSM INITIALISED!'
598  WRITE(lu,*) '------------------------------------------------'
599  ENDIF
600  ENDIF ! HIRANO.AND.VSMTYPE=1
601  ELSE ! BED_MODEL NE 1 OR 2
602  WRITE(lu,*)'ONLY BED_MODEL 1 OR 2 FOR TIS TIME'
603  stop
604  ENDIF
605 
606 ! INITIALIZATION OF ACLADM AND UNLADM
607  IF(nsand.GT.0) THEN
609  ENDif!
610 !
611 ! CALUCLATION INITIAL MASS IN ACTIVE LAYER
612 ! AT THE MOMENT MASS0ACT IS ONLY COMPUTED FOR MASS
613 ! BALANCE
614 !
615  IF(bilma) THEN
616  mass0act = 0.d0
617  IF(nsand.GT.0) THEN
618  DO ipoin=1,npoin
619  DO isand=1,nsand
620  k = num_isand_icla(isand)
621  mass0act(k) = mass0act(k) + mass_sand(isand,1,ipoin)
622  & *volu2d%R(ipoin)
623  ENDDO
624  ENDDO
625  ENDIF
626  IF(nmud.GT.0) THEN
627  DO ipoin=1,npoin
628  DO imud=1,nmud
629  k = num_imud_icla(imud)
630  mass0act(k) = mass0act(k) + mass_mud(imud,1,ipoin)
631  & *volu2d%R(ipoin)
632  ENDDO
633  ENDDO
634  ENDIF
635  IF(ncsize.GT.1) THEN
636  DO icla=1,nsicla
637  mass0act(icla) = p_sum(mass0act(icla))
638  ENDDO
639  ENDIF
640  ENDIF
641 !
642 ! MEAN VELOCITY
643 !
644  CALL os('X=N(Y,Z)',x=unorm,y=u2d,z=v2d)
645 !
646 ! DIRECTION OF CURRENT (° and trigo)
647 ! ATAN2(Y,Z)
648  CALL os('X=A(Y,Z)',x=thetac,y=v2d,z=u2d)
649 ! radians to degres
650  DO i=1,npoin
651  thetac%R(i)=thetac%R(i)*180/pi
652  thetac%R(i)=modulo(thetac%R(i),360.d0)
653  ENDDO
654 !
655 ! WAVE ORBITAL VELOCITY
656 !
657  IF(houle) THEN
658 ! directly if found in hydro file; otherwise compute with CALCUW_GAIA
659  IF(uw%TYPR.NE.'Q') THEN
660  CALL calcuw_gaia(uw%R,hn%R,hw%R,tw%R,grav,npoin, type_houle)
661  ELSE
662  IF(type_houle==1)THEN
663 ! UW is OK!
664  CONTINUE
665  ELSEIF(type_houle==2)THEN
666 ! Case of tomawac (irregular waves)
667 ! UW calculated from sprectum is Urms
668 ! transformation for a Jonswap spectrum (Soulsby 1993)
669  DO i=1,npoin
670  uw%R(i)=sqrt(2.d0)*uw%R(i)
671  ENDDO
672  ELSE
673  WRITE(lu,*)'VALUE OF TYPE OF WAVES IS NOT OK'
674  CALL plante(1)
675  stop
676  ENDIF
677  ENDIF
678  ENDIF
679 !
680 ! TOTAL STRESS AT THE BOTTOM
681 !
682  IF(debug.GT.0) WRITE(lu,*) 'TOB_GAIA'
683  CALL tob_gaia(tob,tobw,tobcw_mean,tobcw_max,thetac, thetaw,mu,
684  & ks,ksp,ksr,cf,fw,
685  & uetcar,cf_tel,ks_tel,code,icr,kspratio,houle,
686  & grav,xmve,xmvs0(1),vce,karman,zero,
687  & hmin,hn,acladm,unorm,uw,tw,npoin,kscalc,iks,
688  & deltar, h_tel)
689  IF(debug.GT.0) WRITE(lu,*) 'END TOB_GAIA'
690 !
691 ! INITIALISATION FOR TRANSPORT
692 !
693  IF(debug.GT.0) WRITE(lu,*) 'INIT_TRANSPORT_GAIA'
694  CALL init_transport_gaia(hiding,nsicla,npoin,
695  & t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t14,
696  & charr,qs_c,qsxc,qsyc,calfa_cl,salfa_cl,coefpn,slopeff,
697  & susp,qs,qscl,qscl_c,qscl_s,
698  & unorm,u2d,v2d,hn,cf,mu,tob,tobw,uw,tw,thetaw,fw,houle,
699  & acladm,unladm,ksp,ksr,
700  & icf,hidfac,xmvs0,xmve,grav,vce,hmin,karman,
701  & zero,pi,ac,cstaeq,seccurrent,bijk,ielmt,mesh,dcla,xwc,sedco,
702  & u3d,v3d,code,h_tel,
703  & hw,thetac,tobcw_mean,tobcw_max)
704  IF(debug.GT.0) WRITE(lu,*) 'END INIT_TRANSPORT_GAIA'
705 !
706 ! PRINT OF INITIAL CONDITION
707 !
708  CALL entete_gaia(1,at0,0)
709 ! PREPARES RESULTS
710  CALL predes_gaia(0,at0,.true.,code,listcount)
711 !
712 ! PRINTS OUT THE RESULTS
713 !
714  IF(debug.GT.0) WRITE(lu,*) 'BIEF_DESIMP'
715  CALL bief_desimp(gai_files(gaires)%FMT,varsor,
716  & npoin,gai_files(gaires)%LU,
717  & at0,0,listcount,grafcount,sorleo,sorimp,maxvar,
718  & texte,ptinig,ptinil,
719  & ileo=yagout)
720  IF(debug.GT.0) WRITE(lu,*) 'END BIEF_DESIMP'
721 !
722 
723 ! COUPLING
724 !
725 ! INITIALISING NESTOR
726  IF(nestor) THEN
727  CALL nestor_interface_gaia(1,grafcount,xmvs0,xkv0(1),volu2d)
728  IF(lng.EQ.1) WRITE(lu,*) 'GAIA COUPLE AVEC NESTOR'
729  IF(lng.EQ.2) WRITE(lu,*) 'GAIA COUPLED WITH NESTOR'
730  ENDIF
731 !
732  WRITE(lu,*) 'GAIA COUPLED WITH: ',code
733 !
734 !------------------------------------------------------------------
735 !
736 ! ADJOINTWARE: ALGORITHMIC DIFFERENTIATION
737 #if defined COMPAD
739 #endif
740 ! ADJOINTWARE
741 !------------------------------------------------------------------
742  RETURN
743  END
integer, pointer nptfr
Number of boundary points.
type(bief_obj), target uw
Orbital wave velocity.
double precision, dimension(nsiclm) summcumucl
Cumulated erosion / deposition mass per class.
logical bilma
Mass balance.
type(bief_obj), target liqbor
Type of boundary conditions on sand transport rate.
subroutine write_mesh(FFORMAT, NFILE, MESH, NPLAN, DATE, TIME, T1, T2, PARALL, NPTIR, NGEO, GEOFORMAT, LATLONG)
Definition: write_mesh.f:8
double precision, dimension(nsiclm) sum_deposition
Cumulated sediment mass deposition.
subroutine ad_gaia_begin
Definition: ad_gaia_begin.f:3
subroutine ov(OP, X, Y, Z, C, DIM1)
Definition: ov.f:7
type(bief_obj), target maskel
Mask.
integer, target gaires
Various files ranks, which are also logical units if no coupling.
double precision hmin
Minimal value of the water height: below this value, the sediment flow rate is set to 0...
type(bief_obj), target u2d
Components of depth-averaged velocity.
double precision, dimension(nsiclm) sumrmascl
Cumulated evolution mass per class.
double precision pi
Pi.
subroutine nestor_interface_gaia(OPTION, GRAFCOUNT, XMVS0, XKV01, VOLU2D)
character(len=72) titca
Title of the case.
integer produc
Matrix-vector product.
subroutine bed1_init_sediment_gaia(NSICLA, ELAY, ZF, ZR, NPOIN, XMVS0, ES, NOMBLAY, DEBU, VOLU2D, NUMSTRAT, MAXVAR)
integer lispr
LISTING PRINTOUT PERIOD.
type(bief_obj), target v2d
logical, target susp
Suspension : yes if there is at least one suspended sediment this is the case if there is mud or if s...
subroutine settling_vel
Definition: settling_vel.f:4
integer, target nomblay
Number of bed load model layers = NUMSTRAT+1 to take the active layer into account.
type(bief_obj), target it3
double precision, dimension(:), pointer x
2d coordinates of the mesh
type(bief_obj), target zr
Non erodable (rigid) bottom elevation.
type(bief_obj), target volu2d
Integral of bases.
double precision, target dt
Time step It may be different from the one in TELEMAC because of the morphological factor...
double precision, dimension(:,:), allocatable, target es
Layer thicknesses as double precision.
type(bief_file), dimension(maxlu_gai), target gai_files
For storing information on files.
logical houle
Include wave effects.
subroutine tob_gaia(TOB, TOBW, TOBCW_MEAN, TOBCW_MAX, THETAC, THETAW, MU, KS, KSP, KSR, CF, FW, UETCAR, CF_TEL, KS_TEL, CODE, ICR, KSPRATIO, HOULE, GRAV, XMVE, XMVS, VCE, KARMAN, ZERO, HMIN, HN, ACLADM, UNORM, UW, TW, NPOIN, KSCALC, IKS, DELTAR, H_TEL)
Definition: tob_gaia.f:11
integer numstrat
Number of layers of initial stratification.
double precision mofac
Morphological factor on the hydrodynamics: distorts the evolution of the hydrodynamics with respect t...
double precision xmve
Water density (from steering file of T2D or T3D)
type(bief_obj), target zrl
Reference level for Nestor.
integer bed_model
Bed model (3 choices, cf. dico)
integer, target gaipre
type(bief_obj), target afbor
Flux condition nu df/dn=afbor * f + bfbor.
integer, dimension(3) martim
Original hour of time.
subroutine init_transport_gaia(HIDING, NSICLA, NPOIN, T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11, T12, T14, CHARR, QS_C, QSXC, QSYC, CALFA_CL, SALFA_CL, COEFPN, SLOPEFF, SUSP, QS, QSCL, QSCL_C, QSCL_S, UNORM, U2D, V2D, HN, CF, MU, TOB, TOBW, UW, TW, THETAW, FW, HOULE, ACLADM, UNLADM, KSP, KSR, ICF, HIDFAC, XMVS0, XMVE, GRAV, VCE, HMIN, KARMAN, ZERO, PI, AC, CSTAEQ, SECCURRENT, BIJK, IELMT, MESH, DCLA, XWC, SEDCO, U3D, V3D, CODE, H_TEL, HW, THETAC, TOBCW_MEAN, TOBCW_MAX)
integer optban
Option for the treatment of tidal flats.
logical hirano
Hirano model used.
double precision, dimension(:,:), allocatable sumbedload_b
Cumulated bedload on boundary for every class (kg): variable for mass balance.
subroutine ad_gaia_initialisation_begin
character(len=32), dimension(maxvar) texte
Name of output variable.
logical, dimension(maxvar) sorleo
Graphical output.
double precision, dimension(:), pointer y
type(bief_obj), target it4
type(bief_obj), target thetaw
Wave direction (deg wrt ox axis) !!!!!some say oy axis!!!!!
type(bief_obj), target numliq
Liquid boundary numbering.
subroutine mean_grain_size_gaia
integer, parameter maxvar
Maximum number of output variables.
type(bief_obj), target it2
integer optass
Matrix storage.
double precision, dimension(nsiclm) xmvs0
Sand density.
integer, parameter kent
double precision, dimension(nsiclm) massnestor
Sediment mass from Nestor per class per time step.
type(bief_obj), target elay
Active layer thickness.
double precision, dimension(nsiclm) sum_erosion
Cumulated sediment mass erosion.
integer nsand
Total number of sand.
double precision zero
Parameter used for clipping variables or testing values against zero.
integer, target gaicou
type(bief_obj), target hw
Significant wave height.
integer icf
Bed-load transport formula.
integer alt_model
CHOOSE A MODEL FOR ESTIMATION OF A DYNAMIC ACTIVE LAYER THICKNESS.
double precision, dimension(nsiclm) summassnestor
Cumulated sediment mass from Nestor per class.
type(bief_obj), target radsec
Curve radius for secondary currents.
integer, pointer nelem
Number of elements in the mesh.
subroutine predes_gaia(LLT, AAT, YAGOUT, CODE, LISTCOUNT)
Definition: predes_gaia.f:7
integer, dimension(maxvar) alire
Variables to read if computation is continued 0 : DISCARD 1 : READ (SEE SUBROUTINE NOMVAR) ...
type(bief_obj), target bfbor
subroutine entete_gaia(IETAPE, AT, LT)
Definition: entete_gaia.f:7
type(bief_obj), target hn
Water depth.
subroutine ad_gaia_initialisation_end
subroutine shields
Definition: shields.f:4
subroutine cvsp_init_gaia
Definition: cvsp_init_gaia.f:4
type(bief_obj), target unsv2d
Inverse of integral of bases.
subroutine write_header(FFORMAT, NRES, TITLE, NVAR, NOMVAR, OUTVAR)
Definition: write_header.f:7
subroutine gaia_init(GRAFCOUNT, LISTCOUNT, TELNIT, U_TEL, V_TEL, H_TEL, ZF_TEL, UETCAR, DELTAR, CF_TEL, KS_TEL, CODE, U3D, V3D, T_TEL, DT_TEL, CHARR_TEL, SUSP_TEL, THETAW_TEL, HW_TEL, TW_TEL, UW_TEL, YAGOUT, XMVE_TEL, GRAV_TEL)
Definition: gaia_init.F:10
type(bief_obj), target q2bor
Imposed solid transport at the boundary In m2/s, total, read in the boundary conditions file...
logical function inclus(C1, C2)
Definition: inclus.f:7
integer, target nsicla
Number of sediment classes of bed material (less than NISCLM)
logical spheri
Work in spherical coordinates (hard-coded)
integer vsmtype
For the Continous Vertical Sorting MODEL.
logical nestor
Coupling with NESTOR.
double precision grav
Gravity acceleration.
type(bief_obj), pointer t2
double precision, dimension(:,:), allocatable, target m2t
conversion mass to thickness
integer, dimension(maxvar) alirh
type(bief_obj), target z
Free surface elevation.
integer cvsmpperiod
Printout Period for Full Vertical Sorting Model: PRO_D & PRO_F.
double precision, dimension(nlaymax), target xkv0
Initial porosity by layers.
integer lvmac
Vector length (for vectorisation) from steering file.
subroutine init_constant_gaia(KARIM_HOLLY_YANG, KARMAN, PI)
subroutine front2(NFRLIQ, LIHBOR, LIUBOR, X, Y, NBOR, KP1BOR, DEJAVU, NPOIN, NPTFR, KLOG, LISTIN, NUMLIQ, MAXFRO)
Definition: front2.f:8
logical msk
Include masking.
type(bief_obj), pointer t1
Aliases for work vectors in tb.
integer nmud
Total number of muds.
type(bief_obj), target mask
Block of masks.
subroutine user_forcing_gaia
subroutine vector(VEC, OP, FORMUL, IELM1, XMUL, F, G, H, U, V, W, MESH, MSK, MASKEL, LEGO, ASSPAR)
Definition: vector.f:7
logical debu
Computation continued.
integer, target gaigeo
integer, target lt
Numero du pas de temps.
integer maxfro
Maximum number of (liquid boundaries, solid boundaries)
subroutine init_zero_gaia
Definition: init_zero_gaia.f:4
subroutine find_variable(FFORMAT, FID, VAR_NAME, RES, N, IERR, TIME, EPS_TIME, RECORD, TIME_RECORD, OFFSET)
Definition: find_variable.f:8
double precision karim_holly_yang
Karim, Holly & Yang constant for hiding (with multiclass beds)
type(bief_obj), target tw
Mean wave period.
integer debug
Debugger.
integer pro_max_max
Maximum Number of Profile SECTIONS.
logical cvsm_out
C-VSM WRITES OUT (OR NOT) IN THIS TIMESTEP.
type(bief_obj), target ebor
Imposed bed evolution at the boundary [m].
type(bief_obj), target varsor
Block of variables for output.
subroutine os(OP, X, Y, Z, C, IOPT, INFINI, ZERO)
Definition: os.f:7
type(bief_obj), target zf
Bottom elevation.
type(bief_mesh), target mesh
Mesh structure.
subroutine parcom(X, ICOM, MESH)
Definition: parcom.f:7
subroutine bed1_update_activelayer_hirano
double precision, dimension(:), allocatable, target mpa2t
conversion mass per area to thickness
character(len=20) equa
Equation solved.
integer ielmt
Missing comment.
logical, target charr
Include bedload in the simulation.
double precision at0
integer leopr
GRAPHIC PRINTOUT PERIOD.
type(bief_obj), target liebor
Type of boundary conditions on bed evolution.
subroutine calcuw_gaia(UW, H, HW, TW, GRAV, NPOIN, TYPE_HOULE)
Definition: calcuw_gaia.f:7
logical second_susp_step
Logical used for coupling with T2D/T3D when suspension activated.
subroutine read_dataset(FFORMAT, FID, VARSOR, NPOIN, RECORD, AT, VAR_LIST, TROUVE, ALIRE, LISTIN, LASTRECORD, MAXVAR)
Definition: read_dataset.f:8
subroutine inbief(LIHBOR, KLOG, IT1, IT2, IT3, LVMAC, IELMX, LAMBD0, SPHERI, MESH, T1, T2, OPTASS, PRODUC, EQUA, MESH2D)
Definition: inbief.f:8
integer, pointer nelmax
Maximum number of elements in the mesh.
subroutine leclim(LIHBOR, LIUBOR, LIVBOR, LITBOR, HBOR, UBOR, VBOR, TBOR, CHBORD, ATBOR, BTBOR, NPTFR, CODE, TRAC, FFORMAT, NGEO, KENT, KENTU, KSORT, KADH, KLOG, KINC, NUMLIQ, MESH, BOUNDARY_COLOUR, NPTFR2)
Definition: leclim.f:10
subroutine rescue_gaia(H, S, ZF, ZR, ES, HW, TW, THETAW, NPOIN, NOMBLAY, TROUVE, ALIRE, PASS, ICF, LISTI, MAXVAR)
Definition: rescue_gaia.f:8
subroutine bed1_update(ZR, ZF, VOLU2D)
Definition: bed1_update.f:7
double precision karman
Karman constant.
integer nfrliq
Number of liquid boundaries.
integer hind_type
Hindered settling method linked to floculation effects, 1:WHITEHOUSE ET AL. (2000), 2:WINTERWERP (1999), at the moment hard-set to 1 in gaia.F.
integer, parameter klog
character(len=path_len), target coupling
type(bief_obj), target it1
Integer working arrays.
subroutine bief_desimp(FORMAT_RES, VARSOR, N, NRES, AT, LT, LISPRD, LEOPRD, SORLEO, SORIMP, MAXVAR, TEXTE, PTINIG, PTINIL, MESH, IIMP, ILEO, COMPGRAPH)
Definition: bief_desimp.f:9
integer nadvar
Number of differentiating arrays, and those with a given name.
type(bief_obj), target lihbor
Type of boundary conditions for h.
integer nvar_ratios
Index in varsor for output variables.
integer, pointer npoin
Number of 2d points in the mesh.
type(bief_obj), target v2dpar
Integral of bases in parallel.
character(len=32), dimension(maxvar) textpr
Name of variable in previous computation file.
type(bief_obj), target boundary_colour
Last column of the boundary condition file.
integer, dimension(3) mardat
Orginal date of time.
Definition: bief.f:3