Changeset 4478 for LMDZ6/trunk/libf/phylmd/cdrag_mod.F90
- Timestamp:
- Mar 27, 2023, 4:25:12 PM (18 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
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
Note: See TracChangeset
for help on using the changeset viewer.