1!------------------------------------------------------------------------------- 2 3! This file is part of Code_Saturne, a general-purpose CFD tool. 4! 5! Copyright (C) 1998-2021 EDF S.A. 6! 7! This program is free software; you can redistribute it and/or modify it under 8! the terms of the GNU General Public License as published by the Free Software 9! Foundation; either version 2 of the License, or (at your option) any later 10! version. 11! 12! This program is distributed in the hope that it will be useful, but WITHOUT 13! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 15! details. 16! 17! You should have received a copy of the GNU General Public License along with 18! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin 19! Street, Fifth Floor, Boston, MA 02110-1301, USA. 20 21!------------------------------------------------------------------------------- 22!> \file atincl.f90 23!> \brief Module for atmospheric models - main variables 24!>- Nota : ippmod(iatmos) = 0 constante density, =1 --> Dry atmosphere, 25!> = 2 --> Humid atmosphere 26!>- A separate vertical grid is used for 1D radiative scheme 27module atincl 28!> \defgroup at_main 29!============================================================================= 30 31use mesh 32use ppppar 33use ppincl 34 35implicit none 36!> \addtogroup at_main 37!> \{ 38 39!============================================================================= 40 41! 1. Arrays specific to the atmospheric physics 42 43! 1.1 Arrays specific to the input meteo profile 44!----------------------------------------------- 45 46! Arrays specific to values read in the input meteo file: 47 48!> time (in sec) of the meteo profile 49double precision, dimension(:), pointer :: tmmet 50 51!> altitudes of the dynamic profiles (read in the input meteo file) 52double precision, allocatable, dimension(:) :: zdmet 53 54!> Pressure drop integrated over a time step (used for automatic open boundaries) 55double precision, allocatable, dimension(:) :: dpdt_met 56 57!> Momentum for each level (used for automatic open boundaries) 58double precision, allocatable, dimension(:,:) :: mom_met 59double precision, allocatable, dimension(:,:) :: mom 60 61!> altitudes of the temperature profile (read in the input meteo file) 62double precision, dimension(:), pointer :: ztmet 63 64!> meteo u profiles (read in the input meteo file) 65double precision, allocatable, dimension(:,:) :: umet 66 67!> meteo v profiles (read in the input meteo file) 68double precision, allocatable, dimension(:,:) :: vmet 69 70!> meteo w profiles - unused 71double precision, allocatable, dimension(:,:) :: wmet 72 73!> meteo turbulent kinetic energy profile (read in the input meteo file) 74double precision, allocatable, dimension(:,:) :: ekmet 75 76!> meteo turbulent dissipation profile (read in the input meteo file) 77double precision, allocatable, dimension(:,:) :: epmet 78 79!> meteo temperature (Celsius) profile (read in the input meteo file) 80double precision, allocatable, dimension(:,:) :: ttmet 81 82!> meteo specific humidity profile (read in the input meteo file) 83double precision, allocatable, dimension(:,:) :: qvmet 84 85!> meteo specific droplet number profile (read in the input meteo file) 86double precision, allocatable, dimension(:,:) :: ncmet 87 88!> Sea level pressure (read in the input meteo file) 89double precision, allocatable, dimension(:) :: pmer 90 91!> X axis coordinates of the meteo profile (read in the input meteo file) 92double precision, allocatable, dimension(:) :: xmet 93 94!> Y axis coordinates of the meteo profile (read in the input meteo file) 95double precision, allocatable, dimension(:) :: ymet 96 97! Arrays specific to values calculated from the meteo file (cf atlecm.f90): 98 99!> density profile 100double precision, allocatable, dimension(:,:) :: rmet 101 102!> potential temperature profile 103double precision, allocatable, dimension(:,:) :: tpmet 104 105!> hydrostatic pressure from Laplace integration 106double precision, dimension(:,:), pointer :: phmet 107 108! 1.2 Pointers for the positions of the variables 109!------------------------------------------------ 110! Variables specific to the atmospheric physics: 111!> total water content (for humid atmosphere) 112integer, save :: iymw 113!> intdrp---> total number of droplets (for humid atmosphere) 114integer, save :: intdrp = -1 115 116! 1.3 Pointers for the positions of the properties for the specific phys. 117!------------------------------------------------------------------------ 118! Properties specific to the atmospheric physics: 119 120!> temperature (in Celsius) 121integer, save :: itempc 122 123!> liquid water content 124integer, save :: iliqwt 125 126!> momentum source term field id (useful when iatmst > 0) 127integer, save :: imomst 128 129!---------------------------------------------------------------------------- 130 131! 2. Data specific to the atmospheric physics 132 133! 2.1 Data specific to the input meteo profile 134!---------------------------------------------- 135 136!>flag for reading the meteo input file 137!> - = 0 -> no reading 138!> - = 1 -> reading 139integer(c_int), pointer, save :: imeteo 140 141!> numbers of altitudes for the dynamics 142integer(c_int), pointer, save :: nbmetd 143 144!> numbers of altitudes for the temperature and specific humidity 145integer(c_int), pointer, save :: nbmett 146 147!> numbers of time steps for the meteo profiles 148integer(c_int), pointer, save :: nbmetm 149 150!> read zone boundary conditions from profile 151integer, save :: iprofm(nozppm) 152 153!> automatic inlet/outlet boundary condition flag 154!> (0: not auto (default); 1,2: auto) 155!> When meteo momentum source terms are activated (iatmst > 0), 156!> iautom = 1 corresponds to a Dirichlet on the pressure and a 157!> Neumann on the velocity, whereas iautom = 2 imposes a Dirichlet 158!> on both pressure and velocity 159integer, allocatable, target, dimension(:) :: iautom 160 161!> use meteo profile for variables initialization 162!> (0: not used; 1: used (default)) 163integer, save :: initmeteo 164 165!> add a momentum source term based on the meteo profile 166!> for automatic open boundaries 167integer(c_int), pointer, save :: iatmst 168 169!> flag for meteo velocity field interpolation 170!> - 0: linear interpolation of the meteo profile 171!> - 1: the user can directly impose the exact meteo velocity 172!> by declaring the 'meteo_velocity' field 173!> Useful for iatmst = 1 174!> Note: deprecated, imeteo=2 can be used instead. 175integer, save :: theo_interp 176 177! 2.1 Constant specific to the physics (defined in atini1.f90) 178!------------------------------------------------------------------------------- 179 180!> reference pressure (to compute potential temp: 1.0d+5) 181real(c_double), pointer, save:: ps 182 183! 2.2. Space and time reference of the run 184!------------------------------------------------------------------------------- 185 186!> starting year 187integer(c_int), pointer, save:: syear 188 189!> starting quantile 190integer(c_int), pointer, save:: squant 191 192!> starting hour 193integer(c_int), pointer, save:: shour 194 195!> starting min 196integer(c_int), pointer, save:: smin 197 198!> starting second 199real(c_double), pointer, save :: ssec 200 201!> longitude of the domain origin 202real(c_double), pointer, save:: xlon 203 204!> latitude of the domain origin 205real(c_double), pointer, save:: xlat 206 207!> x coordinate of the domain origin in Lambert-93 208real(c_double), pointer, save:: xl93 209 210!> y coordinate of the domain origin in Lambert-93 211real(c_double), pointer, save:: yl93 212 213! 2.3 Data specific to the meteo profile above the domain 214!-------------------------------------------------------- 215!> Number of vertical levels (cf. 1D radiative scheme 216integer(c_int), pointer, save:: nbmaxt 217 218!> flag to compute the hydrostastic pressure by Laplace integration 219!> in the meteo profiles 220!> = 0 : bottom to top Laplace integration, based on P(sea level) (default) 221!> = 1 : top to bottom Laplace integration based on P computed for 222!> the standard atmosphere at z(nbmaxt) 223integer, save:: ihpm 224 225! 2.4 Data specific to the 1D vertical grid: 226!------------------------------------------- 227 228!> flag for the definition of the vertical grid 229!> - -1 : no vertical grid (default) 230!> - 0 : automatic definition 231!> - 1 : definition by the user in cs_user_atmospheric_model.f90 232integer, save:: ivert 233 234!> number of vertical arrays 235integer, save:: nvert 236 237!> number of levels (up to the top of the domain) 238integer, save:: kvert 239 240!> Number of levels (up to 11000 m if ray1d used) 241!> (automatically computed) 242integer, save:: kmx 243 244!> Height of the boundary layer 245real(c_double), pointer, save :: meteo_zi 246 247! 2.5 Data specific to the 1d atmospheric radiative module: 248!------------------------------------------------------------------------------- 249!> flag for the use of the 1d atmo radiative model 250!> - 0 no use (default) 251!> - 1 use 252integer, save:: iatra1 253 254!> 1d radiative model pass frequency 255integer, save:: nfatr1 256 257!> flag for the standard atmo humidity profile 258!> - 0: q = 0 (default) 259!> - 1: q = decreasing exponential 260integer, save:: iqv0 261 262!> pointer for 1D infrared profile 263integer, save:: idrayi 264 265!> pointer for 1D solar profile 266integer, save:: idrayst 267 268!> grid formed by 1D profiles 269integer, save:: igrid 270 271!> Flag for the computation of downward and upward infrared radiative fluxes 272!> 0: disabled 273!> 1: enabled 274integer, save:: irdu 275 276! 2.6 Arrays specific to the 1d atmospheric radiative module 277!------------------------------------------------------------------------------- 278 279!> horizontal coordinates of the vertical grid 280double precision, allocatable, dimension(:,:) :: xyvert 281 282!> vertical grid for 1D radiative scheme initialize in 283!> cs_user_atmospheric_model.f90 284double precision, allocatable, dimension(:) :: zvert 285 286!> absorption for CO2 + 03 287double precision, allocatable, dimension(:) :: acinfe 288 289!> differential absorption for CO2 + 03 290double precision, allocatable, dimension(:) :: dacinfe 291 292!> absorption for CO2 only 293double precision, allocatable, dimension(:,:) :: aco2, aco2s 294 295!> differential absorption for CO2 only 296double precision, allocatable, dimension(:,:) :: daco2, daco2s 297 298!> idem acinfe, flux descendant 299double precision, allocatable, dimension(:) :: acsup, acsups 300 301!> internal variable for 1D radiative model 302double precision, allocatable, dimension(:) :: dacsup, dacsups 303 304!> internal variable for 1D radiative model 305double precision, allocatable, dimension(:) :: tauzq 306 307!> internal variable for 1D radiative model 308double precision, allocatable, dimension(:) :: tauz 309 310!> internal variable for 1D radiative model 311double precision, allocatable, dimension(:) :: zq 312 313!> internal variable for 1D radiative model 314double precision, save :: tausup 315 316!> internal variable for 1D radiative model 317double precision, allocatable, dimension(:) :: zray 318double precision, allocatable, dimension(:,:) :: rayi, rayst 319 320!> Upward and downward radiative fluxes (infrared, solar) along each vertical 321double precision, allocatable, dimension(:,:) :: iru, ird, solu, sold 322 323! 3.0 Data specific to the ground model 324!------------------------------------------------------------------------------- 325!> iatsoil --> flag to use the ground model 326integer, save:: iatsoil 327 328!> Water content of the first ground reservoir 329double precision, save:: w1ini 330 331!> Water content of the second ground reservoir 332double precision, save:: w2ini 333 334!> Do we compute z ground every where? 335logical(c_bool), pointer, save :: compute_z_ground 336 337! ------------------------------------------------------------------------------- 338! 4.0 Microphysics parameterization options 339! ------------------------------------------------------------------------------- 340 341!> Option for subgrid models 342!> - modsub = 0 : the simplest parameterization (for numerical verifications) 343!> - modsub = 1 : Bechtold et al. 1995 (Luc Musson-Genon) 344!> - modsub = 2 : Bouzereau et al. 2004 345!> - modsub = 3 : Cuijpers and Duynkerke 1993, Deardorff 1976, Sommeria and 346!> Deardorff 1977 347integer(c_int), pointer, save:: modsub 348 349!> Option for liquid water content distribution models 350!> - moddis = 1 : all or nothing 351!> - moddis = 2 : Gaussian distribution 352integer(c_int), pointer, save:: moddis 353 354!> Option for nucleation 355!> - modnuc = 0 : without nucleation 356!> - modnuc = 1 : Pruppacher and Klett 1997 357!> - modnuc = 2 : Cohard et al. 1998,1999 358!> - modnuc = 3 : Abdul-Razzak et al. 1998,2000 359!> logaritmic standard deviation of the log-normal law of the droplet spectrum 360integer(c_int), pointer, save:: modnuc 361 362!> sedimentation flag 363integer(c_int), pointer, save:: modsedi 364 365!> deposition flag 366integer(c_int), pointer, save:: moddep 367 368!> adimensional : sigc=0.53 other referenced values are 0.28, 0.15 369double precision, save:: sigc 370 371!> force initilization in case of restart (this option is 372!> automatically set in lecamp) 373integer, save :: init_at_chem 374 375!> key id for optimal interpolation 376integer, save :: kopint 377 378!> Aerosol optical properties 379 380! Aerosol optical depth 381!> adimensional : aod_o3_tot=0.2 other referenced values are 0.10, 0.16 382double precision, save:: aod_o3_tot 383!> adimensional : aod_h2o_tot=0.10 other referenced values are 0.06, 0.08 384double precision, save:: aod_h2o_tot 385 386!> Asymmetry factor for O3 (non-dimensional) 387!> climatic value gaero_o3=0.66 388double precision, save:: gaero_o3 389!> Asymmetry factor for H2O (non-dimensional) 390!> climatic value gaero_h2o=0.64 391double precision, save:: gaero_h2o 392 393!> Single scattering albedo for O3 (non-dimensional) 394!> climatic value piaero_o3=0.84, other referenced values are 0.963 395double precision, save:: piaero_o3 396!> Single scattering albedo for H2O (non-dimensional) 397!> climatic value piaero_h2o=0.84, other referenced values are 0.964 398double precision, save:: piaero_h2o 399 400!> Fraction of Black carbon (non-dimensional): black_carbon_frac=1.d-8 for no BC 401double precision, save:: black_carbon_frac 402 403!> Maximal height for aerosol distribution on the vertical 404!> important should be <= zqq(kmray-1); 405!> in meters : referenced value: zaero=6000 406double precision, save:: zaero 407!> \} 408 409!============================================================================= 410 411 interface 412 413 !--------------------------------------------------------------------------- 414 415 ! Interface to C function returning a meteo file name 416 417 subroutine cs_f_atmo_get_meteo_file_name(f_name_max, f_name, f_name_len) & 418 bind(C, name='cs_f_atmo_get_meteo_file_name') 419 use, intrinsic :: iso_c_binding 420 implicit none 421 integer(c_int), value :: f_name_max 422 type(c_ptr), intent(out) :: f_name 423 integer(c_int), intent(out) :: f_name_len 424 end subroutine cs_f_atmo_get_meteo_file_name 425 426 !--------------------------------------------------------------------------- 427 428 !> \brief Return pointers to atmo includes 429 430 subroutine cs_f_atmo_get_pointers(ps, & 431 syear, squant, shour, smin, ssec, & 432 longitude, latitude, & 433 x_l93, y_l93, & 434 compute_z_ground, iatmst, & 435 sedimentation_model, deposition_model, nucleation_model, & 436 subgrid_model, distribution_model, & 437 ichemistry, nespg, nrg, chem_with_photo, & 438 iaerosol, frozen_gas_chem, init_gas_with_lib, & 439 init_aero_with_lib, n_aero, n_sizebin, imeteo, & 440 nbmetd, nbmett, nbmetm, nbmaxt, & 441 meteo_zi) & 442 bind(C, name='cs_f_atmo_get_pointers') 443 use, intrinsic :: iso_c_binding 444 implicit none 445 type(c_ptr), intent(out) :: ps 446 type(c_ptr), intent(out) :: compute_z_ground, iatmst 447 type(c_ptr), intent(out) :: ichemistry, nespg, nrg 448 type(c_ptr), intent(out) :: sedimentation_model, deposition_model 449 type(c_ptr), intent(out) :: nucleation_model 450 type(c_ptr), intent(out) :: subgrid_model, distribution_model 451 type(c_ptr), intent(out) :: syear, squant, shour, smin, ssec 452 type(c_ptr), intent(out) :: longitude, latitude 453 type(c_ptr), intent(out) :: x_l93, y_l93 454 type(c_ptr), intent(out) :: iaerosol, frozen_gas_chem 455 type(c_ptr), intent(out) :: init_gas_with_lib, init_aero_with_lib 456 type(c_ptr), intent(out) :: n_aero, n_sizebin, chem_with_photo 457 type(c_ptr), intent(out) :: imeteo 458 type(c_ptr), intent(out) :: nbmetd, nbmett, nbmetm, nbmaxt 459 type(c_ptr), intent(out) :: meteo_zi 460 end subroutine cs_f_atmo_get_pointers 461 462 !--------------------------------------------------------------------------- 463 464 !> \brief Initialize meteo profiles if no meteo file is given 465 466 subroutine cs_atmo_init_meteo_profiles() & 467 bind(C, name='cs_atmo_init_meteo_profiles') 468 use, intrinsic :: iso_c_binding 469 implicit none 470 end subroutine cs_atmo_init_meteo_profiles 471 472 !--------------------------------------------------------------------------- 473 474 !> \brief Compute meteo profiles if no meteo file is given 475 476 subroutine cs_atmo_compute_meteo_profiles() & 477 bind(C, name='cs_atmo_compute_meteo_profiles') 478 use, intrinsic :: iso_c_binding 479 implicit none 480 end subroutine cs_atmo_compute_meteo_profiles 481 482 !--------------------------------------------------------------------------- 483 484 !> \brief Calculation of the specific enthalpy of liquid water 485 486 !> \return specific enthalpy of liquid water 487 488 !> \param[in] t_l liquid water temperature (in Celsius) 489 490 function cs_liq_t_to_h(t_l) result(h_l) & 491 bind(C, name='cs_liq_t_to_h') 492 use, intrinsic :: iso_c_binding 493 implicit none 494 real(c_double), value :: t_l 495 real(c_double) :: h_l 496 end function cs_liq_t_to_h 497 498 !--------------------------------------------------------------------------- 499 500 !> \brief Calculation of the absolute humidity at saturation 501 !> for a given temperature. 502 503 !> \param[in] t_c temperature (in Celsius) 504 !> \param[in] p pressure 505 506 function cs_air_x_sat(t_c, p) result(x_s) & 507 bind(C, name='cs_air_x_sat') 508 use, intrinsic :: iso_c_binding 509 implicit none 510 real(c_double), value :: t_c, p 511 real(c_double) :: x_s 512 end function cs_air_x_sat 513 514 !--------------------------------------------------------------------------- 515 516 !> \brief Calculation of the air water mass fraction at saturation 517 !> for a given temperature. 518 519 !> \param[in] t_c temperature (in Celsius) 520 !> \param[in] p pressure 521 522 function cs_air_yw_sat(t_c, p) result(x_s) & 523 bind(C, name='cs_air_yw_sat') 524 use, intrinsic :: iso_c_binding 525 implicit none 526 real(c_double), value :: t_c, p 527 real(c_double) :: x_s 528 end function cs_air_yw_sat 529 530 !--------------------------------------------------------------------------- 531 532 !> \brief Computes the saturation water vapor pressure function 533 !> of the temperature (C). 534 535 !> \param[in] t_c temperature (in Celsius) 536 537 function cs_air_pwv_sat(t_c) result(x_s) & 538 bind(C, name='cs_air_pwv_sat') 539 use, intrinsic :: iso_c_binding 540 implicit none 541 real(c_double), value :: t_c 542 real(c_double) :: x_s 543 end function cs_air_pwv_sat 544 545 !--------------------------------------------------------------------------- 546 547 !> \brief Convert the absolute humidity of humid air to the air water mass fraction. 548 549 !> \param[in] x absolute humidity of humid air 550 551 function cs_air_x_to_yw(x) result(qw) & 552 bind(C, name='cs_air_x_to_yw') 553 use, intrinsic :: iso_c_binding 554 implicit none 555 real(c_double), value :: x 556 real(c_double) :: qw 557 end function cs_air_x_to_yw 558 559 !--------------------------------------------------------------------------- 560 561 !> \brief Convert the air water mass fraction to the absolute humidity of humid air. 562 563 !> \param[in] qw air water mass fraction 564 565 function cs_air_yw_to_x(qw) result(x) & 566 bind(C, name='cs_air_yw_to_x') 567 use, intrinsic :: iso_c_binding 568 implicit none 569 real(c_double), value :: qw 570 real(c_double) :: x 571 end function cs_air_yw_to_x 572 573 !--------------------------------------------------------------------------- 574 575 !> \brief Calculation of the density of humid air. 576 577 !> \param[in] ywm air water mass fraction 578 !> \param[in] t_liq pressure 579 !> \param[in] p pressure 580 !> \param[out] yw_liq liquid water mass fraction 581 !> \param[out] t_h temperature of humid air in Celsius 582 !> \param[out] rho_h density of humid air 583 584 subroutine cs_rho_humidair(ywm, t_liq, p, yw_liq, t_h, rho_h) & 585 bind(C, name='cs_rho_humidair') 586 use, intrinsic :: iso_c_binding 587 implicit none 588 real(c_double), value :: ywm, t_liq, p 589 real(c_double), intent(out) :: yw_liq, t_h, rho_h 590 end subroutine cs_rho_humidair 591 592 !============================================================================= 593 594 end interface 595 596contains 597 598 !============================================================================= 599 600 !> \brief Return meteo file name 601 602 !> \param[out] name meteo file name 603 604 subroutine atmo_get_meteo_file_name(name) 605 606 use, intrinsic :: iso_c_binding 607 implicit none 608 609 ! Arguments 610 611 character(len=*), intent(out) :: name 612 613 ! Local variables 614 615 integer :: i 616 integer(c_int) :: name_max, c_name_len 617 type(c_ptr) :: c_name_p 618 character(kind=c_char, len=1), dimension(:), pointer :: c_name 619 620 name_max = len(name) 621 622 call cs_f_atmo_get_meteo_file_name(name_max, c_name_p, c_name_len) 623 call c_f_pointer(c_name_p, c_name, [c_name_len]) 624 625 do i = 1, c_name_len 626 name(i:i) = c_name(i) 627 enddo 628 do i = c_name_len + 1, name_max 629 name(i:i) = ' ' 630 enddo 631 632 return 633 634 end subroutine atmo_get_meteo_file_name 635 636 !============================================================================= 637 638 !> \brief Return pointer to automatic face bc flag array 639 640 !> \return auto_flag pointer to automatic boundary condition array 641 642 function cs_atmo_get_auto_flag() result(auto_flag) & 643 bind(C, name='cs_atmo_get_auto_flag') 644 use, intrinsic :: iso_c_binding 645 implicit none 646 type(c_ptr) :: auto_flag 647 auto_flag = c_loc(iautom) 648 end function cs_atmo_get_auto_flag 649 650 !============================================================================= 651 652 !> \brief Map fortran to C variables 653 654 subroutine atmo_init 655 656 use, intrinsic :: iso_c_binding 657 use atchem, only: nrg, nespg, ichemistry, photolysis 658 use sshaerosol 659 use cs_c_bindings 660 661 implicit none 662 663 ! Local variables 664 type(c_ptr) :: c_ps 665 type(c_ptr) :: c_compute_z_ground, c_iatmst, c_model, c_nrg, c_nespg 666 type(c_ptr) :: c_sedimentation_model, c_deposition_model, c_nucleation_model 667 type(c_ptr) :: c_subgrid_model 668 type(c_ptr) :: c_distribution_model 669 type(c_ptr) :: c_syear, c_squant, c_shour, c_smin, c_ssec 670 type(c_ptr) :: c_longitude, c_latitude 671 type(c_ptr) :: c_xl93, c_yl93 672 type(c_ptr) :: c_modelaero, c_frozen_gas_chem, c_nlayer, c_nsize 673 type(c_ptr) :: c_init_gas_with_lib, c_init_aero_with_lib, c_chem_with_photo 674 type(c_ptr) :: c_imeteo 675 type(c_ptr) :: c_nbmetd, c_nbmett, c_nbmetm, c_nbmaxt 676 type(c_ptr) :: c_meteo_zi 677 678 call cs_f_atmo_get_pointers(c_ps, & 679 c_syear, c_squant, c_shour, c_smin, c_ssec, & 680 c_longitude, c_latitude, & 681 c_xl93, c_yl93, & 682 c_compute_z_ground, c_iatmst, & 683 c_sedimentation_model, c_deposition_model, & 684 c_nucleation_model, c_subgrid_model, & 685 c_distribution_model, & 686 c_model, c_nespg, c_nrg, c_chem_with_photo, & 687 c_modelaero, c_frozen_gas_chem, & 688 c_init_gas_with_lib, & 689 c_init_aero_with_lib, c_nlayer, & 690 c_nsize, c_imeteo, & 691 c_nbmetd, c_nbmett, c_nbmetm, c_nbmaxt, & 692 c_meteo_zi) 693 694 call c_f_pointer(c_ps, ps) 695 call c_f_pointer(c_syear, syear) 696 call c_f_pointer(c_squant, squant) 697 call c_f_pointer(c_shour, shour) 698 call c_f_pointer(c_smin, smin) 699 call c_f_pointer(c_ssec, ssec) 700 701 call c_f_pointer(c_longitude, xlon) 702 call c_f_pointer(c_latitude, xlat) 703 call c_f_pointer(c_xl93, xl93) 704 call c_f_pointer(c_yl93, yl93) 705 706 call c_f_pointer(c_compute_z_ground, compute_z_ground) 707 call c_f_pointer(c_iatmst, iatmst) 708 709 call c_f_pointer(c_sedimentation_model, modsedi) 710 call c_f_pointer(c_deposition_model, moddep) 711 call c_f_pointer(c_nucleation_model, modnuc) 712 call c_f_pointer(c_subgrid_model, modsub) 713 call c_f_pointer(c_distribution_model, moddis) 714 715 call c_f_pointer(c_model, ichemistry) 716 call c_f_pointer(c_nespg, nespg) 717 call c_f_pointer(c_nrg, nrg) 718 call c_f_pointer(c_chem_with_photo, photolysis) 719 call c_f_pointer(c_modelaero, iaerosol) 720 call c_f_pointer(c_frozen_gas_chem, nogaseouschemistry) 721 call c_f_pointer(c_init_gas_with_lib, init_gas_with_lib) 722 call c_f_pointer(c_init_aero_with_lib, init_aero_with_lib) 723 call c_f_pointer(c_nlayer, nlayer_aer) 724 call c_f_pointer(c_nsize, n_aer) 725 call c_f_pointer(c_imeteo, imeteo) 726 call c_f_pointer(c_nbmetd, nbmetd) 727 call c_f_pointer(c_nbmett, nbmett) 728 call c_f_pointer(c_nbmetm, nbmetm) 729 call c_f_pointer(c_nbmaxt, nbmaxt) 730 call c_f_pointer(c_meteo_zi, meteo_zi) 731 732 return 733 734 end subroutine atmo_init 735 736!=============================================================================== 737 738!> \brief Initialisation of meteo data 739subroutine init_meteo 740 741use cs_c_bindings 742use atsoil 743 744implicit none 745 746! Local variables 747integer :: imode, n_level, n_times, n_level_t 748 749type(c_ptr) :: c_z_temp_met, c_time_met 750type(c_ptr) :: c_hyd_p_met 751 752integer(c_int), dimension(2) :: dim_hyd_p_met 753 754if (imeteo.eq.1) then 755 call atlecm(0) 756endif 757if (imeteo.eq.2) then 758 call cs_atmo_init_meteo_profiles() 759endif 760 761call cs_f_atmo_arrays_get_pointers(c_z_temp_met, c_time_met, & 762 c_hyd_p_met, dim_hyd_p_met) 763 764call c_f_pointer(c_z_temp_met, ztmet, [nbmaxt]) 765call c_f_pointer(c_time_met, tmmet, [nbmetm]) 766call c_f_pointer(c_hyd_p_met, phmet, [dim_hyd_p_met]) 767 768! Allocate additional arrays for Water Microphysics 769 770if (imeteo.gt.0) then 771 772 ! NB : only ztmet,ttmet,qvmet,ncmet are extended to 11000m if iatr1=1 773 ! rmet,tpmet,phmet 774 n_level = max(1, nbmetd) 775 n_times = max(1, nbmetm) 776 n_level_t = max(1, nbmaxt) 777 allocate(zdmet(n_level)) 778 779 if (iatmst.ge.1) then 780 allocate(dpdt_met(n_level)) 781 allocate(mom(3, n_level)) 782 allocate(mom_met(3, n_level)) 783 endif 784 785 allocate(umet(n_level,n_times), vmet(n_level,n_times), wmet(n_level,n_times)) 786 allocate(ekmet(n_level,n_times), epmet(n_level,n_times)) 787 allocate(ttmet(n_level_t,n_times), qvmet(n_level_t,n_times), ncmet(n_level_t,n_times)) 788 allocate(pmer(n_times)) 789 allocate(xmet(n_times), ymet(n_times)) 790 allocate(rmet(n_level_t,n_times), tpmet(n_level_t,n_times)) 791 792 ! Allocate additional arrays for 1D radiative model 793 794 if (iatra1.eq.1) then 795 796 imode = 0 797 798 call usatdv(imode) 799 800 allocate(xyvert(nvert,3), zvert(kmx)) 801 allocate(acinfe(kmx), dacinfe(kmx), aco2(kmx,kmx), aco2s(kmx,kmx)) 802 allocate(daco2(kmx,kmx), daco2s(kmx,kmx), acsup(kmx), dacsup(kmx)) 803 allocate(acsups(kmx), dacsups(kmx)) 804 allocate(tauzq(kmx+1), tauz(kmx+1), zq(kmx+1)) 805 allocate(rayi(kmx,nvert), rayst(kmx,nvert), zray(kmx)) 806 807 allocate(soilvert(nvert)) 808 809 call mestcr("rayi", len("rayi"), 1, 0, idrayi) 810 call mestcr("rayst", len("rayst"), 1, 0, idrayst) 811 call gridcr("int_grid", len("int_grid"), igrid) 812 813 if (irdu.eq.1) then 814 allocate(iru(kmx,nvert), ird(kmx,nvert)) 815 endif 816 817 allocate(sold(kmx,nvert), solu(kmx,nvert)) 818 819 endif 820 821endif 822 823end subroutine init_meteo 824 825!=============================================================================== 826 827!> \brief Initialisation of meteo data 828subroutine init_atmo_autom(nfabor) 829 830use cs_c_bindings 831 832implicit none 833 834integer nfabor 835 836! Local variables 837integer :: ifac 838 839! Allocate additional arrays for Water Microphysics 840if (imeteo.gt.0) then 841 842 ! Allocate and initialize auto inlet/outlet flag 843 allocate(iautom(nfabor)) 844 do ifac = 1, nfabor 845 iautom(ifac) = 0 846 enddo 847 848endif 849 850end subroutine init_atmo_autom 851 852!============================================================================= 853 854!> \brief Final step for deallocation 855subroutine finalize_meteo 856 857use cs_c_bindings 858use atsoil 859 860implicit none 861 862if (imeteo.gt.0) then 863 864 deallocate(zdmet) 865 if (allocated(mom)) then 866 deallocate(mom, mom_met, dpdt_met) 867 endif 868 deallocate(umet, vmet, wmet) 869 deallocate(ekmet, epmet) 870 deallocate(ttmet, qvmet, ncmet) 871 deallocate(pmer) 872 deallocate(xmet, ymet) 873 deallocate(rmet, tpmet) 874 875 deallocate(iautom) 876 877 if (iatra1.eq.1) then 878 879 deallocate(xyvert, zvert) 880 deallocate(acinfe, dacinfe, aco2, aco2s) 881 deallocate(daco2, daco2s, acsup, dacsup) 882 deallocate(acsups, dacsups) 883 deallocate(tauzq, tauz, zq) 884 deallocate(rayi, rayst, zray) 885 886 deallocate(soilvert) 887 888 call mestde () 889 890 call grides () 891 892 if (irdu.eq.1) then 893 deallocate(iru, ird) 894 endif 895 896 deallocate(solu, sold) 897 898 endif 899 900endif 901 902end subroutine finalize_meteo 903 904!--------------------------------------------------------------------------- 905 906!> \brief Universal functions of Cheng and Brutsaert 2005, for stable 907!> (derivative function) 908 909!> \param[in] z altitude 910!> \param[in] dlmo inverse Monin Obukhov length 911!> \param[out] coef function 912 913subroutine mo_phim_s (z,dlmo,coef) 914 915 implicit none 916 917 double precision z,dlmo,coef 918 919 double precision a,b,x 920 921 a=6.1d0 922 b=2.5d0 923 x=z * dlmo 924 925 coef=1.d0+a*(x+(x**b)*((1.d0+x**b)**((1.d0-b)/b)))/(x+(1.d0+x**b)**(1.d0/b)) 926 927end subroutine mo_phim_s 928 929subroutine mo_phih_s (z,dlmo,coef) 930 931 implicit none 932 933 double precision z,dlmo,coef 934 935 double precision a,b,x 936 937 a=5.3d0 938 b=1.1d0 939 x=z * dlmo 940 941 coef=1.d0+a*(x+(x**b)*((1.d0+x**b)**((1.d0-b)/b)))/(x+(1.d0+x**b)**(1.d0/b)) 942 943end subroutine mo_phih_s 944 945! Integrated version from z0 to z 946 947subroutine mo_psim_s (z,z0,dlmo,coef) 948 949 implicit none 950 951 double precision z,z0,dlmo,coef 952 953 double precision a,b,x,x0 954 955 a=6.1d0 956 b=2.5d0 957 x=z * dlmo 958 x0=z0 * dlmo 959 960 coef=dlog(z/z0)+a*(dlog(x+(1.d0+x**b)**(1.d0/b))-dlog(x0+(1.d0+x0**b)**(1.d0/b))) 961 962end subroutine mo_psim_s 963 964subroutine mo_psih_s (z,z0,dlmo,coef) 965 966 implicit none 967 968 double precision z,z0,dlmo,coef 969 970 double precision a,b,x,x0 971 972 a=5.3d0 973 b=1.1d0 974 x=z * dlmo 975 x0=z0 * dlmo 976 977 coef=dlog(z/z0)+a*(dlog(x+(1.d0+x**b)**(1.d0/b))-dlog(x0+(1.d0+x0**b)**(1.d0/b))) 978 979end subroutine mo_psih_s 980 981!--------------------------------------------------------------------------- 982 983!> \brief Universal functions of Hogstrom 1988, for unstable 984!> (derivative function) 985 986!> \param[in] z altitude 987!> \param[in] dlmo Monin Obukhov length 988!> \param[out] coef function 989 990subroutine mo_phim_u (z,dlmo,coef) 991 992 implicit none 993 994 double precision z,dlmo,coef 995 996 double precision a,b,e,x 997 998 a=1.d0 999 b=19.3d0 1000 e=-0.25d0 1001 x=z * dlmo 1002 1003 coef=a*(1.d0-b*x)**e 1004 1005end subroutine mo_phim_u 1006 1007subroutine mo_phih_u (z,dlmo,coef) 1008 1009 implicit none 1010 1011 double precision z,dlmo,coef 1012 1013 double precision a,b,e,x 1014 1015 a=0.95d0 1016 b=11.6d0 1017 e=-0.5d0 1018 x=z * dlmo 1019 1020 coef=a*(1.d0-b*x)**e 1021 1022end subroutine mo_phih_u 1023 1024! Integrated version from z0 to z 1025 1026subroutine mo_psim_u (z,z0,dlmo,coef) 1027 1028 implicit none 1029 1030 double precision z,z0,dlmo,coef 1031 1032 double precision a,b,e,x,x0,psi,psi0 1033 1034 a=1.d0 1035 b=19.3d0 1036 e=0.25d0 1037 x=z * dlmo 1038 x0=z0 * dlmo 1039 psi=(1.d0-b*x)**e 1040 psi0=(1.d0-b*x0)**e 1041 1042 coef=a*(dlog(z/z0) & 1043 -2.d0*dlog((1.d0+psi)/(1.d0+psi0)) & 1044 -dlog((1.d0+psi**2.d0)/(1.d0+psi0**2.d0)) & 1045 +2.d0*(datan(psi)-datan(psi0))) 1046 1047end subroutine mo_psim_u 1048 1049subroutine mo_psih_u (z,z0,dlmo,coef) 1050 1051 implicit none 1052 1053 double precision z,z0,dlmo,coef 1054 1055 double precision a,b,e,x,x0,psi,psi0 1056 1057 a=0.95d0 1058 b=11.6d0 1059 e=0.5d0 1060 x=z * dlmo 1061 x0=z0 * dlmo 1062 psi=(1.d0-b*x)**e 1063 psi0=(1.d0-b*x0)**e 1064 1065 coef=a*(dlog(z/z0) & 1066 -2.d0*dlog((1.d0+psi)/(1.d0+psi0))) 1067end subroutine mo_psih_u 1068 1069!--------------------------------------------------------------------------- 1070 1071!> \brief Universal functions, for neutral 1072!> (derivative function) 1073 1074!> \param[in] z altitude 1075!> \param[in] dlmo Inverse Monin Obukhov length 1076!> \param[out] coef function 1077 1078subroutine mo_phim_n (z,dlmo,coef) 1079 1080 implicit none 1081 1082 double precision z,dlmo,coef 1083 1084 coef=1.d0 1085 1086end subroutine mo_phim_n 1087 1088subroutine mo_phih_n (z,dlmo,coef) 1089 1090 implicit none 1091 1092 double precision z,dlmo,coef 1093 1094 coef=1.d0 1095 1096end subroutine mo_phih_n 1097 1098! Integrated version from z0 to z 1099 1100subroutine mo_psim_n (z,z0,dlmo,coef) 1101 1102 implicit none 1103 1104 double precision z,z0,dlmo,coef 1105 1106 coef=dlog(z/z0) 1107 1108end subroutine mo_psim_n 1109 1110subroutine mo_psih_n (z,z0,dlmo,coef) 1111 1112 implicit none 1113 1114 double precision z,z0,dlmo,coef 1115 1116 coef=dlog(z/z0) 1117 1118end subroutine mo_psih_n 1119 1120! Switch universal functions 1121 1122! Derivative function 1123 1124function cs_mo_phim(z,dlmo) result(coef) & 1125 bind(C, name='cs_mo_phim') 1126 use, intrinsic :: iso_c_binding 1127 1128 implicit none 1129 1130 real(c_double), value :: z,dlmo 1131 real(c_double) :: coef 1132 double precision :: dlmoneutral = 1.d-12 1133 1134 if (abs(dlmo).lt.dlmoneutral) then 1135 call mo_phim_n(z,dlmo,coef) 1136 elseif (dlmo.ge.0.d0) then 1137 call mo_phim_s(z,dlmo,coef) 1138 else 1139 call mo_phim_u(z,dlmo,coef) 1140 endif 1141 1142end function cs_mo_phim 1143 1144function cs_mo_phih(z,dlmo) result(coef) & 1145 bind(C, name='cs_mo_phih') 1146 use, intrinsic :: iso_c_binding 1147 1148 implicit none 1149 1150 real(c_double), value :: z,dlmo 1151 real(c_double) :: coef 1152 double precision :: dlmoneutral = 1.d-12 1153 1154 if (abs(dlmo).lt.dlmoneutral) then 1155 call mo_phih_n(z,dlmo,coef) 1156 elseif (dlmo.ge.0.d0) then 1157 call mo_phih_s(z,dlmo,coef) 1158 else 1159 call mo_phih_u(z,dlmo,coef) 1160 endif 1161 1162end function cs_mo_phih 1163 1164! Integrated version from z0 to z 1165 1166function cs_mo_psim(z,z0,dlmo) result(coef) & 1167 bind(C, name='cs_mo_psim') 1168 use, intrinsic :: iso_c_binding 1169 implicit none 1170 real(c_double), value :: z,z0,dlmo 1171 real(c_double) :: coef 1172 1173 ! Local variable 1174 double precision :: dlmoneutral = 1.d-12 1175 1176 if (abs(dlmo).lt.dlmoneutral) then 1177 call mo_psim_n(z,z0,dlmo,coef) 1178 elseif (dlmo.ge.0.d0) then 1179 call mo_psim_s(z,z0,dlmo,coef) 1180 else 1181 call mo_psim_u(z,z0,dlmo,coef) 1182 endif 1183 1184end function cs_mo_psim 1185 1186function cs_mo_psih (z,z0,dlmo) result(coef) & 1187 bind(C, name='cs_mo_psih') 1188 use, intrinsic :: iso_c_binding 1189 1190 implicit none 1191 1192 real(c_double), value :: z,z0,dlmo 1193 real(c_double) :: coef 1194 double precision :: dlmoneutral = 1.d-12 1195 1196 if (abs(dlmo).lt.dlmoneutral) then 1197 call mo_psih_n(z,z0,dlmo,coef) 1198 elseif (dlmo.ge.0.d0) then 1199 call mo_psih_s(z,z0,dlmo,coef) 1200 else 1201 call mo_psih_u(z,z0,dlmo,coef) 1202 endif 1203 1204end function cs_mo_psih 1205 1206!--------------------------------------------------------------------------- 1207 1208!> \brief Compute LMO, friction velocity ustar, friction temperature 1209!> tstar from a thermal flux using Monin Obukhov 1210 1211!> \param[in] z altitude 1212!> \param[in] z0 1213!> \param[in] du velocity difference 1214!> \param[in] flux thermal flux 1215!> \param[in] tm 1216!> \param[in] gredu 1217!> \param[out] dlmo Inverse Monin Obukhov length 1218!> \param[out] ustar friction velocity 1219 1220subroutine mo_compute_from_thermal_flux(z,z0,du,flux,tm,gredu,dlmo,ustar) 1221 1222 use cstphy 1223 use cstnum 1224 use entsor 1225 1226 implicit none 1227 1228 ! Arguments 1229 double precision z,z0,du,tm,gredu,dlmo,ustar,flux 1230 1231 ! Local variables 1232 double precision tstar 1233 double precision coef_mom 1234 double precision coef_mom_old 1235 double precision prec_ustar 1236 double precision num, denom 1237 double precision :: dlmoclip = 1.0d0 1238 integer icompt 1239 1240 ! Precision initialisation 1241 prec_ustar=1.d-2 1242 1243 icompt=0 1244 1245 ! Initial inverse LMO (neutral) 1246 dlmo = 0.d0 1247 1248 ! Call universal functions 1249 coef_mom = cs_mo_psim((z+z0),z0,dlmo) 1250 1251 ! Initial ustar and tstar 1252 ustar = xkappa * du / coef_mom 1253 tstar = flux / ustar 1254 1255123 continue 1256 1257 icompt=icompt+1 1258 1259 ! Storage previous values 1260 coef_mom_old = coef_mom 1261 1262 ! Update LMO 1263 num = coef_mom**3.d0 * gredu * flux 1264 denom = du**3.d0 * xkappa**2.d0 * tm 1265 if (abs(denom).gt.(epzero*num)) then 1266 dlmo = num / denom 1267 else 1268 dlmo = 0.d0 !FIXME other clipping ? 1269 endif 1270 1271 ! Clipping dlmo (we want |LMO| > 1 m ie 1/|LMO| < 1 m^-1) 1272 if (abs(dlmo).ge.dlmoclip) then 1273 if (dlmo.ge.0.d0) dlmo = dlmoclip 1274 if (dlmo.le.0.d0) dlmo = - dlmoclip 1275 endif 1276 1277 ! Evaluate universal functions 1278 coef_mom = cs_mo_psim((z+z0),z0,dlmo) 1279 1280 ! Update ustar,tstar 1281 ustar = xkappa*du/coef_mom 1282 tstar = flux/ustar 1283 1284 ! Convergence test 1285 if (icompt.le.1000) then 1286 if (abs((coef_mom-coef_mom_old)).ge.prec_ustar) go to 123 1287 endif 1288 1289 return 1290 1291end subroutine mo_compute_from_thermal_flux 1292 1293!--------------------------------------------------------------------------- 1294 1295!> \brief Compute LMO, friction velocity ustar, friction temperature 1296!> tstar from a thermal difference using Monin Obukhov 1297 1298!> \param[in] z altitude 1299!> \param[in] z0 1300!> \param[in] du velocity difference 1301!> \param[in] dt thermal difference 1302!> \param[in] tm 1303!> \param[in] gredu 1304!> \param[out] dlmo Inverse Monin Obukhov length 1305!> \param[out] ustar friction velocity 1306 1307subroutine mo_compute_from_thermal_diff(z,z0,du,dt,tm,gredu,dlmo,ustar) 1308 1309 use cstphy 1310 use cstnum 1311 use entsor 1312 1313 implicit none 1314 1315 ! Arguments 1316 double precision z,z0,du,dt,tm,gredu,dlmo,ustar 1317 1318 ! Local variables 1319 double precision tstar 1320 double precision coef_mom,coef_moh 1321 double precision coef_mom_old,coef_moh_old 1322 double precision prec_ustar,prec_tstar 1323 double precision num, denom 1324 real(c_double) :: zref 1325 double precision :: dlmoclip = 1.0d0 1326 integer icompt 1327 1328 ! Precision initialisation 1329 prec_ustar=1.d-2 1330 prec_tstar=1.d-2 1331 1332 icompt=0 1333 1334 ! Initial LMO 1335 dlmo = 0.d0 1336 1337 ! Call universal functions 1338 zref = z+z0 1339 coef_mom = cs_mo_psim(zref,z0,dlmo) 1340 coef_moh = cs_mo_psih(zref,z0,dlmo) 1341 1342 ! Initial ustar and tstar 1343 ustar = xkappa * du / coef_mom 1344 if (abs(coef_moh).gt.epzero) then 1345 tstar = xkappa*dt/coef_moh 1346 else 1347 tstar = 0.d0 1348 endif 1349 1350 1351123 continue 1352 1353 icompt=icompt+1 1354 1355 ! Storage previous values 1356 coef_mom_old = coef_mom 1357 coef_moh_old = coef_moh 1358 1359 ! Update LMO 1360 num = coef_mom**2.d0 * gredu * dt 1361 denom = (du**2.d0)*tm*coef_moh 1362 if (abs(denom).gt.(epzero* abs(num))) then 1363 dlmo = num / denom 1364 else 1365 dlmo = 0.d0 !FIXME 1366 endif 1367 1368 ! Clipping dlmo (we want |LMO| > 1 m ie 1/|LMO| < 1 m^-1) 1369 if (abs(dlmo).ge.dlmoclip) then 1370 if (dlmo.ge.0.d0) dlmo = dlmoclip 1371 if (dlmo.le.0.d0) dlmo = - dlmoclip 1372 endif 1373 1374 ! Evaluate universal functions 1375 coef_mom = cs_mo_psim((z+z0),z0,dlmo) 1376 coef_moh = cs_mo_psih(z+z0,z0,dlmo) 1377 1378 ! Update ustar,tstar 1379 ustar = xkappa*du/coef_mom 1380 if (abs(coef_moh).gt.epzero) then 1381 tstar=xkappa*dt/coef_moh 1382 else 1383 tstar=0.d0 1384 endif 1385 1386 ! Convergence test 1387 if (icompt.le.1000) then !FIXME compteur max 1000 a mettre en param 1388 if (abs((coef_mom-coef_mom_old)).ge.prec_ustar) go to 123 1389 if (abs((coef_moh-coef_moh_old)).ge.prec_tstar) go to 123 1390 endif 1391 1392 return 1393 1394end subroutine mo_compute_from_thermal_diff 1395 1396end module atincl 1397