The TELEMAC-MASCARET system  trunk
t3d_aed2.F90
Go to the documentation of this file.
1 ! ***************
2  MODULE t3d_aed2
3 ! ***************
4 !
5 !***********************************************************************
6 ! WAQTEL V8P1
7 !***********************************************************************
8 !
9 !brief Interface for TELEMAC3D to AED2 modules (libaed2)
10 !+ This is the main interface module that manages the connection
11 !+ with the host hydrodynamic model;
12 !+ done through the 4 PUBLIC functions listed below.
13 !
14 !history
15 !+ 08/2013
16 !+ Interface for FV (Finite Volume) Model to AED2 modules (libaed2)
17 !+ Designed for TUFLOW-FV, released by BMT-WBM:
18 !+ http://www.tuflow.com/Tuflow%20FV.aspx
19 !+ Originally fv_aed2.F90
20 !+
21 !+ Developed by:
22 !+ AquaticEcoDynamics (AED) Group
23 !+ School of Earth & Environment
24 !+ (C) The University of Western Australia
25 !+ Copyright by the AED-team @ UWA under the GNU Public License
26 !+ www.gnu.org
27 !
28 !history
29 !+ 03/2016
30 !+ Add new env variables and feedback links
31 !
32 !history M. Jodeau (EDF)
33 !+ 04/2016
34 !+ Modifications from fv_aed2.F90 and adaptation to TELEMAC3D
35 !
36 !history C.-T. PHAM (LNHE)
37 !+ 15/11/2018
38 !+ V8P0
39 !+ Variation of salinity, temperature, density with the order expected
40 !+ by AED2 (from top to bottom, then 2D number) + real wind velocity
41 !
42 #if defined HAVE_AED2
43 #include "aed2.h"
44 #endif
45 
46 !#define FV_AED_VERS "0.9.21"
47 
48 !#ifndef DEBUG
49 !#define DEBUG 0
50 !#endif
51 !#define _NO_ODE_ 1
52 #if defined HAVE_AED2
53  USE aed2_common
54  USE fv_zones
57  USE bief
58 
59  IMPLICIT NONE
60 
61  PUBLIC init_aed2_models, &
67 
68 !#--------------------------------------------------------------------------#
69 !# MODULE DATA
70 
71  aed_real :: kw, ksed
72  INTEGER :: solution_method
73 
74 !# MAIN ARRAYS STORING/POINTING TO THE STATE AND DIAGNOSTIC VARIABLES
75  aed_real,DIMENSION(:,:),POINTER :: cc, cc_diag
76 
77 !# ARRAY POINTING TO THE LAGRANGIAN PARTICLE MASSES AND DIAGNOSTIC PROPERTIES
78  aed_real,DIMENSION(:,:),POINTER :: pp
79 
80 !# NAME OF FILE BEING USED TO LOAD INITIAL VALUES FOR BENTHIC OR BENTHIC_DIAG VARS
81  CHARACTER(LEN=128) :: init_values_file = ''
82 
83 !# MAPS OF SURFACE, BOTTOM AND WET/DRY (ACTIVE) CELLS
84  LOGICAL,DIMENSION(:),POINTER :: active
85  INTEGER,DIMENSION(:),POINTER :: surf_map, benth_map
86 
87 !# ARRAYS FOR WORK, VERTICAL MOVEMENT, AND CROSS-BOUNDARY FLUXES
88  aed_real,ALLOCATABLE,DIMENSION(:,:) :: flux
89  aed_real,ALLOCATABLE,DIMENSION(:) :: ws, total
90  aed_real,ALLOCATABLE,DIMENSION(:) :: min_, max_
91 
92 !# ARRAYS FOR ENVIRONMENTAL VARIABLES (USED IF THEY ARE NOT SUPPLIED EXTERNALLY)
93  aed_real,ALLOCATABLE,DIMENSION(:),TARGET :: nir
94  aed_real,ALLOCATABLE,DIMENSION(:),TARGET :: par
95  aed_real,ALLOCATABLE,DIMENSION(:),TARGET :: uva
96  aed_real,ALLOCATABLE,DIMENSION(:),TARGET :: uvb
97  aed_real,DIMENSION(:),POINTER :: lpar
98 
99  aed_real,TARGET :: col_taub ! A TEMP VAR FOR THE TAUB FOR COLUMN (COMPUTED FROM USTAR_BED)
100 
101 !# EXTERNAL VARIABLES
102  aed_real :: dtaed2 ! MODIF MAG !DTAED2
103  aed_real :: eps_aed2=1.d-9 ! MODIF MAG !DTAED2
104  aed_real,DIMENSION(:), POINTER :: temp_aed2, salt, rho_aed2
105  aed_real,DIMENSION(:), POINTER :: h_aed2, z_aed2
106  aed_real,DIMENSION(:), POINTER :: extcoeff, tss, bio_drag
107  aed_real,DIMENSION(:), POINTER :: i_0, wnd, air_temp, rainaed2 ! MODIF MAG !RAINAED2
108  aed_real,DIMENSION(:), POINTER :: area, bathy, shadefrac
109  aed_real,DIMENSION(:), POINTER :: rainloss
110  aed_real,DIMENSION(:), POINTER :: fsed_setl, ustar_bed
111  INTEGER, DIMENSION(:,:),POINTER :: mat
112 
113 !# SWITCHES FOR CONFIGURING MODEL OPERATION AND ACTIVE LINKS WITH THE HOST MODEL
119 
120  LOGICAL :: old_zones = .true.
121  LOGICAL :: do_zone_averaging = .false.
122 
123 !# INTEGERS STORING NUMBER OF VARIABLES BEING SIMULATED
125  INTEGER :: n_vars_diag_sheet
126 
127  CONTAINS
128 
129 ! ***************************
130  SUBROUTINE init_aed2_models &
131 ! ***************************
132 !
133  (namlst,dname,nwq_var,nben_var,ndiag_var)
134 !
135 !***********************************************************************
136 ! WAQTEL V7P3
137 !***********************************************************************
138 !
139 !brief This subroutine is called by TELEMAC-3D to define numbers of
140 !+ variables. TELEMAC-3D will allocate the variables after return
141 !+ from this subroutine
142 !
143 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
144 !| NAMLST |-->| UNIT FOR AED2 STEERING FILE
145 !| NBEN_VAR |<--| NUMBER OF BENTHIC VARIABLES
146 !| NDIAG_VAR |-->| NUMBER OF DIAGNOSIS VARIABLES
147 !| NWQ_VAR |<--| NUMBER OF WATER QUALITY VARIABLES
148 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
149 !
150 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
151 !
152  INTEGER, INTENT(IN) :: NAMLST
153  INTEGER, INTENT(OUT) :: NWQ_VAR,NBEN_VAR,NDIAG_VAR
154  CHARACTER(LEN=*), INTENT(IN) :: DNAME
155 !
156 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
157 !
158  CHARACTER(LEN=30), ALLOCATABLE :: NAMES(:)
159  CHARACTER(LEN=30), ALLOCATABLE :: BENNAMES(:)
160  CHARACTER(LEN=30), ALLOCATABLE :: DIAGNAMES(:)
161  CHARACTER(LEN=128) :: AED2_NML_FILE
162  CHARACTER(LEN=200) :: TNAME
163  aed_real :: base_par_extinction, tss_par_extinction
164  INTEGER :: STATUS
165  INTEGER :: N_SD, I, J
166  TYPE(aed2_variable_t),POINTER :: TVAR
167 
168  CHARACTER(LEN=64) :: MODELS(64)
169  namelist /aed2_models/ models
170 
171  namelist /aed2_bio/ solution_method, link_bottom_drag, &
173  link_water_clarity, aed2_nml_file, &
174  link_ext_par, base_par_extinction, &
175  ext_tss_extinction, tss_par_extinction, &
179 !
180 !-----------------------------------------------------------------------
181 !
182 !BEGIN
183  WRITE(lu,*) "*** USING AED2 COUPLING ***"
184 
185  old_zones = .true.
186  do_zone_averaging = .false.
187  aed2_nml_file = 'aed2.nml'
188  solution_method = 1
189  link_bottom_drag = .false.
190  link_surface_drag = .false.
191  link_water_density = .false.
192  link_water_clarity = .false.
193  link_solar_shade = .true.
194  link_rain_loss = .false.
195  link_ext_par = .false.
196  base_par_extinction = 0.1d0
197  ext_tss_extinction = .false.
198  tss_par_extinction = 0.02d0
199  do_2d_atm_flux = .true.
200  do_limiter = .false.
201  do_particle_bgc = .false.
202  dtaed2 = 0.d0
203 
204  IF ( aed2_init_core(dname).NE.0 ) THEN
205  stop "INITIALISATION OF AED2_CORE FAILED"
206  ENDIF
207  tname = trim(dname)//'aed2.nml'
208 ! TNAME = TRIM(DNAME)
209 
210  WRITE(lu,*)"READING AED2_MODELS CONFIG FROM ",trim(tname)
211  WRITE(lu,*) namlst
212 
213  OPEN(namlst,file=tname,action='READ',status='OLD',iostat=status)
214  IF ( status.NE.0 ) CALL stopit("CANNOT OPEN FILE " // trim(tname))
215  READ(namlst,nml=aed2_bio,iostat=status)
216  IF ( status.NE.0 ) stop "CANNOT READ NAMELIST ENTRY AED2_BIO"
217  kw = base_par_extinction
218  ksed = tss_par_extinction
219 
220  models = ''
221  READ(namlst, nml=aed2_models, iostat=status)
222  IF ( status.NE.0 ) stop "CANNOT READ NAMELIST ENTRY AED2_MODELS"
223 
224  IF ( do_zone_averaging ) old_zones = .false.
225 
226  DO i=1,SIZE(models)
227  IF (models(i).EQ.'') EXIT
228  CALL aed2_define_model(models(i), namlst)
229  ENDDO
230 
231  n_aed2_vars = aed2_core_status(nwq_var, nben_var, ndiag_var, n_sd)
232 
233 !#IF DEBUG
234  DO i=1,n_aed2_vars
235  IF ( aed2_get_var(i, tvar) ) THEN
236 ! WRITE(LU,*)"AED2 VAR ", I, TVAR%SHEET, TVAR%DIAG, TVAR%EXTERN, TVAR%FOUND, TRIM(TVAR%NAME),' ', TRIM(TVAR%MODEL_NAME)
237  ELSE
238  WRITE(lu,*)"AED2 VAR ", i, " IS EMPTY"
239  ENDIF
240  ENDDO
241 !#endif
242 
243  ndiag_var = ndiag_var + n_sd
244  n_vars = nwq_var
245  n_vars_ben = nben_var
246  n_vars_diag = ndiag_var
247  n_vars_diag_sheet = n_sd
248 
249 !#IF DEBUG
250 ! WRITE(LU,*)'AED2 INIT_AED2_MODELS : N_AED2_VARS = ',N_AED2_VARS,' NWQ_VAR = ',NWQ_VAR,' NBEN_VAR ',NBEN_VAR
251 !#endif
252 
253  CALL check_data
254 
255 !# NAMES = GRAB THE NAMES FROM INFO
256  ALLOCATE(names(1:nwq_var),stat=status)
257  IF (status.NE.0) THEN
258  stop 'ALLOCATE_MEMORY(): ERROR ALLOCATING (NAMES)'
259  ENDIF
260  ALLOCATE(bennames(1:nben_var),stat=status)
261  IF (status.NE.0) THEN
262  stop 'ALLOCATE_MEMORY(): ERROR ALLOCATING (BENNAMES)'
263  ENDIF
264  IF ( .NOT. ALLOCATED(diagnames) ) ALLOCATE(diagnames(ndiag_var))
265  IF (status.NE.0) THEN
266  stop 'ALLOCATE_MEMORY(): ERROR ALLOCATING (DIAGNAMES)'
267  ENDIF
268  ALLOCATE(min_(1:nwq_var+nben_var))
269  ALLOCATE(max_(1:nwq_var+nben_var))
270 
271  j = 0
272  DO i=1,n_aed2_vars
273  IF ( aed2_get_var(i, tvar) ) THEN
274  IF ( .NOT. (tvar%SHEET .OR. tvar%DIAG .OR. tvar%EXTERN) ) THEN
275  j = j + 1
276  names(j) = trim(tvar%NAME)
277  min_(j) = tvar%MINIMUM
278  max_(j) = tvar%MAXIMUM
279 ! WRITE(LU,*)"AED2 VAR NAME(",J,") : ", TRIM(NAMES(J))
280  ENDIF
281  ENDIF
282  ENDDO
283 
284  j = 0
285  DO i=1,n_aed2_vars
286  IF ( aed2_get_var(i, tvar) ) THEN
287  IF ( tvar%SHEET .AND. .NOT.(tvar%DIAG.OR.tvar%EXTERN) ) THEN
288  j = j + 1
289  bennames(j) = trim(tvar%NAME)
290  min_(nwq_var+j) = tvar%MINIMUM
291  max_(nwq_var+j) = tvar%MAXIMUM
292 ! WRITE(LU,*)"AED2 VAR_BEN NAME(",J,") : ", TRIM(BENNAMES(J))
293  ENDIF
294  ENDIF
295  ENDDO
296 
297  j = 0
298  DO i=1,n_aed2_vars
299  IF ( aed2_get_var(i, tvar) ) THEN
300  IF ( tvar%DIAG ) THEN
301  j = j + 1
302  diagnames(j) = trim(tvar%NAME)
303 ! WRITE(LU,*)"AED2 DIAG NAME(",J,") : ", TRIM(DIAGNAMES(J))
304  ENDIF
305  ENDIF
306  ENDDO
307 
308  CLOSE(namlst)
309  ! Intialising all pointers to null
310  cc => null()
311  cc_diag => null()
312  pp => null()
313  active => null()
314  surf_map => null()
315  benth_map => null()
316  lpar => null()
317  temp_aed2 => null()
318  salt => null()
319  rho_aed2 => null()
320  h_aed2 => null()
321  z_aed2 => null()
322  extcoeff => null()
323  tss => null()
324  bio_drag => null()
325  i_0 => null()
326  wnd => null()
327  air_temp => null()
328  rainaed2 => null()
329  area => null()
330  bathy => null()
331  shadefrac => null()
332  rainloss => null()
333  fsed_setl => null()
334  ustar_bed => null()
335  mat => null()
336 !
337 !-----------------------------------------------------------------------
338 !
339  RETURN
340  END SUBROUTINE init_aed2_models
341 
342 ! *************************
343  SUBROUTINE names_var_aed2 &
344 ! *************************
345 !
346  (nv_aed2,namesaed2)
347 !
348 !***********************************************************************
349 ! WAQTEL V7P3
350 !***********************************************************************
351 !
352 !brief Assigns the names of AED2 variables
353 !
354 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
355 !| NAMESAED2 |<--| NAMES OF AED2 VARIABLES
356 !| NAMLST |-->| UNIT FOR AED2 STEERING FILE
357 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
358 !
359 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
360 !
361  INTEGER, INTENT(IN) :: NV_AED2
362  CHARACTER(LEN=16),ALLOCATABLE,INTENT(OUT) :: NAMESAED2(:)
363 !
364 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
365 !
366  INTEGER I, J, STATUS
367  TYPE(aed2_variable_t),POINTER :: TVAR
368 !
369 !-----------------------------------------------------------------------
370 !
371  ALLOCATE(namesaed2(1:nv_aed2),stat=status)
372  IF(status.NE.0) stop 'ALLOCATE_MEMORY(): ERROR ALLOCATING (NAMES)'
373 
374  j = 0
375  DO i=1,n_aed2_vars
376  IF ( aed2_get_var(i, tvar) ) THEN
377  IF ( .NOT. (tvar%SHEET .OR. tvar%DIAG .OR. tvar%EXTERN) ) THEN
378  j = j + 1
379  namesaed2(j) = trim(tvar%NAME)
380  ENDIF
381  ENDIF
382  ENDDO
383 !
384 !-----------------------------------------------------------------------
385 !
386  RETURN
387  END SUBROUTINE names_var_aed2
388 
389 ! *******************************
390  SUBROUTINE init_var_aed2_models &
391 ! *******************************
392 !
393  (ncells, cc_, cc_diag_, nwq, nwqben, sm, bm)
394 !
395 !***********************************************************************
396 ! WAQTEL V7P3
397 !***********************************************************************
398 !
399 !brief Allocate memory for AED2 WQ variables
400 !
401 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
402 !| BM |-->| NUMBER OF BENTHIC NODES
403 !| CC_ |-->| VALUES OF WQ VARIABLES
404 !| CC_DIAG |-->| VALUES OF WQ DIAG VARIABLES
405 !| NCELLS |-->| NUMBER OF 3d NODES
406 !| NWQ |<->| NUMBER OF WQ VARIABLES
407 !| NWQBEN |<->| NUMER OF BENTHIC WQ VARIABLES
408 !| SM |-->| NUMBER OF SURFACE NODES
409 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
410 !
411 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
412 !
413  INTEGER,INTENT(IN) :: NCELLS
414  aed_real,POINTER,DIMENSION(:,:),INTENT(IN) :: cc_, cc_diag_
415  INTEGER,INTENT(INOUT) :: NWQ, NWQBEN
416  INTEGER,POINTER,DIMENSION(:),INTENT(IN) :: SM, BM
417 !
418 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
419 !
420  INTEGER RC, AV, V, SV, D, SD
421  TYPE(aed2_variable_t),POINTER :: TV
422 !
423 !-----------------------------------------------------------------------
424 !
425 !BEGIN
426  nwq = n_vars
427  nwqben = n_vars_ben
428 
429 ! WRITE(LU,*)'INIT_VAR_AED2_MODELS : NWQ = ',NWQ,' NWQBEN = ',NWQBEN
430 
431  cc => cc_
432  cc_diag => cc_diag_
433  surf_map => sm
434  benth_map => bm
435 
436 ! ALLOCATE STATE VARIABLE ARRAY
437  IF ( .NOT.ASSOCIATED(cc) ) stop ' ERROR : NO ASSOCIATION FOR (CC)'
438  cc = 0.d0
439 
440  IF (.NOT. ASSOCIATED(cc_diag) ) THEN
441  stop ' ERROR : NO ASSOCIATION FOR (CC_DIAG)'
442  ENDIF
443  cc_diag = 0.d0
444 ! ALLOCATE ARRAY WITH VERTICAL MOVEMENT RATES (M/S, POSITIVE FOR UPWARDS),
445 ! AND SET THESE TO THE VALUES PROVIDED BY THE MODEL.
446 
447  ALLOCATE(ws(1:ncells),stat=rc)
448  IF (rc.NE.0) stop 'ALLOCATE_MEMORY(): ERROR ALLOCATING (WS)'
449  ws = 0.d0
450 
451 !!# PLACE HOLDER FOR LAGRANIGAN PARTICLES
452 !IF(DO_PARTICLE_BGC) THEN
453 ! PP => PP_
454 !END IF
455 
456 ! ALLOCATE ARRAY FOR PHOTOSYNTHETICALLY ACTIVE RADIATION (PAR).
457 ! THIS WILL BE CALCULATED INTERNALLY DURING EACH TIME STEP.
458  ALLOCATE(par(1:ncells),stat=rc)
459  IF (rc.NE.0) stop 'ALLOCATE_MEMORY(): ERROR ALLOCATING (PAR)'
460  par = 0.d0
461  ALLOCATE(nir(1:ncells),stat=rc)
462  IF (rc.NE.0) stop 'ALLOCATE_MEMORY(): ERROR ALLOCATING (NIR)'
463  nir = 0.d0
464  ALLOCATE(uva(1:ncells),stat=rc)
465  IF (rc.NE.0) stop 'ALLOCATE_MEMORY(): ERROR ALLOCATING (UVA)'
466  uva = 0.d0
467  ALLOCATE(uvb(1:ncells),stat=rc)
468  IF (rc.NE.0) stop 'ALLOCATE_MEMORY(): ERROR ALLOCATING (UVB)'
469  uvb = 0.d0
470 
471 !# ALLOCATE ARRAY FOR SEDIMENTATION FLUXES AND INITIALIZE THESE TO ZERO (NO FLUX).
472  ALLOCATE(fsed_setl(1:ncells),stat=rc)
473  IF(rc.NE.0) stop 'ALLOCATE_MEMORY(): ERROR ALLOCATING (FSED_SETL)'
474  fsed_setl = 0.d0
475 
476 !# NOW SET INITIAL VALUES
477  v = 0
478  sv = 0
479  DO av=1,n_aed2_vars
480  IF ( .NOT. aed2_get_var(av, tv) ) THEN
481  stop "ERROR GETTING VARIABLE INFO"
482  ENDIF
483  IF ( .NOT. ( tv%EXTERN .OR. tv%DIAG) ) THEN !# NEITHER GLOBAL NOR DIAGNOSTIC VARIABLE
484  IF ( tv%SHEET ) THEN
485  sv = sv + 1
486  cc(n_vars+sv, :) = tv%INITIAL
487  ELSE
488  v = v + 1
489  cc(v,:) = tv%INITIAL
490  ENDIF
491  ENDIF
492  ENDDO
493 
494  IF ( init_values_file.NE.'' ) CALL set_initial_from_file
495 
496 ! CALL FV_INITIALIZE()
497 
498 
499  ALLOCATE(flux(n_vars, ncells),stat=rc)
500  IF (rc.NE.0) stop 'ALLOCATE_MEMORY(): ERROR ALLOCATING (FLUX)'
501 !#if !_NO_ODE_
502 ! ALLOCATE(FLUX2(N_VARS, NCELLS),STAT=RC) ; IF (RC.NE.0) STOP 'ALLOCATE_MEMORY(): ERROR ALLOCATING (FLUX2)'
503 ! ALLOCATE(FLUX3(N_VARS, NCELLS),STAT=RC) ; IF (RC.NE.0) STOP 'ALLOCATE_MEMORY(): ERROR ALLOCATING (FLUX3)'
504 ! ALLOCATE(FLUX4(N_VARS, NCELLS),STAT=RC) ; IF (RC.NE.0) STOP 'ALLOCATE_MEMORY(): ERROR ALLOCATING (FLUX4)'
505 ! ALLOCATE(CC1(N_VARS, NCELLS),STAT=RC) ; IF (RC.NE.0) STOP 'ALLOCATE_MEMORY(): ERROR ALLOCATING (CC1)'
506 !#endif
507 !
508 !-------------------------------------------------------------------------------
509  CONTAINS
510 
511 ! **************************
512  CHARACTER FUNCTION tolower &
513 ! **************************
514 !
515  (c)
516 !
517 !***********************************************************************
518 ! WAQTEL V7P3
519 !***********************************************************************
520 !
521 !brief
522 !
523 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
524 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
525 !
526 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
527 !
528  CHARACTER, INTENT(IN) :: C
529 !
530 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
531 !
532  INTEGER IC
533 !
534 !-----------------------------------------------------------------------
535 !
536 !BEGIN
537 !----------------------------------------------------------------------------
538  ic = ichar(c)
539  IF (ic.GE.65 .AND. ic.LT.90) ic = (ic+32)
540  tolower = char(ic)
541  END FUNCTION tolower
542 
543 ! *****************************************
544  FUNCTION same_str_icase(A, B) RESULT(RES)
545 ! *****************************************
546 !
547 !
548 !***********************************************************************
549 ! WAQTEL V7P3
550 !***********************************************************************
551 !
552 !brief Tests if A and B are the same string
553 !
554 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
555 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
556 !
557 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
558 !
559  CHARACTER(LEN=*), INTENT(IN) :: A,B
560 !
561 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
562 !
563  INTEGER LEN, I
564  LOGICAL RES
565 !
566 !-----------------------------------------------------------------------
567 !
568 !BEGIN
569  res = .false.
570  len = len_trim(a)
571  IF ( len.NE.len_trim(b) ) RETURN
572  DO i=1, len
573  IF (tolower(a(i:i)).NE.tolower(b(i:i)) ) RETURN
574  ENDDO
575  res = .true.
576 
577  END FUNCTION same_str_icase
578 
579 ! ********************************
580  SUBROUTINE set_initial_from_file
581 ! ********************************
582 !
583 !
584 !***********************************************************************
585 ! WAQTEL V7P3
586 !***********************************************************************
587 !
588 !brief Set initial values of AED2 variables from files
589 !
590 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
591 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
592 !
593 !
594 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
595 !
596 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
597 !
598  INTEGER UNIT, NCCOLS, CCOL
599  CHARACTER(LEN=32),POINTER,DIMENSION(:) :: CSVNAMES
600  aed_real,DIMENSION(:),ALLOCATABLE :: values
601  INTEGER :: IDX_COL = 0, numv = 0, numd = 0, t
602  INTEGER,DIMENSION(:),ALLOCATABLE :: VARS, VMAP
603  INTEGER,DIMENSION(:),ALLOCATABLE :: DVAR, DMAP
604  LOGICAL,DIMENSION(:),ALLOCATABLE :: VSHEET, DSHEET
605  LOGICAL MEH
606 !
607 !-----------------------------------------------------------------------
608 !
609 !BEGIN
610 !----------------------------------------------------------------------------
611  unit = aed_csv_read_header(init_values_file, csvnames, nccols)
612  WRITE(lu,*)'BENTHIC VARS INITIALISED FROM FILE : ', csvnames
613  IF (unit.LE.0) RETURN !# NO FILE FOUND
614  DO ccol=1,nccols
615  IF ( csvnames(ccol).EQ."ID" ) THEN
616  idx_col = ccol
617  EXIT
618  ENDIF
619  ENDDO
620 
621  ALLOCATE(vars(nccols)) ; ALLOCATE(vmap(nccols))
622  ALLOCATE(dvar(nccols)) ; ALLOCATE(dmap(nccols))
623  ALLOCATE(vsheet(nccols)) ; ALLOCATE(dsheet(nccols))
624  ALLOCATE(values(nccols))
625 
626  IF ( idx_col.GT.0 ) THEN
627  v = 0 ; sv = 0; d = 0; sd = 0
628  DO av=1,n_aed2_vars
629  IF ( .NOT. aed2_get_var(av, tv) ) THEN
630  stop "ERROR GETTING VARIABLE INFO"
631  ENDIF
632  IF ( .NOT. ( tv%EXTERN ) ) THEN !# DONT DO ENVIRONMENT VARS
633  DO ccol=1,nccols
634  IF ( same_str_icase(tv%NAME, csvnames(ccol)) ) THEN
635  IF (tv%DIAG) THEN
636  numd = numd + 1
637  dmap(numd) = ccol
638  dvar(numd) = numd
639  dsheet(numd) = tv%SHEET
640  ELSE
641  numv = numv + 1
642  vmap(numv) = ccol
643  IF ( tv%SHEET ) THEN
644  sv = sv + 1
645  vars(numv) = n_vars + sv
646  ELSE
647  v = v + 1
648  vars(numv) = v
649  ENDIF
650  vsheet(numd) = tv%SHEET
651  ENDIF
652  ENDIF
653  ENDDO
654  ENDIF
655  ENDDO
656 
657  DO WHILE ( aed_csv_read_row(unit, values) )
658  t = int(values(idx_col))
659  DO v=1,numv
660  IF ( vsheet(v) ) THEN
661  cc(vars(v), bm(t)) = values(vmap(v))
662  ELSE
663  cc(vars(v), sm(t):bm(t)) = values(vmap(v))
664  ENDIF
665  ENDDO
666  DO v=1,numd
667  IF ( vsheet(v) ) THEN
668  cc_diag(dvar(v), bm(t)) = values(dmap(v))
669  ELSE
670  cc_diag(dvar(v), sm(t):bm(t)) = values(dmap(v))
671  ENDIF
672  ENDDO
673  ENDDO
674  ENDIF
675 
676  meh = aed_csv_close(unit)
677 !# DO NOT CARE IF CLOSE FAILS
678 
679  IF (ASSOCIATED(csvnames)) DEALLOCATE(csvnames)
680  IF (ALLOCATED(values)) DEALLOCATE(values)
681  IF (ALLOCATED(vars)) DEALLOCATE(vars)
682  IF (ALLOCATED(vmap)) DEALLOCATE(vmap)
683  IF (ALLOCATED(dvar)) DEALLOCATE(dvar)
684  IF (ALLOCATED(dmap)) DEALLOCATE(dmap)
685 !
686 !-----------------------------------------------------------------------
687 !
688  RETURN
689  END SUBROUTINE set_initial_from_file
690 
691  END SUBROUTINE init_var_aed2_models
692 
693 ! ******************************
694  SUBROUTINE set_env_aed2_models &
695 ! ******************************
696 !
697  (dt_, & !
698 ! 3D ENV VARIABLES
699  temp_, &
700  salt_, &
701  rho_, &
702  h_, &
703  tss_, &
704  rad_, &
705 ! 3D FEEDBACK ARRAYS
706  extcoeff_, &
707 ! 2D ENV VARIABLES
708  area_, &
709  i_0_, &
710  wnd_, &
711  rain_, &
712  air_temp_, &
713  ustar_bed_, &
714  ustar_surf_, &
715  z_, &
716  bathy_, &
717  mat_id_, &
718  active_, &
719  benth_map_, &
720 ! 2D FEEDBACK ARRAYS
721  biodrag_, &
722  solarshade_, &
723  rainloss_ &
724  )
725 !
726 !***********************************************************************
727 ! WAQTEL V8P0
728 !***********************************************************************
729 !
730 !brief Provide information about TELEMAC-3D data not declared by
731 !+ AED2_BIO
732 !
733 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
734 !| ACTIVE_ |-->| COLUMN ACTIVE STATUS
735 !| AIR_TEMP_ |-->| AIR TEMPERATURE
736 !| AREA_ |-->| CELL AREA
737 !| BATHY_ |-->| HEIGHT OF COLUMN BOTTOM
738 !| BENTH_MAP_ |-->| NUMBER OF BENTHIC NODES
739 !| BIODRAG_ |-->| BIODRAG
740 !| DT_ |-->| TIME STEP
741 !| EXTCOEFF_ |-->| EXTINCTION COEFFICIENT
742 !| H_ |-->| CELL THICKNESS
743 !| I_0_ |-->| NET SURFACE IRRADIANCE
744 !| MAT_ID_ |-->| MATERIAL ID
745 !| RAD_ |-->| NET SHORT WAVE RADIATION
746 !| RAIN_ |-->| RAIN
747 !| RAINLOSS_ |-->| RAIN LOSS
748 !| RHO_ |-->| DENSITY
749 !| SALT_ |-->| SALINITY
750 !| SOLARSHADE_ |-->| SOLAR SHADE
751 !| TEMP_ |-->| TEMPERATURE
752 !| TSS_ |-->| TSS
753 !| USTAR_BED_ |-->| BED FRICTION VELOCITY
754 !| USTAR_SURF_ |-->| SURFACE FRICTION VELOCITY
755 !| WND_ |-->| 10M WIND SPEED
756 !| Z_ |-->| DEPTH
757 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
758 !
759 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
760 !
761  doubletype, INTENT(IN) :: dt_
762  aed_real, INTENT(IN), DIMENSION(:), POINTER :: temp_,salt_,rho_
763  aed_real, INTENT(IN), DIMENSION(:), POINTER :: h_,area_,tss_
764  aed_real, INTENT(IN), DIMENSION(:), POINTER :: extcoeff_, z_
765  aed_real, INTENT(IN), DIMENSION(:,:), POINTER :: rad_
766  aed_real, INTENT(IN), DIMENSION(:), POINTER :: i_0_,wnd_
767  aed_real, INTENT(IN), DIMENSION(:), POINTER :: ustar_bed_
768  aed_real, INTENT(IN), DIMENSION(:), POINTER :: rain_, bathy_ ! JC BATHY,
769  aed_real, INTENT(IN), DIMENSION(:), POINTER :: biodrag_
770  aed_real, INTENT(IN), DIMENSION(:), POINTER :: ustar_surf_
771  aed_real, INTENT(IN), DIMENSION(:), POINTER :: solarshade_
772  aed_real, INTENT(IN), DIMENSION(:), POINTER :: rainloss_ ! JC BATHY,
773  aed_real, INTENT(IN), DIMENSION(:), POINTER :: air_temp_
774  INTEGER, INTENT(IN), DIMENSION(:,:), POINTER :: MAT_ID_
775  INTEGER, INTENT(IN), DIMENSION(:), POINTER :: BENTH_MAP_
776  LOGICAL, INTENT(IN), DIMENSION(:), POINTER :: ACTIVE_
777 !
778 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
779 !
780  INTEGER I
781 !
782 !-----------------------------------------------------------------------
783 !
784 !BEGIN
785 !# PROVIDE POINTERS TO ARRAYS WITH ENVIRONMENTAL VARIABLES TO AED2.
786 
787 !# 2D (SHEET) VARIABLES BEING POINTED TO
788  area => area_
789  i_0 => i_0_
790  wnd => wnd_
791  ustar_bed => ustar_bed_
792  mat => mat_id_
793  bathy => bathy_
794  rainaed2 => rain_
795  shadefrac => solarshade_
796  rainloss => rainloss_
797  bio_drag => biodrag_
798  air_temp => air_temp_
799 
800 !# 3D VARIABLES BEING POINTED TO
801  h_aed2 => h_ !# LAYER HEIGHTS [1D ARRAY] NEEDED FOR ADVECTION, DIFFUSION
802  z_aed2 => z_ !# DEPTH [1D ARRAY], USED TO CALCULATE LOCAL PRESSURE
803  extcoeff => extcoeff_ !# BIOGEOCHEMICAL LIGHT ATTENUATION COEFFICIENTS [1D ARRAY],
804  !# OUTPUT OF BIOGEOCHEMISTRY, INPUT FOR PHYSICS
805  salt => salt_
806  temp_aed2 => temp_
807 
808  rho_aed2 => rho_
809  tss => tss_
810  active => active_
811 
812  benth_map => benth_map_
813 
814 !IF ( .NOT. ASSOCIATED(AREA) ) WRITE(LU,*) " NO ASSOCIATION FOR AREA"
815 !IF ( .NOT. ASSOCIATED(I_0) ) WRITE(LU,*) " NO ASSOCIATION FOR I_0"
816 !IF ( .NOT. ASSOCIATED(WND) ) WRITE(LU,*) " NO ASSOCIATION FOR WND"
817 !IF ( .NOT. ASSOCIATED(USTAR_BED) ) WRITE(LU,*) " NO ASSOCIATION FOR USTAR_BED"
818 !IF ( .NOT. ASSOCIATED(MAT) ) WRITE(LU,*) " NO ASSOCIATION FOR MAT"
819 !IF ( .NOT. ASSOCIATED(BATHY) ) WRITE(LU,*) " NO ASSOCIATION FOR BATHY"
820 !IF ( .NOT. ASSOCIATED(RAINAED2) ) WRITE(LU,*) " NO ASSOCIATION FOR RAINAED2"
821 !IF ( .NOT. ASSOCIATED(SHADEFRAC) ) WRITE(LU,*) " NO ASSOCIATION FOR SHADEFRAC"
822 !IF ( .NOT. ASSOCIATED(RAINLOSS) ) WRITE(LU,*) " NO ASSOCIATION FOR RAINLOSS"
823 !IF ( .NOT. ASSOCIATED(BIO_DRAG) ) WRITE(LU,*) " NO ASSOCIATION FOR BIO_DRAG"
824 !IF ( .NOT. ASSOCIATED(AIR_TEMP) ) WRITE(LU,*) " NO ASSOCIATION FOR AIR_TEMP"
825 !IF ( .NOT. ASSOCIATED(H) ) WRITE(LU,*) " NO ASSOCIATION FOR H"
826 !IF ( .NOT. ASSOCIATED(Z) ) WRITE(LU,*) " NO ASSOCIATION FOR Z"
827 !IF ( .NOT. ASSOCIATED(EXTCOEFF) ) WRITE(LU,*) " NO ASSOCIATION FOR EXTCOEFF"
828 !IF ( .NOT. ASSOCIATED(SALT) ) WRITE(LU,*) " NO ASSOCIATION FOR SALT"
829 !IF ( .NOT. ASSOCIATED(TEMP) ) WRITE(LU,*) " NO ASSOCIATION FOR TEMP"
830 !IF ( .NOT. ASSOCIATED(RHO) ) WRITE(LU,*) " NO ASSOCIATION FOR RHO"
831 !IF ( .NOT. ASSOCIATED(TSS) ) WRITE(LU,*) " NO ASSOCIATION FOR TSS"
832 !IF ( .NOT. ASSOCIATED(ACTIVE) ) WRITE(LU,*) " NO ASSOCIATION FOR ACTIVE"
833 
834  dtaed2 = dt_
835 
836  IF (link_ext_par) THEN
837  lpar => rad_(1,:)
838  ENDIF
839 
840  IF (old_zones) THEN
841 !# WE ALLOCATE THE FULL SIZE ARRAY BECAUSE THAT'S HOW THE INDICES ARE PRESENTED
842  ALLOCATE(zone(ubound(mat_id_,2))) ! JC
843  zone = 1.
844  DO i=1, ubound(mat_id_,2) ! JC
845 !# USE THE BOTTOM INDEX TO FILL THE ARRAY
846  zone(i) = mat_id_(1,i)! JC
847  ENDDO
848  ELSE
849  CALL init_zones(ubound(mat_id_, 2),mat_id_,n_aed2_vars,n_vars, &
850  n_vars_ben)
851  ENDIF
852 
853 ! WRITE(LU,*)'AED2 RE_INITIALIZE SKIPED STARTING ... '
854 !CALL RE_INITIALIZE() !
855 ! WRITE(LU,*)'AED2 RE_INITIALIZE SKIPED COMPLETED'
856 
857  CONTAINS
858 
859 !! ************************
860 ! SUBROUTINE RE_INITIALIZE
861 !! ************************
862 !!
863 !!***********************************************************************
864 !! WAQTEL V8P0
865 !!***********************************************************************
866 !!
867 !!brief Re initialize variables
868 !!
869 !!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
870 !!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
871 !!
872 !!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
873 !!
874 !!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
875 !! AED_REAL :: FLUX_BEN(N_VARS+N_VARS_BEN)
876 !! AED_REAL :: FLUX_ATM(N_VARS+N_VARS_BEN)
877 !! TYPE (AED2_COLUMN_T) :: COLUMN(N_AED2_VARS)
878 !
879 ! INTEGER COL, LEV, TOP, BOT, COUNT, NCOLS
880 !!
881 !!-----------------------------------------------------------------------
882 !!
883 !!BEGIN
884 ! NCOLS = UBOUND(ACTIVE, 1)
885 ! DO COL=1, NCOLS
886 ! TOP = SURF_MAP(COL)
887 ! BOT = BENTH_MAP(COL)
888 ! COUNT = TOP-BOT+1
889 !!!$ CALL DEFINE_COLUMN(COLUMN,COL,CC,CC_DIAG,FLUX,FLUX_ATM,FLUX_BEN)
890 !!! CTP TEMPORARY COMMENTED BECAUSE RE_INITIALIZE NOT CALLED
891 ! DO LEV=1, COUNT
892 !! CALL AED2_INITIALIZE(COLUMN, LEV) !!MJ
893 ! ENDDO
894 ! ENDDO
895 !!
896 !!-----------------------------------------------------------------------
897 !!
898 ! RETURN
899 ! END SUBROUTINE RE_INITIALIZE
900 !
901 !-----------------------------------------------------------------------
902 !
903  END SUBROUTINE set_env_aed2_models
904 
905 ! *********************
906  SUBROUTINE check_data
907 ! *********************
908 !
909 !
910 !***********************************************************************
911 ! WAQTEL V7P3
912 !***********************************************************************
913 !
914 !brief Check that all variable dependencies have been met
915 !
916 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
917 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
918 !
919 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
920 !
921 !-------------------------------------------------------------------------------
922 ! CHECK THAT ALL VARIABLE DEPENDENCIES HAVE BEEN MET
923 !-------------------------------------------------------------------------------
924 ! ARGUMENTS
925 !
926 ! LOCALS
927  INTEGER AV, V, D, SV, SD, EV, ERR_COUNT
928  TYPE(aed2_variable_t),POINTER :: TVAR
929 !-------------------------------------------------------------------------------
930 !BEGIN
931  v = 0 ; d = 0; sv = 0; sd = 0 ; ev = 0
932  err_count = 0
933 
934  DO av=1,n_aed2_vars
935  IF ( .NOT. aed2_get_var(av, tvar) ) THEN
936  stop "ERROR GETTING VARIABLE INFO"
937  ENDIF
938 
939  IF ( tvar%EXTERN ) THEN !# GLOBAL VARIABLE
940  ev = ev + 1
941  SELECT CASE (tvar%NAME)
942  CASE ( 'temperature' ) ; tvar%FOUND = .true.
943  CASE ( 'salinity' ) ; tvar%FOUND = .true.
944  CASE ( 'density' ) ; tvar%FOUND = .true.
945  CASE ( 'layer_ht' ) ; tvar%FOUND = .true.
946  CASE ( 'layer_area' ) ; tvar%FOUND = .true.
947  CASE ( 'rainaed2' ) ; tvar%FOUND = .true.
948  CASE ( 'rainloss' ) ; tvar%FOUND = .true.
949  CASE ( 'material' ) ; tvar%FOUND = .true.
950  CASE ( 'bathy' ) ; tvar%FOUND = .true.
951  CASE ( 'extc_coef' ) ; tvar%FOUND = .true.
952  CASE ( 'tss' ) ; tvar%FOUND = .true.
953  CASE ( 'par' ) ; tvar%FOUND = .true.
954  CASE ( 'nir' ) ; tvar%FOUND = .true.
955  CASE ( 'uva' ) ; tvar%FOUND = .true.
956  CASE ( 'uvb' ) ; tvar%FOUND = .true.
957  CASE ( 'sed_zone' ) ; tvar%FOUND = .true.
958  CASE ( 'wind_speed' ) ; tvar%FOUND = .true.
959  CASE ( 'par_sf' ) ; tvar%FOUND = .true.
960  CASE ( 'taub' ) ; tvar%FOUND = .true.
961  CASE ( 'air_temp' ) ; tvar%FOUND = .true.
962 ! CASE DEFAULT ; CALL STOPIT("ERROR: EXTERNAL VARIABLE "//TRIM(TVAR%NAME)//" NOT FOUND.")
963  END SELECT
964  ELSEIF ( tvar%DIAG ) THEN !# DIAGNOSTIC VARIABLE
965  IF ( tvar%SHEET ) THEN
966  sd = sd + 1
967  ELSE
968  d = d + 1
969  ENDIF
970  ELSE !# STATE VARIABLE
971  IF ( tvar%SHEET ) THEN
972  sv = sv + 1
973  ELSE
974  v = v + 1
975  ENDIF
976  ENDIF
977  IF ( .NOT. tvar%FOUND ) THEN
978  WRITE(lu,*) 'AV & VAR NAME ', av, tvar%NAME, tvar%FOUND
979  WRITE(lu,*) "ERROR: UNDEFINED VARIABLE ", tvar%NAME
980  err_count = err_count + 1
981  ENDIF
982  ENDDO
983 
984  IF ( n_vars.LT.v ) WRITE(lu,*)"MORE VARS THAN EXPECTED"
985  IF ( n_vars_ben.LT.sv ) WRITE(lu,*)"MORE SHEET VARS THAN EXPECTED"
986  IF(n_vars_diag.LT.sd+d) WRITE(lu,*)"MORE DIAG VARS THAN EXPECTED"
987  IF ( n_vars_diag_sheet.LT.sd ) THEN
988  WRITE(lu,*) "MORE SHEET DIAG VARS THAN EXPECTED"
989  ENDIF
990 
991  IF ( err_count.GT.0 ) CALL stopit("*** ERRORS IN CONFIGURATION")
992 !
993 !-----------------------------------------------------------------------
994 !
995  RETURN
996  END SUBROUTINE check_data
997 
998 ! ************************
999  SUBROUTINE define_column &
1000 ! ************************
1001 !
1002  (column, col, cc, cc_diag, flux_pel, flux_atm, flux_ben, nplan)
1004 !***********************************************************************
1005 ! WAQTEL V8P0
1006 !***********************************************************************
1007 !
1008 !brief Set column data structure from global arrays
1009 !
1010 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1011 !| CC |-->| WQ VARIABLES
1012 !| CC_DIAG |-->| WQ DIAG VARIABLES
1013 !| COL |-->| NODE
1014 !| COLUMN |<->| COLUMN DATA STRUCTURE
1015 !| FLUX_ATM |<->| ATMOSPHERIC FLUXES
1016 !| FLUX_BEN |<->| BENTHIC FLUXES
1017 !| FLUX_PEL |<->| PELAGIC FLUXES
1018 !| NPLAN |-->| NUMBER OF HORIZONTAL PLANES
1019 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1020 !
1021 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1022 !
1023  type(aed2_column_t), INTENT(INOUT) :: column(:)
1024  INTEGER, INTENT(IN) :: COL,NPLAN
1025  aed_real, TARGET, INTENT(IN) :: cc(:,:) !# (N_VARS, N_LAYERS)
1026  aed_real, TARGET, INTENT(IN) :: cc_diag(:,:) !# (N_VARS, N_LAYERS)
1027  aed_real, TARGET, INTENT(INOUT) :: flux_pel(:,:) !# (N_VARS, N_LAYERS)
1028  aed_real, TARGET, INTENT(INOUT) :: flux_atm(:) !# (N_VARS)
1029  aed_real, TARGET, INTENT(INOUT) :: flux_ben(:) !# (N_VARS)
1030 !
1031 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1032 !
1033  INTEGER AV, TOP, BOT, V, D, SV, SD, EV
1034  TYPE(aed2_variable_t),POINTER :: TVAR
1035 !
1036 !-----------------------------------------------------------------------
1037 !
1038 !BEGIN
1039  top = surf_map(col)
1040  bot = benth_map(col)
1041 
1042  v = 0 ; d = 0; sv = 0; sd = 0 ; ev = 0
1043  DO av=1,n_aed2_vars
1044 
1045  IF (.NOT.aed2_get_var(av,tvar)) THEN
1046  stop "ERROR GETTING VARIABLE INFO"
1047  ENDIF
1048 
1049  IF ( tvar%EXTERN ) THEN !# GLOBAL VARIABLE
1050  ev = ev + 1
1051  SELECT CASE (tvar%NAME)
1052  CASE ( 'temperature' )
1053  column(av)%CELL => temp_aed2(top+(col-1)*nplan:bot+(col-1)*nplan)
1054  CASE ( 'salinity' )
1055  column(av)%CELL => salt(top+(col-1)*nplan:bot+(col-1)*nplan)
1056  CASE ( 'density' )
1057  column(av)%CELL => rho_aed2(top+(col-1)*nplan:bot+(col-1)*nplan)
1058  CASE ( 'layer_ht' )
1059  column(av)%CELL => h_aed2(top+(col-1)*nplan:bot+(col-1)*nplan)
1060 !!$ CASE ( 'temperature' )
1061 !!$ COLUMN(AV)%CELL => TEMP_AED2(TOP:BOT)
1062 !!$ CASE ( 'salinity' ) ; COLUMN(AV)%CELL => SALT(TOP:BOT)
1063 !!$ CASE ( 'densigy' ) ; COLUMN(AV)%CELL => RHO_AED2(TOP:BOT)
1064 !!$ CASE ( 'layer_ht' ) ; COLUMN(AV)%CELL => H_AED2(TOP:BOT)
1065  CASE ( 'layer_area' ) ; column(av)%CELL_SHEET => area(col)
1066  CASE ( 'rainaed2' )
1067  column(av)%CELL_SHEET => rainaed2(col)! DOES IT MEAN THIS VAR IS ALREADY POINTED TO THE CORRECT VARIABLE OF TFFV? JC
1068  CASE ( 'rainloss' )
1069  column(av)%CELL_SHEET => rainloss(col)! JC
1070  CASE ( 'material' ) ; column(av)%CELL_SHEET => zone(col)
1071  CASE ( 'bathy' ) ; column(av)%CELL_SHEET => bathy(col)! DOES IT MEAN THIS VAR IS ALREADY POINTED TO THE CORRECT VARIABLE OF TFFV? JC
1072  CASE ( 'extc_coef' )
1073  column(av)%CELL => extcoeff(top:bot)
1074  CASE ( 'tss' ) ; column(av)%CELL => tss(top:bot)
1075  CASE ( 'par' ) ; IF (link_ext_par) THEN
1076  column(av)%CELL => lpar(top:bot)
1077  ELSE
1078  column(av)%CELL => par(top+(col-1)*nplan:bot+(col-1)*nplan)
1079 !!$ COLUMN(AV)%CELL => PAR(TOP:BOT)
1080  ENDIF
1081  CASE ( 'nir' )
1082  column(av)%CELL => nir(top+(col-1)*nplan:bot+(col-1)*nplan)
1083  CASE ( 'uva' )
1084  column(av)%CELL => uva(top+(col-1)*nplan:bot+(col-1)*nplan)
1085  CASE ( 'uvb' )
1086  column(av)%CELL => uvb(top+(col-1)*nplan:bot+(col-1)*nplan)
1087 !!$ CASE ( 'nir' ) ; COLUMN(AV)%CELL => NIR(TOP:BOT)
1088 !!$ CASE ( 'uva' ) ; COLUMN(AV)%CELL => UVA(TOP:BOT)
1089 !!$ CASE ( 'uvb' ) ; COLUMN(AV)%CELL => UVB(TOP:BOT)
1090  CASE ( 'sed_zone' ) ; column(av)%CELL_SHEET => zone(col)
1091  CASE ( 'wind_speed' ) ; column(av)%CELL_SHEET => wnd(col)
1092  CASE ( 'par_sf' ) ; column(av)%CELL_SHEET => i_0(col)
1093  CASE ( 'taub' ) ; column(av)%CELL_SHEET => col_taub
1094  CASE ( 'air_temp' )
1095  column(av)%CELL_SHEET => air_temp(col)
1096  CASE DEFAULT
1097  CALL stopit("ERROR: EXTERNAL VARIABLE " &
1098  //trim(tvar%NAME)//" NOT FOUND.")
1099  END SELECT
1100  ELSEIF ( tvar%DIAG ) THEN !# DIAGNOSTIC VARIABLE
1101  d = d + 1
1102  IF ( tvar%SHEET ) THEN
1103  column(av)%CELL_SHEET => cc_diag(d, bot)
1104  ELSE
1105 ! ABR TO FILL THE WHOLE ARRAY, SO THAT IT CAN BE USED IN OUTPUT
1106  column(av)%CELL => cc_diag(d,top+(col-1)*nplan:bot+(col-1)*nplan)
1107  ENDIF
1108  ELSE !# STATE VARIABLE
1109  IF ( tvar%SHEET ) THEN
1110  sv = sv + 1
1111  IF ( tvar%BOT ) THEN
1112  column(av)%CELL_SHEET => cc(n_vars+sv, bot)
1113  ELSEIF ( tvar%TOP ) THEN
1114  column(av)%CELL_SHEET => cc(n_vars+sv, top)
1115  ENDIF
1116  column(av)%FLUX_BEN => flux_ben(n_vars+sv)
1117  column(av)%FLUX_ATM => flux_atm(n_vars+sv)
1118  ELSE
1119  v = v + 1
1120  column(av)%CELL => cc(v,top:bot)
1121  column(av)%FLUX_PEL => flux_pel(v,top:bot)
1122  column(av)%FLUX_BEN => flux_ben(v)
1123  column(av)%FLUX_ATM => flux_atm(v)
1124  ENDIF
1125  ENDIF
1126  ENDDO
1127 !
1128 !-----------------------------------------------------------------------
1129 !
1130  RETURN
1131  END SUBROUTINE define_column
1132 
1133 ! ****************************
1134  SUBROUTINE calculate_fluxes &
1135 ! ****************************
1136 !
1137  (column, count, flux_pel, flux_atm, flux_ben, h_aed2)
1139 !***********************************************************************
1140 ! WAQTEL V8P0
1141 !***********************************************************************
1142 !
1143 !brief Checks the current values of all state variables and repairs
1144 !+ these
1145 !
1146 !history M. JODEAU + J. VIDAL HURTADO (LNHE)
1147 !+ 20/02/2018
1148 !+ V8P0
1149 !+ Limitation of fluxes to prevent from negative concentrations
1150 !+ Code similar to glm_aed2.F90
1151 !
1152 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1153 !| COLUMN |<->| COLUMN DATA STRUCTURE
1154 !| COUNT |-->| NUMBER OF LAYERS
1155 !| FLUX_ATM |<->| ATMOSPHERIC FLUXES
1156 !| FLUX_BEN |<->| BENTHIC FLUXES
1157 !| FLUX_PEL |<->| PELAGIC FLUXES
1158 !| H_AED2 |-->| CELL THICKNESS
1159 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1160 !
1161 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1162 !
1163  type(aed2_column_t), INTENT(INOUT) :: column(:)
1164  INTEGER, INTENT(IN) :: COUNT
1165  aed_real, INTENT(INOUT) :: flux_pel(:,:) !# (N_VARS, N_LAYERS)
1166  aed_real, INTENT(INOUT) :: flux_atm(:) !# (N_VARS)
1167  aed_real, INTENT(INOUT) :: flux_ben(:) !# (N_VARS)
1168  aed_real, INTENT(INOUT) :: h_aed2(:) !# (N_LAYERS)
1169 !
1170 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1171 !
1172  INTEGER I
1173 !
1174 !-----------------------------------------------------------------------
1175 !
1176 !BEGIN
1177  flux_pel = zero_
1178  flux_atm = zero_
1179  flux_ben = zero_
1180 
1181 !# CALCULATE TEMPORAL DERIVATIVES DUE TO AIR-WATER EXCHANGE.
1182  CALL aed2_calculate_surface(column, 1)
1183 
1184 !# DISTRIBUTE THE FLUXES INTO PELAGIC SURFACE LAYER
1185  IF ( do_2d_atm_flux .OR. count.GT.1 ) THEN
1186  IF (h_aed2(1).GT.eps_aed2) THEN !!ajout mj
1187  flux_pel(:,1) = flux_pel(:,1) + flux_atm(:)/h_aed2(1)
1188  ENDIF
1189  ENDIF
1190 
1191  IF ( .NOT. do_zone_averaging ) THEN
1192 !# CALCULATE TEMPORAL DERIVATIVES DUE TO BENTHIC EXCHANGE PROCESSES.
1193  CALL aed2_calculate_benthic(column, count)
1194 !# DISTRIBUTE BOTTOM FLUX INTO PELAGIC OVER BOTTOM BOX (I.E., DIVIDE BY LAYER HEIGHT).
1195  IF (h_aed2(count).GT.eps_aed2) THEN !!ajout mj
1196 ! FLUX_PEL(:,COUNT) = FLUX_PEL(:,COUNT)/H_AED2(COUNT)
1197 ! LIMIT FLUXES DEPENDING ON THE CONCENTRATION, CODE SIMILAR TO glm_aed2.F90
1198  flux_pel(:,count) = max( -1.0*cc(:,count) , flux_pel(:,count)/h_aed2(count) )
1199  ENDIF
1200  ENDIF
1201 
1202 ! LIMIT FLUXES DEPENDING ON THE CONCENTRATION, CODE SIMILAR TO glm_aed2.F90
1203 ! DONE ABOVE
1204 ! DO I=2,COUNT
1205 ! IF (H_AED2(I).GT.0.D0) THEN
1206 ! FLUX_PEL(:,I) = MAX ( -1.0*CC(:,I) , FLUX_PEL(:,I)/H_AED2(I) ) ! JVH -->
1207 ! FLUX_PEL(:,I) = MAX ( -1.0*CC(:,I) , FLUX_PEL(:,I) )
1208 ! ENDIF
1209 ! ENDDO
1210 
1211 !# ADD PELAGIC SINK AND SOURCE TERMS FOR ALL DEPTH LEVELS.
1212  DO i=1,count
1213  CALL aed2_calculate(column, i)
1214  ENDDO
1215 !
1216 !-----------------------------------------------------------------------
1217 !
1218  RETURN
1219  END SUBROUTINE calculate_fluxes
1220 
1221 ! ***********************
1222  SUBROUTINE check_states &
1223 ! ***********************
1224 !
1225  (column, top, bot)
1227 !***********************************************************************
1228 ! WAQTEL V7P3
1229 !***********************************************************************
1230 !
1231 !brief Checks states
1232 !
1233 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1234 !| BOT |-->| NUMBER OF THE BOTTOM LAYER
1235 !| COLUMN |<->| COLUMN DATA STRUCTURE
1236 !| TOP |-->| NUMBER OF THE SURFACE LAYER
1237 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1238 !
1239 ! USES
1240 ! USE IEEE_ARITHMETIC ! MODIF MJ
1241 !
1242 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1243 !
1244  type(aed2_column_t),INTENT(INOUT) :: column(:)
1245  INTEGER,INTENT(IN) :: TOP, BOT
1246 !
1247 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1248 !
1249  TYPE(aed2_variable_t),POINTER :: TV
1250  INTEGER I,V,LEV
1251 !
1252 !-----------------------------------------------------------------------
1253 !
1254 !BEGIN
1255  DO lev=top, bot
1256 !CALL AED2_EQUILIBRATE(COLUMN, LEV)
1257  v = 0
1258  DO i=1,n_aed2_vars
1259  IF ( aed2_get_var(i, tv) ) THEN
1260  IF ( .NOT. (tv%DIAG .OR. tv%EXTERN) ) THEN
1261  v = v + 1
1262  IF ( do_limiter ) THEN
1263  IF ( .NOT. isnan(min_(v)) ) THEN
1264  IF ( cc(v, lev).LT.min_(v) ) cc(v, lev) = min_(v)
1265  ENDIF
1266  IF ( .NOT. isnan(max_(v)) ) THEN
1267  IF ( cc(v, lev).GT.max_(v) ) cc(v, lev) = max_(v)
1268  ENDIF
1269  ENDIF
1270  ENDIF
1271  ENDIF
1272  ENDDO
1273  ENDDO
1274 !
1275 !-----------------------------------------------------------------------
1276 !
1277  RETURN
1278  END SUBROUTINE check_states
1279 
1280 ! *************************
1281  SUBROUTINE do_aed2_models &
1282 ! *************************
1283 !
1284  (ncells, ncols, nplan, fluxaed2, extcaed2, ta)
1286 !***********************************************************************
1287 ! WAQTEL V8P0
1288 !***********************************************************************
1289 !
1290 !brief Calculates source terms
1291 !
1292 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1293 !| EXTCAED2 |<--| EXTINCTION COEFFICIENT
1294 !| FLUXAED2 |<--| SOURCES TERMS
1295 !| NCELLS |-->| NUMBER OF 3D NODES
1296 !| NCOLS |-->| NUMBER OF 2D NODES
1297 !| NPLAN |-->| NUMBER OF HORIZONTAL PLANES
1298 !| TA |-->| TRACERS
1299 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1300 !
1301 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1302 !
1303  INTEGER, INTENT(IN) :: NCELLS, NCOLS, NPLAN
1304 ! NO DEFAULT INITIALISATION FOR USER TYPE COMPONENTS ALLOWED
1305  DOUBLE PRECISION , INTENT(INOUT) :: FLUXAED2(:,:,:)
1306  DOUBLE PRECISION , INTENT(OUT ) :: EXTCAED2(:,:)
1307  TYPE(bief_obj), INTENT(IN) :: TA
1308 !
1309 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1310 !
1311  TYPE(aed2_variable_t),POINTER :: TV
1312 
1313  aed_real :: flux_ben(n_vars+n_vars_ben)
1314  aed_real :: flux_atm(n_vars+n_vars_ben)
1315  type(aed2_column_t) :: column(n_aed2_vars)
1316 
1317  INTEGER I, J, COL, LEV, TOP, BOT, V
1318 !
1319 !-----------------------------------------------------------------------
1320  fluxaed2(:,:,:)= 0.d0
1321  extcaed2(:,:) = kw
1322 !
1323 !BEGIN
1324 
1325 ! DO_ZONE_AVERAGING = .FALSE.
1326  IF ( do_zone_averaging ) THEN
1327  IF (link_ext_par) THEN
1328  CALL calc_zone_areas(ncols,temp_aed2,salt,h_aed2,area,wnd, &
1329  rho_aed2,extcoeff,i_0,par,tss,active, &
1330  rainaed2)
1331  ELSE
1332  CALL calc_zone_areas(ncols,temp_aed2,salt,h_aed2,area,wnd, &
1333  rho_aed2,extcoeff,i_0,lpar,tss,active, &
1334  rainaed2)
1335  ENDIF
1336  ENDIF
1337 
1338 !OMP DO
1339 !#--------------------------------------------------------------------
1340 !# LOOP THROUGH COLUMNS DOING JOBS PRIOR TO THE KINETICS BEING SOLVED
1341 
1342  DO col=1, ncols
1343  IF (.NOT. active(col)) cycle
1344 
1345  top = surf_map(col)
1346  bot = benth_map(col)
1347 
1348  CALL define_column(column,col,cc,cc_diag,flux,flux_atm, &
1349  flux_ben,nplan)
1350 
1351 !# FIRSTLY RUN THROUGH ALL VARS AND SELECT STATE VARS
1352  v = 0
1353  DO i=1,n_aed2_vars
1354  IF ( aed2_get_var(i, tv) ) THEN
1355  IF ( .NOT. (tv%SHEET .OR. tv%DIAG .OR. tv%EXTERN) ) THEN
1356  v = v + 1
1357 !# ONLY FOR STATE_VARS THAT ARE NOT SHEET
1358  IF ( .NOT. isnan(tv%MOBILITY) ) THEN
1359  ws(top:bot) = tv%MOBILITY
1360  CALL settling(bot-top+1,dtaed2, &
1361  h_aed2(top+(col-1)*nplan: &
1362  bot+(col-1)*nplan), &
1363  ws(top:bot),fsed_setl(col),column(i)%CELL)
1364 ! WS TO BE DEPENDENT ON COL!
1365  ENDIF
1366  ENDIF
1367  ENDIF
1368  ENDDO
1369  CALL check_states(column, top, bot)
1370 
1371 !# POPULATE LOCAL LIGHT/EXTC ARRAYS ONE COLUMN AT A TIME
1372  IF (.NOT. link_ext_par) & !#MH CHECK LINK_EXT_PAR LOGIC
1373  CALL light(column,bot-top+1,i_0(col),extcoeff(top:bot), &
1374  par(top+(col-1)*nplan:bot+(col-1)*nplan), &
1375  h_aed2(top+(col-1)*nplan:bot+(col-1)*nplan))
1376 
1377 !# NON-PAR LIGHT FUDGE
1378  DO i=top,bot
1379  nir(i+(col-1)*nplan) = (par(i+(col-1)*nplan)/0.45d0)*0.510d0
1380  uva(i+(col-1)*nplan) = (par(i+(col-1)*nplan)/0.45d0)*0.035d0
1381  uvb(i+(col-1)*nplan) = (par(i+(col-1)*nplan)/0.45d0)*0.005d0
1382  ENDDO
1383 !!$ NIR(TOP:BOT) = (PAR(TOP:BOT)/0.45D0) * 0.510D0
1384 !!$ UVA(TOP:BOT) = (PAR(TOP:BOT)/0.45D0) * 0.035D0
1385 !!$ UVB(TOP:BOT) = (PAR(TOP:BOT)/0.45D0) * 0.005D0
1386  ENDDO
1387 !OMP END DO
1388 
1389  IF ( do_zone_averaging ) THEN
1390  CALL copy_to_zone(ncols, cc, area, active, benth_map)
1392  ENDIF
1393 
1394 !OMP DO
1395 
1396 !#--------------------------------------------------------------------
1397 !# THIS IS THE MAIN WQ SOLUTION LOOP
1398  DO col=1, ncols
1399 
1400 !# FIND TOP AND BOTTOM CELL INDICES BASED ON MAPS PROVIDED BY THE HOST
1401  top = surf_map(col)
1402  bot = benth_map(col)
1403 
1404 !# COMPUTE BOTTOM SHEAR STRESS FOR THIS COLUMN BASED ON USTAR FROM HOST
1405  col_taub = rho_aed2(bot)*(ustar_bed(col)*ustar_bed(col))
1406 
1407  DO j=1,n_vars
1408  DO lev=top,bot
1409 !# IF BOT .NE. NPLAN, IT IS NOT THE TOP FOR TA%ADR()%P%R
1410 !!$ CC(J,LEV) = TA%ADR(J+2)%P%R(COL+(BOT-LEV)*NCOLS)
1411  cc(j,lev) = ta%ADR(j+2)%P%R(col+(nplan-lev)*ncols)
1412  ENDDO
1413  ENDDO
1414 
1415 !# SET COLUMN DATA STRUCTURE FROM GLOBAL ARRAYS
1416  CALL define_column(column,col,cc,cc_diag,flux,flux_atm, &
1417  flux_ben,nplan)
1418 
1419 ! IF (.NOT. ACTIVE(COL)) THEN
1420 ! CALL AED2_CALCULATE_RIPARIAN(COLUMN, BOT-TOP+1, ZERO_);
1421 ! CALL AED2_CALCULATE_DRY(COLUMN, BOT-TOP+1);
1422 ! CALL AED2_RAIN_LOSS(COLUMN, BOT-TOP+1, RAIN_LOSS);! JC
1423 ! CYCLE
1424 ! ELSE
1425 ! CALL AED2_CALCULATE_RIPARIAN(COLUMN, BOT-TOP+1, 1.D0);
1426 ! ENDIF
1427 
1428 !#IF _NO_ODE_
1429 !# FOR THIS COLUMN, DO THE MAIN KINETIC/BGC FLUX CALCULATION
1430 !# (THIS INCLUDES WATER COLUMN, SURFACE AND BENTHIC INTERFACES)
1431  CALL calculate_fluxes(column,bot-top+1,flux(:,top:bot), &
1432  flux_atm,flux_ben, &
1433  h_aed2(top+(col-1)*nplan: &
1434  bot+(col-1)*nplan))
1435 
1436 !# FIND THE PARTICLES IN THIS COLUMN AND UPDATE PARTICLE BGC
1437  IF(do_particle_bgc) THEN
1438  CALL particles(column,bot-top+1, &
1439  h_aed2(top+(col-1)*nplan:bot+(col-1)*nplan))
1440  ENDIF
1441 
1442 !# DO RIPARIAN INTERFACES FOR THIS COLUMN AND UPDATE FLUXES
1443 ! IF ( .NOT. RIPARIAN(COLUMN, ACTIVE(COL), SHADEFRAC(COL), RAINLOSS(COL)) ) &
1444 ! CYCLE !!! COMMENTAIRE MJ
1445 
1446 !# NOW GO FORTH AND SOLVE THE ODE (EULER! - LINK TO FV_ODE WOULD BE NICE)
1447  DO lev = top, bot
1448  DO i = 1, n_vars
1449 ! EXPLICIT IF FLUX IS POSITIVE
1450 ! IMPLICIT (USING PATANKAR, 1980) IF THE FLUX IS NEGATIVE
1451  IF(flux(i,lev).GE.0.d0) THEN
1452  cc(i,lev)=cc(i,lev)+dtaed2*flux(i,lev)
1453  ELSE
1454  cc(i,lev)=cc(i,lev)/(1.d0-dtaed2*flux(i,lev)/max(cc(i,lev),1.d-15))
1455  ENDIF
1456 !#IF DEBUG>1
1457 !# CHECK FOR NANS
1458  IF ( isnan(cc(i,lev)) ) THEN
1459  WRITE(lu,*)'NAN AT I = ', i, ' LEV = ', lev
1460  WRITE(lu,*)'H_AED2(LEV) = ', h_aed2(lev+(col-1)*nplan), &
1461  ' FLUX(I,LEV) = ', flux(i,lev)
1462  WRITE(lu,*)'TOP OF COLUMN @ ',top,' BOTTOM OF COLUMN @ ',&
1463  bot
1464  CALL stopit('NAN VALUE')
1465  ENDIF
1466 
1467 !#endif _NO_ODE
1468 
1469  ENDDO ! VARS
1470  ENDDO ! LEVELS
1471 
1472  IF ( do_zone_averaging ) THEN ! UNTESTED
1473  DO i = n_vars+1, n_vars+n_vars_ben
1474  cc(i,bot)=cc(i,bot)+dtaed2*flux(i,bot)
1475 
1476 !#IF DEBUG>1
1477 !# CHECK FOR NANS
1478  IF ( isnan(cc(i,bot)) ) THEN
1479  WRITE(lu,*)'NAN AT I = ', i, ' BOT = ', bot
1480  WRITE(lu,*)'H_AED2(BOT) = ',h_aed2(bot+(col-1)*nplan), &
1481  ' FLUX(I,BOT) = ', flux(i,bot)
1482  CALL stopit('NAN VALUE')
1483  ENDIF
1484 !#endif
1485  ENDDO ! BEN VARS
1486  ENDIF
1487 !#endif
1488 
1489 !# DO NON-KINETIC UPDATES TO BGC VARIABLES (EQ EQUILIBRATION)
1490  CALL update(column, bot-top+1)
1491 
1492 !# NOW THE BGC UPDATES ARE COMPLETE, UPDATE LINKS TO HOST MODEL
1493  CALL biodrag(column, bot-top+1, bio_drag(col))
1494  CALL bioextinction(column, bot-top+1, extcoeff(top:bot))
1495 !CALL BIODENSITY()
1496 
1497  CALL check_states(column, top, bot)
1498 
1499 ! EXPORT OF SOURCE TERMS
1500  DO lev = top, bot
1501  DO i = 1, n_vars
1502 ! EXPLICIT IF FLUX IS POSITIVE
1503 ! IMPLICIT (USING PATANKAR, 1980) IF THE FLUX IS NEGATIVE
1504  IF(flux(i,lev).GE.0.d0) THEN
1505  fluxaed2(col,i,nplan-lev+1) = flux(i,lev)
1506  ELSE
1507  fluxaed2(col,i,nplan-lev+1) = flux(i,lev)/max(cc(i,lev),1.d-15)
1508  ENDIF
1509  ENDDO !N_VARS
1510  extcaed2(col,nplan-lev+1) = kw + extcoeff(lev)
1511  ENDDO !LEV
1512  ENDDO ! COLS
1513 !OMP END DO
1514 
1515  IF ( do_zone_averaging ) THEN
1516  CALL copy_to_zone(ncols, cc, area, active, benth_map)
1517  ENDIF
1518 !
1519 !-----------------------------------------------------------------------
1520 !
1521  RETURN
1522  END SUBROUTINE do_aed2_models
1523 
1524 ! ****************************
1525  SUBROUTINE clean_aed2_models
1526 ! ****************************
1527 !
1528 !
1529 !***********************************************************************
1530 ! WAQTEL V7P3
1531 !***********************************************************************
1532 !
1533 !brief Cleans AED2 models
1534 !
1535 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1536 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1537 !
1538 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1539 !
1540 !BEGIN
1541 ! DEALLOCATE INTERNAL ARRAYS
1542  IF (ALLOCATED(ws)) DEALLOCATE(ws)
1543  IF (ALLOCATED(total)) DEALLOCATE(total)
1544  IF (ALLOCATED(nir)) DEALLOCATE(nir)
1545  IF (ALLOCATED(par)) DEALLOCATE(par)
1546  IF (ALLOCATED(uva)) DEALLOCATE(uva)
1547  IF (ALLOCATED(uvb)) DEALLOCATE(uvb)
1548  IF (ALLOCATED(min_)) DEALLOCATE(min_)
1549  IF (ALLOCATED(max_)) DEALLOCATE(max_)
1550  IF (ALLOCATED(flux)) DEALLOCATE(flux)
1551  IF (ALLOCATED(zone)) DEALLOCATE(zone)
1552 !
1553 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1554 !
1555 !
1556 !-----------------------------------------------------------------------
1557 !
1558  RETURN
1559  END SUBROUTINE clean_aed2_models
1560 
1561 ! ****************
1562  SUBROUTINE light &
1563 ! ****************
1564  (column, count, io, extc, par_, h_)
1565 !-------------------------------------------------------------------------------
1566 !
1567 !brief Calculate photosynthetically active radiation over the entire
1568 !+ domain based on surface radiation, and background and biotic
1569 !+ extention.
1570 !
1571 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1572 !| COLUMN |<->| COLUMN DATA STRUCTURE
1573 !| COUNT |-->| NUMBER OF PLANES
1574 !| EXTC |<->| EXTINCTION COEFFICIENT
1575 !| H_ |<->| LAYER THICKNESS
1576 !| I0 |-->| NET SURFACE IRRADIANCE
1577 !| PAR_ |<->| NET SHORT WAVE RADIATION
1578 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1579 !
1580 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1581 !
1582  type(aed2_column_t), INTENT(INOUT) :: column(:)
1583  INTEGER, INTENT(IN) :: COUNT
1584  aed_real, INTENT(IN) :: io
1585  aed_real, INTENT(INOUT) :: extc(:)
1586  aed_real, INTENT(INOUT) :: par_(:)
1587  aed_real, INTENT(INOUT) :: h_(:)
1588 !
1589 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1590 !
1591  INTEGER I
1592  aed_real zz, localext
1593 !
1594 !-----------------------------------------------------------------------
1595 !
1596 !BEGIN
1597  zz = zero_
1598  localext = zero_
1599 
1600  CALL bioextinction(column,count,extc)
1601 
1602  localext = extc(1)
1603  zz = 0.001d0 !0.5*H_(1) !MH: ASSUME TOP OF LAYER
1604  par_(1) = 0.45d0 * io * exp( -(kw+localext) * zz )
1605 
1606  IF (count.LE.1) RETURN
1607 
1608  DO i = 2, count
1609  localext = extc(i)
1610 
1611 !ZZ = ZZ + 0.5D0*H_(I)
1612  zz = h_(i)
1613  par_(i) = par_(i-1) * exp( -(kw+localext) * zz )
1614  ENDDO
1615 !
1616 !-----------------------------------------------------------------------
1617 !
1618  RETURN
1619  END SUBROUTINE light
1620 
1621 ! *******************
1622  SUBROUTINE settling &
1623 ! *******************
1624  (n,dtaed2,h_aed2,ww,fsed,y)
1626 !***********************************************************************
1627 ! WAQTEL V7P3
1628 !***********************************************************************
1629 !
1630 !brief Update settling of AED2 state variables in a given column
1631 !
1632 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1633 !| DTAED2 |-->| TIME STEP
1634 !| FSED |<->| VALUE OF SEDIMENT FLUX DUE TO SETTLING
1635 !| H_AED2 |-->| LAYER THICKNESS
1636 !| N |-->| NUMBER OF VERTICAL LAYERS
1637 !| WW |-->| VERTICAL ADVECTION SPEED
1638 !| Y |<->|
1639 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1640 !
1641 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1642 !
1643  INTEGER,INTENT(IN) :: N !# NUMBER OF VERTICAL LAYERS
1644  aed_real,INTENT(IN) :: dtaed2 !# TIME STEP (S)
1645  aed_real,INTENT(IN) :: h_aed2(:) !# LAYER THICKNESS (M)
1646  aed_real,INTENT(IN) :: ww(:) !# VERTICAL ADVECTION SPEED
1647  aed_real,INTENT(INOUT) :: fsed !# VALUE OF SEDIMENT FLUX DUE TO SETTLING
1648  aed_real,INTENT(INOUT) :: y(:)
1649 !
1650 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1651 !
1652 !CONSTANTS
1653  INTEGER,PARAMETER :: ITMAX=100
1654 !
1655  INTEGER I,K,IT
1656  aed_real step_dt,yc,c,cmax
1657  aed_real cu(0:n)
1658 !
1659 !-----------------------------------------------------------------------
1660 !
1661 !BEGIN
1662  fsed = 0.d0 !# INITIALIZE SEDIMENT SETTLING FLUXES WITH ZERO
1663  cu = 0.d0 !# INITIALIZE INTERFACE FLUXES WITH ZERO
1664  cmax = 0.d0 !# INITIALIZE MAXIMUM COURANT NUMBER
1665 
1666 !# COMPUTE MAXIMUM COURANT NUMBER
1667 !# CALCULATED AS NUMBER OF LAYERS THAT THE PARTICLES WILL TRAVEL BASED ON SETTLING OR
1668 !# BUOYANCY VELOCITY
1669 !# THIS NUMBER IS THEN USED TO SPLIT THE VERTICAL MOVEMENT CALCULATIONS TO LIMIT
1670 !# MOVEMENT ACROSS A SINGLE LAYER
1671  DO k=1,n-1
1672 !# SINKING PARTICLES
1673  IF (h_aed2(k+1)+h_aed2(k).GT.eps_aed2) THEN !!ajout mj
1674  c=abs(ww(k+1))*dtaed2/(0.5d0*(h_aed2(k+1)+h_aed2(k)))
1675  IF (c.GT.cmax) cmax=c
1676 !# RISING PARTICLES
1677  c=abs(ww(k))*dtaed2/(0.5d0*(h_aed2(k+1)+h_aed2(k)))
1678  IF (c.GT.cmax) cmax=c
1679  ENDIF
1680  ENDDO
1681 
1682  it=min(itmax,int(cmax)+1)
1683  step_dt = dtaed2 / float(it);
1684 
1685 !# SPLITTING LOOP
1686  DO i=1, it
1687 !# VERTICAL LOOP
1688  DO k=1,n-1
1689 !# COMPUTE THE SLOPE RATION
1690  IF (ww(k).GT.0.d0) THEN !# PARTICLE IS RISING
1691  yc=y(k ) !# CENTRAL VALUE
1692  ELSE !# NEGATIVE SPEED PARTICLE IS SINKING
1693  yc=y(k+1) !# CENTRAL VALUE
1694  ENDIF
1695 !# COMPUTE THE LIMITED FLUX
1696  cu(k)=ww(k) * yc
1697  ENDDO
1698 
1699 !# DO THE UPPER BOUNDARY CONDITIONS
1700  cu(n) = zero_ !# FLUX INTO THE DOMAIN FROM ATMOSPHERE
1701 
1702 !# DO THE LOWER BOUNDARY CONDITIONS
1703  IF (ww(1).GT.0.d0) THEN !# PARTICLE IS RISING
1704  cu(0) = 0.d0 !FLUX FROM BENTHOS IS ZERO
1705  ELSE !# PARTICLE IS SETTLING
1706  cu(0) = ww(1)*y(1)
1707  fsed = cu(0) * step_dt !# FLUX SETTLED INTO THE SEDIMENTS PER SUB TIME STEP
1708  ENDIF
1709 !# DO THE VERTICAL ADVECTION STEP INCLUDING POSITIVE MIGRATION
1710 !# AND SETTLING OF SUSPENDED MATTER.
1711  DO k=1,n
1712  IF (h_aed2(k).GT.eps_aed2) THEN !!ajout mj
1713  y(k)=y(k) - step_dt * ((cu(k) - cu(k-1)) / h_aed2(k))
1714  ENDIF
1715  ENDDO
1716  ENDDO !# END OF THE ITERATION LOOP
1717  fsed = fsed / dtaed2 !# AVERAGE FLUX RATE FOR FULL TIME STEP USED IN AED2
1718 !
1719 !-----------------------------------------------------------------------
1720 !
1721  RETURN
1722  END SUBROUTINE settling
1723 
1724 ! *************************
1725  LOGICAL FUNCTION riparian &
1726 ! *************************
1727  (column, actv, shade_frac, rain_loss)
1729 !***********************************************************************
1730 ! WAQTEL V7P3
1731 !***********************************************************************
1732 !
1733 !brief Do Riparian functionality, including operations in dry and
1734 !+ fringing celles populate feedback arrays to the hoist model
1735 !+ associatied with Riparian effects
1736 !
1737 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1738 !
1739  type(aed2_column_t), INTENT(INOUT) :: column(:)
1740  LOGICAL, INTENT(IN) :: ACTV
1741  aed_real, INTENT(INOUT) :: shade_frac, rain_loss
1742 !
1743 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1744 !
1745  aed_real localshade, localrainl
1746  aed_real :: one = 1.d0
1747 !
1748 !-----------------------------------------------------------------------
1749 !
1750 !BEGIN
1751 !# COMPUTE THE METHODS RELEVANT TO EITHER DRY OR WET CELLS
1752  IF (.NOT. actv) THEN
1753  CALL aed2_calculate_riparian(column, 1, zero_);
1754  CALL aed2_calculate_dry(column, 1);
1755  riparian = .false.
1756  ELSE
1757  CALL aed2_calculate_riparian(column, 1, one);
1758  ENDIF
1759 
1760 !# UPDATE FEEDBACK ARRAYS TO HOST MODEL, TO REDUCE RAIN (OR IF -VE THEN ADD FLOW)
1761 !CALL AED2_RAIN_LOSS(COLUMN, 1, LOCALRAINL); ! COMMENTAIRE MJ
1762  IF (link_rain_loss) rain_loss = localrainl
1763 
1764 !# UPDATE FEEDBACK ARRAYS TO SHADE THE WATER (IE REDUCE INCOMING LIGHT, IO)
1765 !CALL AED2_LIGHT_SHADING(COLUMN, 1, LOCALSHADE)
1766  IF (link_solar_shade) shade_frac = localshade
1767 
1768  riparian = .true.
1769  END FUNCTION riparian
1770 
1771 ! *****************
1772  SUBROUTINE update &
1773 ! *****************
1774  (column,count)
1776 !***********************************************************************
1777 ! WAQTEL V7P3
1778 !***********************************************************************
1779 !
1780 !brief Do non-kinetic (eg equilibrium) updates to state variables in
1781 !+ AED2 modules
1782 !
1783 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1784 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1785 !
1786 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1787 !
1788  type(aed2_column_t), INTENT(INOUT) :: column(:)
1789  INTEGER, INTENT(IN) :: COUNT
1790 !
1791 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1792 !
1793  INTEGER LEV
1794 !
1795 !-----------------------------------------------------------------------
1796 !BEGIN
1797  DO lev=1,count
1798  CALL aed2_equilibrate(column, lev)
1799  ENDDO
1800 !
1801 !-----------------------------------------------------------------------
1802 !
1803  RETURN
1804  END SUBROUTINE update
1805 
1806 ! ********************
1807  SUBROUTINE particles &
1808 ! ********************
1809  (column, count, h_)
1811 !***********************************************************************
1812 ! WAQTEL V7P3
1813 !***********************************************************************
1814 !
1815 !brief Calculate biogeochemical transformations on particles
1816 !+ !to be completed!
1817 !
1818 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1819 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1820 !
1821 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1822 !
1823  type(aed2_column_t), INTENT(INOUT) :: column(:)
1824  INTEGER, INTENT(IN) :: COUNT
1825  aed_real, INTENT(INOUT) :: h_(:)
1826 !
1827 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1828 !
1829  INTEGER I
1830  aed_real zz
1831 !
1832 !-----------------------------------------------------------------------
1833 !
1834 !BEGIN
1835  zz = zero_
1836 
1837 ! CALL AED2_PARTICLE_BGC(COLUMN,1,PPID ...)
1838  IF (count.LE.1) RETURN
1839 
1840  DO i = 2, count
1841 ! CALL AED2_PARTICLE_BGC(COLUMN,COUNT,PPID ...)
1842  ENDDO
1843 !
1844 !-----------------------------------------------------------------------
1845 !
1846  RETURN
1847  END SUBROUTINE particles
1848 
1849 ! ************************
1850  SUBROUTINE bioextinction &
1851 ! ************************
1852  (column,count,extc)
1854 !***********************************************************************
1855 ! WAQTEL V7P3
1856 !***********************************************************************
1857 !
1858 !brief Calculate the specific light attenuation additions due to AED2
1859 !+ modules
1860 !
1861 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1862 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1863 !
1864 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1865 !
1866  type(aed2_column_t), INTENT(INOUT) :: column(:)
1867  INTEGER, INTENT(IN) :: COUNT
1868  aed_real, INTENT(INOUT) :: extc(:)
1869 !
1870 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1871 !
1872  INTEGER I
1873  aed_real localext
1874 !
1875 !-----------------------------------------------------------------------
1876 !
1877 !BEGIN
1878  localext = zero_
1879 
1880  CALL aed2_light_extinction(column, 1, localext)
1881  IF (link_water_clarity) extc(1) = localext
1882  IF (count.LE.1) RETURN
1883 
1884  DO i = 2, count
1885  CALL aed2_light_extinction(column, i, localext)
1886  IF (link_water_clarity) extc(i) = localext
1887  ENDDO
1888 !
1889 !-----------------------------------------------------------------------
1890 !
1891  RETURN
1892  END SUBROUTINE bioextinction
1893 
1894 ! ******************
1895  SUBROUTINE biodrag &
1896 ! ******************
1897  (column,count,bdrag)
1899 !***********************************************************************
1900 ! WAQTEL V7P3
1901 !***********************************************************************
1902 !
1903 !brief Calculate the drag addition to be returned to the host model
1904 !+ due to vegetation
1905 !
1906 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1907 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1908 !
1909 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1910 !
1911  type(aed2_column_t), INTENT(INOUT) :: column(:)
1912  INTEGER, INTENT(IN) :: COUNT
1913  aed_real, INTENT(INOUT) :: bdrag
1914 !
1915 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1916 !
1917  aed_real localdrag
1918 !
1919 !-----------------------------------------------------------------------
1920 !
1921 ! BEGIN
1922 ! CALL AED2_BIO_DRAG(COLUMN, COUNT, LOCALDRAG) ! MODIF MJ DESACTIVE
1923 
1924  IF (link_bottom_drag) bdrag = localdrag
1925 !
1926 !-----------------------------------------------------------------------
1927 !
1928  RETURN
1929  END SUBROUTINE biodrag
1930 
1931 ! *********************
1932  SUBROUTINE biodensity &
1933 ! *********************
1934 !
1935  (column,count,bio_density)
1937 !***********************************************************************
1938 ! WAQTEL V8P0
1939 !***********************************************************************
1940 !
1941 !brief Calculate the density addition to be returned to the host
1942 !+ model due to WQ
1943 !
1944 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1945 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1946 !
1947 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1948 !
1949  type(aed2_column_t), INTENT(INOUT) :: column(:)
1950  INTEGER, INTENT(IN) :: COUNT
1951  aed_real, INTENT(INOUT) :: bio_density(:)
1952 !
1953 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
1954 !
1955 !
1956 !-----------------------------------------------------------------------
1957 !BEGIN
1958 !
1959 !-----------------------------------------------------------------------
1960 !
1961  RETURN
1962  END SUBROUTINE biodensity
1963 #endif
1964 
1965  END MODULE t3d_aed2
logical old_zones
Definition: t3d_aed2.F90:120
subroutine, public init_aed2_models(NAMLST, DNAME, NWQ_VAR, NBEN_VAR, NDIAG_VAR)
Definition: t3d_aed2.F90:134
subroutine biodrag(COLUMN, COUNT, BDRAG)
Definition: t3d_aed2.F90:1898
integer n_vars
Definition: t3d_aed2.F90:124
subroutine light(COLUMN, COUNT, IO, EXTC, PAR_, H_)
Definition: t3d_aed2.F90:1565
integer, dimension(:,:), pointer mat
Definition: t3d_aed2.F90:111
logical link_water_density
Definition: t3d_aed2.F90:114
subroutine bioextinction(COLUMN, COUNT, EXTC)
Definition: t3d_aed2.F90:1853
integer n_vars_diag
Definition: t3d_aed2.F90:124
integer n_vars_ben
Definition: t3d_aed2.F90:124
logical link_rain_loss
Definition: t3d_aed2.F90:114
subroutine calculate_fluxes(COLUMN, COUNT, FLUX_PEL, FLUX_ATM, FLUX_BEN, H_AED2)
Definition: t3d_aed2.F90:1138
subroutine set_initial_from_file
Definition: t3d_aed2.F90:581
subroutine biodensity(COLUMN, COUNT, BIO_DENSITY)
Definition: t3d_aed2.F90:1936
subroutine, public names_var_aed2(NV_AED2, NAMESAED2)
Definition: t3d_aed2.F90:347
logical link_water_clarity
Definition: t3d_aed2.F90:114
logical link_surface_drag
Definition: t3d_aed2.F90:114
integer, dimension(:), pointer benth_map
Definition: t3d_aed2.F90:85
logical link_ext_par
Definition: t3d_aed2.F90:114
subroutine check_data
Definition: t3d_aed2.F90:907
logical function same_str_icase(A, B)
Definition: t3d_aed2.F90:545
logical do_particle_bgc
Definition: t3d_aed2.F90:114
logical, dimension(:), pointer active
Definition: t3d_aed2.F90:84
subroutine, public do_aed2_models(NCELLS, NCOLS, NPLAN, FLUXAED2, EXTCAED2, TA)
Definition: t3d_aed2.F90:1285
logical ext_tss_extinction
Definition: t3d_aed2.F90:114
logical do_zone_averaging
Definition: t3d_aed2.F90:121
integer solution_method
Definition: t3d_aed2.F90:72
logical do_2d_atm_flux
Definition: t3d_aed2.F90:114
character(len=128) init_values_file
Definition: t3d_aed2.F90:81
integer nwq_var
Definition: fv_zones.F90:73
subroutine settling(N, DTAED2, H_AED2, WW, FSED, Y)
Definition: t3d_aed2.F90:1625
logical function riparian(COLUMN, ACTV, SHADE_FRAC, RAIN_LOSS)
Definition: t3d_aed2.F90:1728
character function tolower(C)
Definition: t3d_aed2.F90:516
logical link_solar_shade
Definition: t3d_aed2.F90:114
subroutine define_column(COLUMN, COL, CC, CC_DIAG, FLUX_PEL, FLUX_ATM, FLUX_BEN, NPLAN)
Definition: t3d_aed2.F90:1003
subroutine particles(COLUMN, COUNT, H_)
Definition: t3d_aed2.F90:1810
integer, dimension(:), pointer surf_map
Definition: t3d_aed2.F90:85
integer nben_var
Definition: fv_zones.F90:73
logical link_bottom_drag
Definition: t3d_aed2.F90:114
subroutine check_states(COLUMN, TOP, BOT)
Definition: t3d_aed2.F90:1226
integer n_vars_diag_sheet
Definition: t3d_aed2.F90:125
logical do_limiter
Definition: t3d_aed2.F90:114
integer n_aed2_vars
Definition: t3d_aed2.F90:124
subroutine, public clean_aed2_models
Definition: t3d_aed2.F90:1526
subroutine update(COLUMN, COUNT)
Definition: t3d_aed2.F90:1775
Definition: bief.f:3
subroutine, public set_env_aed2_models(DT_, TEMP_, SALT_, RHO_, H_, TSS_, RAD_, EXTCOEFF_, AREA_, I_0_, WND_, RAIN_, AIR_TEMP_, USTAR_BED_, USTAR_SURF_, Z_, BATHY_, MAT_ID_, ACTIVE_, BENTH_MAP_, BIODRAG_, SOLARSHADE_, RAINLOSS_)
Definition: t3d_aed2.F90:725
subroutine, public init_var_aed2_models(NCELLS, CC_, CC_DIAG_, NWQ, NWQBEN, SM, BM)
Definition: t3d_aed2.F90:394