Changeset 4478
- Timestamp:
- Mar 27, 2023, 4:25:12 PM (15 months ago)
- Location:
- LMDZ6/trunk/libf
- Files:
-
- 8 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90
r4463 r4478 7 7 subroutine atke_compute_km_kh(ngrid,nlay, & 8 8 wind_u,wind_v,temp,play,pinterf, & 9 tke,Km ,Kh)9 tke,Km_out,Kh_out) 10 10 11 11 !======================================================================== … … 16 16 ! collective and collaborative workshop, 17 17 ! the so-called 'Atelier TKE (ATKE)' with 18 ! K. Arjdal, L. Raillard, C. Dehondt, P. Tiengou, A. Spiga, F. Cheruy, 18 ! K. Arjdal, L. Raillard, C. Dehondt, P. Tiengou, A. Spiga, F. Cheruy, T Dubos, 19 19 ! M. Coulon-Decorzens, S. Fromang, G. Riviere, A. Sima, F. Hourdin, E. Vignon 20 20 ! … … 26 26 27 27 28 USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cn, cinf, rpi 29 USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, rg, rd 28 USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cn, cinf, rpi, rcpd 29 USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, rg, rd 30 USE atke_turbulence_ini_mod, ONLY : viscom, viscoh 30 31 31 32 implicit none … … 47 48 REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT) :: tke ! turbulent kinetic energy at interface between layers 48 49 49 REAL, DIMENSION(ngrid,nlay +1), INTENT(OUT) :: Km !Exchange coefficient for momentum at interface between layers50 REAL, DIMENSION(ngrid,nlay +1), INTENT(OUT) :: Kh !Exchange coefficient for heat flux at interface between layers50 REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: Km_out ! output: Exchange coefficient for momentum at interface between layers 51 REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: Kh_out ! output: Exchange coefficient for heat flux at interface between layers 51 52 52 53 ! Local variables 53 REAL, DIMENSION(ngrid,nlay+1) :: l_exchange ! Length of exchange 54 REAL, DIMENSION(ngrid,nlay+1) :: z_interf ! Altitude at the interface 55 REAL, DIMENSION(ngrid,nlay+1) :: dz_interf ! distance between two consecutive interfaces 56 REAL, DIMENSION(ngrid,nlay+1) :: Ri ! Richardson's number 57 REAL, DIMENSION(ngrid,nlay+1) :: Prandtl ! Turbulent Prandtl's number 58 REAL, DIMENSION(ngrid,nlay+1) :: Sm ! Stability function for momentum 59 REAL, DIMENSION(ngrid,nlay+1) :: Sh ! Stability function for heat 54 REAL, DIMENSION(ngrid,nlay+1) :: Km ! Exchange coefficient for momentum at interface between layers 55 REAL, DIMENSION(ngrid,nlay+1) :: Kh ! Exchange coefficient for heat flux at interface between layers 56 REAL, DIMENSION(ngrid,nlay) :: theta ! Potential temperature 57 REAL, DIMENSION(ngrid,nlay+1) :: l_exchange ! Length of exchange (at interface) 58 REAL, DIMENSION(ngrid,nlay+1) :: z_interf ! Altitude at the interface 59 REAL, DIMENSION(ngrid,nlay) :: z_lay ! Altitude of layers 60 REAL, DIMENSION(ngrid,nlay) :: dz_interf ! distance between two consecutive interfaces 61 REAL, DIMENSION(ngrid,nlay) :: dz_lay ! distance between two layer middles (NB: first and last are half layers) 62 REAL, DIMENSION(ngrid,nlay+1) :: Ri ! Richardson's number (at interface) 63 REAL, DIMENSION(ngrid,nlay+1) :: Prandtl ! Turbulent Prandtl's number (at interface) 64 REAL, DIMENSION(ngrid,nlay+1) :: Sm ! Stability function for momentum (at interface) 65 REAL, DIMENSION(ngrid,nlay+1) :: Sh ! Stability function for heat (at interface) 60 66 61 67 INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid) 62 REAL :: Ri0 ! parameter for Sm stability function 63 68 REAL :: Ri0,Ri1 ! parameter for Sm stability function and Prandlt 69 REAL :: preff ! reference pressure for potential temperature calculations 70 REAL :: thetam ! mean potential temperature at interface 64 71 65 72 … … 68 75 69 76 DO igrid=1,ngrid 77 dz_interf(igrid,1) = 0.0 70 78 z_interf(igrid,1) = 0.0 71 Ri(igrid,1) = 0.1! TO DO 72 Sm(igrid,1) = 0.0 ! TO DO 79 ! Km(igrid,1) = 0.0 80 ! Kh(igrid,1) = 0.0 81 ! tke(igrid,1) = 0.0 82 ! Km(igrid,nlay+1) = 0.0 83 ! Kh(igrid,nlay+1) = 0.0 84 ! tke(igrid,nlay+1) = 0.0 85 END DO 86 87 ! Calculation of potential temperature: (if vapor -> todo virtual potential temperature) 88 !===================================== 89 90 preff=100000. 91 ! The result should not depend on the choice of preff 92 DO ilay=1,nlay 93 DO igrid = 1, ngrid 94 theta(igrid,ilay)=temp(igrid,ilay)*(preff/play(igrid,ilay))**(rd/rcpd) 95 END DO 73 96 END DO 74 97 75 98 76 99 77 ! Calculation of altitude of layers' bottom interfaces:78 !===================================================== 100 ! Calculation of altitude of layers' middle and bottom interfaces: 101 !================================================================= 79 102 80 103 DO ilay=2,nlay+1 81 104 DO igrid=1,ngrid 82 dz_interf(igrid,ilay ) = rg * temp(igrid,ilay) / rg * log(play(igrid,ilay-1)/play(igrid,ilay))83 z_interf(igrid,ilay) = z_interf(igrid,ilay-1) + dz_interf(igrid,ilay )105 dz_interf(igrid,ilay-1) = rd*temp(igrid,ilay-1)/rg/play(igrid,ilay-1)*(pinterf(igrid,ilay-1)-pinterf(igrid,ilay)) 106 z_interf(igrid,ilay) = z_interf(igrid,ilay-1) + dz_interf(igrid,ilay-1) 84 107 ENDDO 108 ENDDO 109 110 DO ilay=1,nlay 111 DO igrid=1,ngrid 112 z_lay(igrid,ilay)=0.5*(z_interf(igrid, ilay+1) + z_interf(igrid, ilay)) 113 ENDDO 85 114 ENDDO 86 115 87 116 88 117 ! Computing the mixing length: 89 !============================= 118 ! so far, we have neglected the effect of local stratification 119 !============================================================== 120 90 121 91 122 DO ilay=2,nlay+1 … … 98 129 !=================================================================== 99 130 100 DO ilay=2,nlay +1131 DO ilay=2,nlay 101 132 DO igrid=1,ngrid 102 ! Warning: sure that gradient of temp and wind should be calculated with dz_interf and not with dz_lay?103 Ri(igrid,ilay) = rg * (temp(igrid,ilay) - temp(igrid,ilay-1)) / dz_interf(igrid,ilay) / &104 (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_interf(igrid,ilay))**2 +&105 ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_interf(igrid,ilay))**2 )106 133 dz_lay(igrid,ilay)=z_lay(igrid,ilay)-z_lay(igrid,ilay-1) 134 thetam=0.5*(theta(igrid,ilay) + theta(igrid,ilay-1)) 135 Ri(igrid,ilay) = rg * (theta(igrid,ilay) - theta(igrid,ilay-1))/thetam / dz_lay(igrid,ilay) / & 136 MAX(((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + & 137 ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2,1E-10) 107 138 ! calculation of Ri0 such that continuity in slope of Sm at Ri=0 108 Ri0=2.*ric*cinf/rpi/cn 139 Ri0=2./rpi*(cinf - cn)*ric/cn 140 ! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0 141 Ri1 = -2./rpi * (pr_asym - pr_neut) / pr_slope 109 142 110 143 IF (Ri(igrid,ilay) < 0.) THEN 111 Sm(igrid,ilay) = 2./rpi * cinf* atan(-Ri(igrid,ilay)/Ri0) + cn112 Prandtl(igrid,ilay) = 1. + Ri(igrid,ilay) * pr_slope144 Sm(igrid,ilay) = 2./rpi * (cinf-cn) * atan(-Ri(igrid,ilay)/Ri0) + cn 145 Prandtl(igrid,ilay) = -2./rpi * (pr_asym - pr_neut) * atan(Ri(igrid,ilay)/Ri1) + pr_neut 113 146 ELSE 114 147 Sm(igrid,ilay) = max(0.,cn*(1.-Ri(igrid,ilay)/Ric)) 115 Prandtl(igrid,ilay) = 2./rpi * atan(-Ri(igrid,ilay)) - 1. + pr_asym ! exp(Ri(igrid,ilay))148 Prandtl(igrid,ilay) = pr_neut + Ri(igrid,ilay) * pr_slope 116 149 END IF 117 150 … … 126 159 127 160 ! stationary solution neglecting the vertical transport of TKE by turbulence 128 DO ilay=2,nlay +1161 DO ilay=2,nlay 129 162 DO igrid=1,ngrid 130 163 tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * & 131 (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_ interf(igrid,ilay))**2 + &132 ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_ interf(igrid,ilay))**2 ) * &164 (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + & 165 ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 ) * & 133 166 (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay)) 134 167 ENDDO … … 145 178 ! Computing eddy diffusivity coefficients: 146 179 !======================================== 147 DO ilay=2,nlay +1180 DO ilay=2,nlay ! TODO: also calculate for nlay+1 ? 148 181 DO igrid=1,ngrid 149 Km(igrid,ilay) = l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5 150 Kh(igrid,ilay) = l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5 182 ! we add the molecular viscosity to Km,h 183 Km(igrid,ilay) = viscom + l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5 184 Kh(igrid,ilay) = viscoh + l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5 151 185 END DO 152 186 END DO 153 187 154 188 ! for output: 189 !=========== 190 Km_out(1:ngrid,2:nlay)=Km(1:ngrid,2:nlay) 191 Kh_out(1:ngrid,2:nlay)=Kh(1:ngrid,2:nlay) 155 192 156 193 end subroutine atke_compute_km_kh -
LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90
r4449 r4478 9 9 real :: kappa = 0.4 ! Von Karman constant 10 10 !$OMP THREADPRIVATE(kappa) 11 real :: l0, ric, ri0, cn, cinf, cepsilon, pr_slope, pr_asym 12 !$OMP THREADPRIVATE(l0, ric, cn, cinf, cepsilon, pr_slope, pr_asym )11 real :: l0, ric, ri0, cn, cinf, cepsilon, pr_slope, pr_asym, pr_neut 12 !$OMP THREADPRIVATE(l0, ric, cn, cinf, cepsilon, pr_slope, pr_asym, pr_neut) 13 13 integer :: lunout,prt_level 14 14 !$OMP THREADPRIVATE(lunout,prt_level) 15 real :: rg, rd, rpi 16 !$OMP THREADPRIVATE(rg, rd, rpi) 15 real :: rg, rd, rpi, rcpd 16 !$OMP THREADPRIVATE(rg, rd, rpi, rcpd) 17 18 real :: viscom, viscoh 19 !$OMP THREADPRIVATE(viscom,viscoh) 20 17 21 18 22 19 23 CONTAINS 20 24 21 SUBROUTINE atke_ini(prt_level_in, lunout_in, rg_in, rd_in, rpi_in )25 SUBROUTINE atke_ini(prt_level_in, lunout_in, rg_in, rd_in, rpi_in, rcpd_in) 22 26 23 27 USE ioipsl_getin_p_mod, ONLY : getin_p 24 28 25 29 integer, intent(in) :: lunout_in,prt_level_in 26 real, intent(in) :: rg_in, rd_in, rpi_in 30 real, intent(in) :: rg_in, rd_in, rpi_in, rcpd_in 27 31 28 32 … … 32 36 rg=rg_in 33 37 rpi=rpi_in 38 rcpd=rcpd_in 39 40 viscom=1.46E-5 41 viscoh=2.06E-5 34 42 35 43 ! flag that controls options in atke_compute_km_kh … … 65 73 CALL getin_p('atke_pr_asym',pr_asym) 66 74 75 ! value of turbulent prandtl number in neutral conditions (Ri=0) 76 ! to guarantee tke>0 in stationary case, pr_neut has to be >= 1 77 pr_neut=1.0 78 CALL getin_p('atke_pr_neut',pr_neut) 79 67 80 68 81 RETURN -
LMDZ6/trunk/libf/phylmd/cdrag_mod.F90
r3817 r4478 16 16 speed, t1, q1, zgeop1, & 17 17 psol, tsurf, qsurf, z0m, z0h, & 18 ri_in, iri_in, & 18 19 cdm, cdh, zri, pref) 19 20 … … 22 23 USE print_control_mod, ONLY: lunout, prt_level 23 24 USE ioipsl_getin_p_mod, ONLY : getin_p 25 USE atke_turbulence_ini_mod, ONLY : ric, cn, cinf, cepsilon, pr_slope, pr_asym, pr_neut 24 26 25 27 IMPLICIT NONE … … 84 86 REAL, DIMENSION(klon), INTENT(IN) :: qsurf ! Surface humidity (Kg/Kg) 85 87 REAL, DIMENSION(klon), INTENT(IN) :: z0m, z0h ! Rugosity at surface (m) 88 REAL, DIMENSION(klon), INTENT(IN) :: ri_in ! Input Richardson 1st layer for first guess calculations of screen var. 89 INTEGER, INTENT(IN) :: iri_in! iflag to activate cdrag calculation using ri1 86 90 REAL, DIMENSION(klon), INTENT(IN) :: t1 ! Temperature au premier niveau (K) 87 91 REAL, DIMENSION(klon), INTENT(IN) :: q1 ! humidite specifique au premier niveau (kg/kg) … … 123 127 REAL, DIMENSION(klon) :: cdmn, cdhn ! Drag coefficient in neutral conditions 124 128 REAL zzzcd 129 REAL, DIMENSION(klon) :: sm, prandtl ! Stability function from atke turbulence scheme 130 REAL ri0, ri1 ! to have dimensionless term in sm and prandtl 125 131 126 132 LOGICAL, SAVE :: firstcall = .TRUE. … … 162 168 163 169 ALPHA=5.0 164 170 165 171 166 172 ! ================================================================= c … … 224 230 *(1.+RETV*max(q1(i),0.0)) ! negative q1 set to zero 225 231 zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd) 232 IF (iri_in.EQ.1) THEN 233 zri(i) = ri_in(i) 234 ENDIF 226 235 227 236 … … 280 289 FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min) 281 290 282 291 CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023) 292 293 ri0=2./rpi*(cinf - cn)*ric/cn 294 ri1=-2./rpi * (pr_asym - pr_neut)/pr_slope 295 sm(i) = 2./rpi * (cinf - cn) * atan(-zri(i)/ri0) + cn 296 prandtl(i) = -2./rpi * (pr_asym - pr_neut) * atan(zri(i)/ri1) + pr_neut 297 FM(i) = MAX((sm(i)**(3/2) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1/2)),f_ri_cd_min) 298 FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min) 283 299 284 300 CASE default ! Louis 1982 … … 376 392 endif 377 393 378 379 380 394 CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023) 395 sm(i) = MAX(0.,cn*(1.-zri(i)/ric)) 396 prandtl(i) = pr_neut + zri(i) * pr_slope 397 FM(i) = MAX((sm(i)**(3/2) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1/2)),f_ri_cd_min) 398 FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min) 381 399 382 400 CASE default ! Louis 1982 … … 414 432 END SUBROUTINE cdrag 415 433 416 !417 SUBROUTINE cdragn_ri1(knon, nsrf, &418 speed, t1, q1, zgeop1, &419 psol, tsurf, qsurf, z0m, z0h, &420 ri1, iri1, &421 cdm, cdh, zri, pref)422 423 USE dimphy424 USE indice_sol_mod425 USE print_control_mod, ONLY: lunout, prt_level426 USE ioipsl_getin_p_mod, ONLY : getin_p427 428 IMPLICIT NONE429 ! ================================================================= c430 !431 ! Objet : calcul des cdrags pour le moment (pcfm) et432 ! les flux de chaleur sensible et latente (pcfh) d'apr??s433 ! Louis 1982, Louis 1979, King et al 2001434 ! ou Zilitinkevich et al 2002 pour les cas stables, Louis 1979435 ! et 1982 pour les cas instables436 !437 ! Modified history:438 ! writting on the 20/05/2016439 ! modified on the 13/12/2016 to be adapted to LMDZ6440 !441 ! References:442 ! Louis, J. F., 1979: A parametric model of vertical eddy fluxes in the443 ! atmosphere. Boundary-Layer Meteorology. 01/1979; 17(2):187-202.444 ! Louis, J. F., Tiedtke, M. and Geleyn, J. F., 1982: `A short history of the445 ! operational PBL parametrization at ECMWF'. Workshop on boundary layer446 ! parametrization, November 1981, ECMWF, Reading, England.447 ! Page: 19. Equations in Table 1.448 ! Mascart P, Noilhan J, Giordani H 1995.A MODIFIED PARAMETERIZATION OF FLUX-PROFILE RELATIONSHIPS449 ! IN THE SURFACE LAYER USING DIFFERENT ROUGHNESS LENGTH VALUES FOR HEAT AND MOMENTUM450 ! Boundary-Layer Meteorology 72: 331-344451 ! Anton Beljaars. May 1992. The parametrization of the planetary boundary layer.452 ! European Centre for Medium-Range Weather Forecasts.453 ! Equations: 110-113. Page 40.454 ! Miller,M.J., A.C.M.Beljaars, T.N.Palmer. 1992. The sensitivity of the ECMWF455 ! model to the parameterization of evaporation from the tropical oceans. J.456 ! Climate, 5:418-434.457 ! King J.C, Connolley, W.M ad Derbyshire S.H. 2001, Sensitivity of Modelled Antarctic climate458 ! to surface and boundary-layer flux parametrizations459 ! QJRMS, 127, pp 779-794460 !461 ! ================================================================= c462 ! ================================================================= c463 ! On choisit le couple de fonctions de correction avec deux flags:464 ! Un pour les cas instables, un autre pour les cas stables465 !466 ! iflag_corr_insta:467 ! 1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h)468 ! 2: Louis 1982469 ! 3: Laurent Li470 !471 ! iflag_corr_sta:472 ! 1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h)473 ! 2: Louis 1982474 ! 3: Laurent Li475 ! 4: King 2001 (SHARP)476 ! 5: MO 1st order theory (allow collapse of turbulence)477 !478 !479 !*****************************************************************480 ! Parametres d'entree481 !*****************************************************************482 483 INTEGER, INTENT(IN) :: knon, nsrf ! nombre de points de grille sur l'horizontal + type de surface484 REAL, DIMENSION(klon), INTENT(IN) :: speed ! module du vent au 1er niveau du modele485 REAL, DIMENSION(klon), INTENT(IN) :: zgeop1! geopotentiel au 1er niveau du modele486 REAL, DIMENSION(klon), INTENT(IN) :: tsurf ! Surface temperature (K)487 REAL, DIMENSION(klon), INTENT(IN) :: qsurf ! Surface humidity (Kg/Kg)488 REAL, DIMENSION(klon), INTENT(IN) :: z0m, z0h ! Rugosity at surface (m)489 REAL, DIMENSION(klon), INTENT(IN) :: ri1 ! Richardson 1ere couche490 INTEGER, INTENT(IN) :: iri1 !491 REAL, DIMENSION(klon), INTENT(IN) :: t1 ! Temperature au premier niveau (K)492 REAL, DIMENSION(klon), INTENT(IN) :: q1 ! humidite specifique au premier niveau (kg/kg)493 REAL, DIMENSION(klon), INTENT(IN) :: psol ! pression au sol494 495 496 497 ! Parametres de sortie498 !******************************************************************499 REAL, DIMENSION(klon), INTENT(OUT) :: cdm ! Drag coefficient for heat flux500 REAL, DIMENSION(klon), INTENT(OUT) :: cdh ! Drag coefficient for momentum501 REAL, DIMENSION(klon), INTENT(OUT) :: zri ! Richardson number502 REAL, DIMENSION(klon), INTENT(OUT) :: pref ! Pression au niveau zgeop/RG503 ! Variables Locales504 !******************************************************************505 506 507 INCLUDE "YOMCST.h"508 INCLUDE "YOETHF.h"509 INCLUDE "clesphys.h"510 511 512 REAL, PARAMETER :: CKAP=0.40, CKAPT=0.42513 REAL CEPDU2514 REAL ALPHA515 REAL CB,CC,CD,C2,C3516 REAL MU, CM, CH, B, CMstar, CHstar517 REAL PM, PH, BPRIME518 REAL C519 INTEGER ng_q1 ! Number of grids that q1 < 0.0520 INTEGER ng_qsurf ! Number of grids that qsurf < 0.0521 INTEGER i522 REAL zdu2, ztsolv523 REAL ztvd, zscf524 REAL zucf, zcr525 REAL friv, frih526 REAL, DIMENSION(klon) :: FM, FH ! stability functions527 REAL, DIMENSION(klon) :: cdmn, cdhn ! Drag coefficient in neutral conditions528 REAL zzzcd529 530 LOGICAL, SAVE :: firstcall = .TRUE.531 !$OMP THREADPRIVATE(firstcall)532 INTEGER, SAVE :: iflag_corr_sta533 !$OMP THREADPRIVATE(iflag_corr_sta)534 INTEGER, SAVE :: iflag_corr_insta535 !$OMP THREADPRIVATE(iflag_corr_insta)536 537 !===================================================================c538 ! Valeurs numeriques des constantes539 !===================================================================c540 541 542 ! Minimum du carre du vent543 544 CEPDU2 = (0.1)**2545 ! Louis 1982546 547 CB=5.0548 CC=5.0549 CD=5.0550 551 552 ! King 2001553 554 C2=0.25555 C3=0.0625556 557 558 ! Louis 1979559 560 BPRIME=4.7561 B=9.4562 563 564 !MO565 566 ALPHA=5.0567 568 569 ! ================================================================= c570 ! Tests avant de commencer571 ! Fuxing WANG, 04/03/2015572 ! To check if there are negative q1, qsurf values.573 !====================================================================c574 ng_q1 = 0 ! Initialization575 ng_qsurf = 0 ! Initialization576 DO i = 1, knon577 IF (q1(i).LT.0.0) ng_q1 = ng_q1 + 1578 IF (qsurf(i).LT.0.0) ng_qsurf = ng_qsurf + 1579 ENDDO580 IF (ng_q1.GT.0 .and. prt_level > 5) THEN581 WRITE(lunout,*)" *** Warning: Negative q1(humidity at 1st level) values in cdrag.F90 !"582 WRITE(lunout,*)" The total number of the grids is: ", ng_q1583 WRITE(lunout,*)" The negative q1 is set to zero "584 ! abort_message="voir ci-dessus"585 ! CALL abort_physic(modname,abort_message,1)586 ENDIF587 IF (ng_qsurf.GT.0 .and. prt_level > 5) THEN588 WRITE(lunout,*)" *** Warning: Negative qsurf(humidity at surface) values in cdrag.F90 !"589 WRITE(lunout,*)" The total number of the grids is: ", ng_qsurf590 WRITE(lunout,*)" The negative qsurf is set to zero "591 ! abort_message="voir ci-dessus"592 ! CALL abort_physic(modname,abort_message,1)593 ENDIF594 595 596 597 !=============================================================================c598 ! Calcul du cdrag599 !=============================================================================c600 601 ! On choisit les fonctions de stabilite utilisees au premier appel602 !**************************************************************************603 IF (firstcall) THEN604 iflag_corr_sta=2605 iflag_corr_insta=2606 607 CALL getin_p('iflag_corr_sta',iflag_corr_sta)608 CALL getin_p('iflag_corr_insta',iflag_corr_insta)609 610 firstcall = .FALSE.611 ENDIF612 613 !xxxxxxxxxxxxxxxxxxxxxxx614 DO i = 1, knon ! Boucle sur l'horizontal615 !xxxxxxxxxxxxxxxxxxxxxxx616 617 618 ! calculs preliminaires:619 !***********************620 621 622 zdu2 = MAX(CEPDU2, speed(i)**2)623 pref(i) = EXP(LOG(psol(i)) - zgeop1(i)/(RD*t1(i)* &624 (1.+ RETV * max(q1(i),0.0)))) ! negative q1 set to zero625 ztsolv = tsurf(i) * (1.0+RETV*max(qsurf(i),0.0)) ! negative qsurf set to zero626 ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*max(q1(i),0.0))) &627 *(1.+RETV*max(q1(i),0.0)) ! negative q1 set to zero628 zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)629 630 IF (iri1.EQ.1) THEN631 zri(i) = ri1(i)632 ENDIF633 634 ! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)):635 !********************************************************************636 637 zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*z0m(i)))638 cdmn(i) = zzzcd*zzzcd639 cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i))))640 641 642 ! Calcul des fonctions de stabilit?? FMs, FHs, FMi, FHi :643 !*******************************************************644 645 !''''''''''''''646 ! Cas instables647 !''''''''''''''648 649 IF (zri(i) .LT. 0.) THEN650 651 652 SELECT CASE (iflag_corr_insta)653 654 CASE (1) ! Louis 1979 + Mascart 1995655 656 MU=LOG(MAX(z0m(i)/z0h(i),0.01))657 CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3)658 PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3)659 CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3)660 PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3)661 CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &662 & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i))) &663 & * ((zgeop1(i)/(RG*z0h(i)))**PH)664 CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &665 & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &666 & * ((zgeop1(i)/(RG*z0m(i)))**PM)667 668 669 670 671 FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i))))672 FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i))))673 674 CASE (2) ! Louis 1982675 676 zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &677 *(1.0+zgeop1(i)/(RG*z0m(i)))))678 FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)679 FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)680 681 682 CASE (3) ! Laurent Li683 684 685 FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)686 FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)687 688 689 690 CASE default ! Louis 1982691 692 zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &693 *(1.0+zgeop1(i)/(RG*z0m(i)))))694 FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)695 FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)696 697 698 END SELECT699 700 701 702 ! Calcul des drags703 704 705 cdm(i)=cdmn(i)*FM(i)706 cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)707 708 709 ! Traitement particulier des cas oceaniques710 ! on applique Miller et al 1992 en l'absence de gustiness711 712 IF (nsrf == is_oce) THEN713 ! cdh(i)=f_cdrag_oce*cdhn(i)*FH(i)714 715 IF(iflag_gusts==0) THEN716 zcr = (0.0016/(cdmn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)717 cdh(i) =f_cdrag_oce* cdhn(i)*(1.0+zcr**1.25)**(1./1.25)718 ENDIF719 720 721 cdm(i)=MIN(cdm(i),cdmmax)722 cdh(i)=MIN(cdh(i),cdhmax)723 724 END IF725 726 727 728 ELSE729 730 !'''''''''''''''731 ! Cas stables :732 !'''''''''''''''733 zri(i) = MIN(20.,zri(i))734 735 SELECT CASE (iflag_corr_sta)736 737 CASE (1) ! Louis 1979 + Mascart 1995738 739 FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min)740 FH(i)=FM(i)741 742 743 CASE (2) ! Louis 1982744 745 zscf = SQRT(1.+CD*ABS(zri(i)))746 FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)747 FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )748 749 750 CASE (3) ! Laurent Li751 752 FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min)753 FH(i)=FM(i)754 755 756 CASE (4) ! King 2001757 758 if (zri(i) .LT. C2/2.) then759 FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min)760 FH(i)= FM(i)761 762 763 else764 FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min)765 FH(i)= FM(i)766 endif767 768 769 CASE (5) ! MO770 771 if (zri(i) .LT. 1./alpha) then772 773 FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min)774 FH(i)=FM(i)775 776 else777 778 779 FM(i)=MAX(1E-7,f_ri_cd_min)780 FH(i)=FM(i)781 782 endif783 784 785 786 787 788 CASE default ! Louis 1982789 790 zscf = SQRT(1.+CD*ABS(zri(i)))791 FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)792 FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )793 794 795 796 END SELECT797 798 ! Calcul des drags799 800 801 cdm(i)=cdmn(i)*FM(i)802 cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)803 804 IF(nsrf.EQ.is_oce) THEN805 806 cdh(i)=f_cdrag_oce*cdhn(i)*FH(i)807 cdm(i)=MIN(cdm(i),cdmmax)808 cdh(i)=MIN(cdh(i),cdhmax)809 810 ENDIF811 812 813 ENDIF814 815 !xxxxxxxxxxx816 END DO ! Fin de la boucle sur l'horizontal817 !xxxxxxxxxxx818 ! ================================================================= c819 820 END SUBROUTINE cdragn_ri1821 822 434 END MODULE cdrag_mod -
LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90
r4449 r4478 876 876 REAL, DIMENSION(klon) :: yrmu0 877 877 ! Martin 878 REAL, DIMENSIOn(klon) :: yri0 878 879 879 880 REAL, DIMENSION(klon):: ydelta_sst, ydelta_sal, yds_ns, ydt_ns, ydter, & … … 1060 1061 1061 1062 ytke=0. 1063 yri0(:)=0. 1062 1064 !FC 1063 1065 y_treedrg=0. … … 1569 1571 CALL cdrag(knon, nsrf, & 1570 1572 speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1),& 1571 yts, yqsurf, yz0m, yz0h, &1573 yts, yqsurf, yz0m, yz0h, yri0, 0, & 1572 1574 ycdragm, ycdragh, zri1, pref ) 1573 1575 … … 1603 1605 CALL cdrag(knon, nsrf, & 1604 1606 speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),& 1605 yts_x, yqsurf_x, yz0m, yz0h, &1607 yts_x, yqsurf_x, yz0m, yz0h, yri0, 0, & 1606 1608 ycdragm_x, ycdragh_x, zri1_x, pref_x ) 1607 1609 … … 1630 1632 CALL cdrag(knon, nsrf, & 1631 1633 speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),& 1632 yts_w, yqsurf_w, yz0m, yz0h, &1634 yts_w, yqsurf_w, yz0m, yz0h, yri0, 0, & 1633 1635 ycdragm_w, ycdragh_w, zri1_w, pref_w ) 1634 1636 ! -
LMDZ6/trunk/libf/phylmd/physiq_mod.F90
r4466 r4478 1748 1748 CALL wake_ini(rg,rd,rv,prt_level) 1749 1749 CALL yamada_ini(klon,lunout,prt_level) 1750 CALL atke_ini(prt_level, lunout, RG, RD, RPI )1750 CALL atke_ini(prt_level, lunout, RG, RD, RPI, RCPD) 1751 1751 CALL thermcell_ini(iflag_thermals,prt_level,tau_thermals,lunout, & 1752 1752 & RG,RD,RCPD,RKAPPA,RLVTT,RETV) -
LMDZ6/trunk/libf/phylmd/screenc_mod.F90
r3821 r4478 70 70 ! Variables locales 71 71 INTEGER :: i 72 REAL, dimension(klon) :: cdram, cdrah, cdran, zri1, gref,ycdragm 72 REAL, dimension(klon) :: cdram, cdrah, cdran, zri1, gref,ycdragm,zri_zero 73 73 ! 74 74 !------------------------------------------------------------------------- … … 89 89 speed, temp, q_zref, gref, & 90 90 psol, ts, qsurf, z0m, z0h, & 91 zri_zero,0, & 91 92 cdram, cdrah, zri1, pref) 92 93 DO i = 1, knon … … 177 178 ! Richardson at reference level 178 179 ! 179 CALL cdrag n_ri1(knon, nsrf, &180 CALL cdrag(knon, nsrf, & 180 181 speed, temp, q_zref, gref, & 181 182 psol, ts, qsurf, z0m, z0h, & -
LMDZ6/trunk/libf/phylmd/stdlevvar_mod.F90
r3839 r4478 103 103 REAL, dimension(klon) :: te_zref_con, q_zref_con 104 104 REAL, dimension(klon) :: u_zref_c, te_zref_c, temp_c, q_zref_c 105 REAL, dimension(klon) :: ok_pred, ok_corr 105 REAL, dimension(klon) :: ok_pred, ok_corr, zri_zero 106 106 ! REAL, dimension(klon) :: conv_te, conv_q 107 107 !------------------------------------------------------------------------- … … 121 121 & speed, t1, q1, z1, & 122 122 & psol, ts1, qsurf, z0m, z0h, & 123 & zri_zero, 0, & 123 124 & cdram, cdrah, zri1, pref) 124 125 … … 413 414 REAL, dimension(klon) :: cdrm10m, cdrh10m, ri10m 414 415 REAL, dimension(klon) :: cdmn10m, cdhn10m, fm10m, fh10m 415 REAL, dimension(klon) :: cdn2m, cdn1 416 REAL, dimension(klon) :: cdn2m, cdn1, zri_zero 416 417 REAL :: CEPDUE,zdu2 417 418 INTEGER :: nzref, nz1 … … 444 445 & speed, t1, q1, z1, & 445 446 & psol, ts1, qsurf, z0m, z0h, & 447 & zri_zero, 0, & 446 448 & cdram, cdrah, zri1, pref) 447 449 -
LMDZ6/trunk/libf/phylmdiso/pbl_surface_mod.F90
r4449 r4478 1030 1030 REAL, DIMENSION(klon) :: yrmu0 1031 1031 ! Martin 1032 1032 REAL, DIMENSION(klon) :: yri0 1033 1033 REAL, DIMENSION(klon):: ydelta_sst, ydelta_sal, yds_ns, ydt_ns, ydter, & 1034 1034 ydser, ydt_ds, ytkt, ytks, ytaur, ysss … … 1279 1279 1280 1280 ytke=0. 1281 yri0(:)=0. 1281 1282 !FC 1282 1283 y_treedrg=0. … … 1849 1850 CALL cdrag(knon, nsrf, & 1850 1851 speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1),& 1851 yts, yqsurf, yz0m, yz0h, &1852 yts, yqsurf, yz0m, yz0h, yri0, 0, & 1852 1853 ycdragm, ycdragh, zri1, pref ) 1853 1854 … … 1883 1884 CALL cdrag(knon, nsrf, & 1884 1885 speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),& 1885 yts_x, yqsurf_x, yz0m, yz0h, &1886 yts_x, yqsurf_x, yz0m, yz0h, yri0, 0, & 1886 1887 ycdragm_x, ycdragh_x, zri1_x, pref_x ) 1887 1888 … … 1910 1911 CALL cdrag(knon, nsrf, & 1911 1912 speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),& 1912 yts_w, yqsurf_w, yz0m, yz0h, &1913 yts_w, yqsurf_w, yz0m, yz0h, yri0, 0, & 1913 1914 ycdragm_w, ycdragh_w, zri1_w, pref_w ) 1914 1915 !
Note: See TracChangeset
for help on using the changeset viewer.