Changeset 4738
- Timestamp:
- Oct 24, 2023, 12:58:23 PM (13 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/add_phys_tend_mod.F90
r4523 r4738 71 71 zzdt(i, k) = hthturb_gcssold(k)*dtime_frcg 72 72 zzdq(i, k) = hqturb_gcssold(k)*dtime_frcg 73 END 74 END 73 ENDDO 74 ENDDO 75 75 PRINT *, ' add_pbl_tend, dtime_frcg ', dtime_frcg 76 76 PRINT *, ' add_pbl_tend, zzdt ', zzdt … … 79 79 ELSE 80 80 CALL add_phys_tend(zdu, zdv, zdt, zdq, zdql, zdqi, zdqbs, paprs, text,abortphy,flag_inhib_tend, itap, 0) 81 END IF 82 81 ENDIF 83 82 84 83 RETURN … … 110 109 & , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_qbs_col, d_h_col 111 110 IMPLICIT none 112 include"YOMCST.h"113 include"clesphys.h"111 INCLUDE "YOMCST.h" 112 INCLUDE "clesphys.h" 114 113 115 114 ! Arguments : … … 148 147 LOGICAL done(klon) 149 148 150 integerdebug_level151 logical, save:: first=.true.149 INTEGER debug_level 150 LOGICAL, SAVE :: first=.true. 152 151 !$OMP THREADPRIVATE(first) 153 152 ! … … 179 178 REAL zh_qw_col(klon,2), zh_ql_col(klon,2), zh_qs_col(klon,2), zh_qbs_col(klon,2) 180 179 REAL zh_col(klon,2) 181 182 180 REAL zcpvap, zcwat, zcice 183 181 … … 185 183 ! Initialisations 186 184 187 IF (prt_level >= 5) then 188 write (*,*) "In add_phys_tend, after ",text 189 call flush 190 end if 191 192 ! if flag_inhib_tend != 0, tendencies are not added 193 IF (flag_inhib_tend /= 0) then 194 ! If requiered, diagnostics are shown 195 IF (flag_inhib_tend > 0) then 196 ! print some diagnostics if xxx_seri have changed 197 call cmp_seri(flag_inhib_tend,text) 198 END IF 199 RETURN ! on n ajoute pas les tendance 200 END IF 201 202 IF (abortphy==1) RETURN ! on n ajoute pas les tendance si le modele 203 ! a deja plante. 204 205 debug_level=10 206 if (first) then 207 print *,"TestJLD rcpv, rcw, rcs",rcpv, rcw, rcs 208 first=.false. 209 endif 185 IF (prt_level >= 5) THEN 186 write (*,*) "In add_phys_tend, after ",text 187 CALL flush 188 ENDIf 189 190 ! if flag_inhib_tend != 0, tendencies are not added 191 IF (flag_inhib_tend /= 0) THEN 192 ! If requiered, diagnostics are shown 193 IF (flag_inhib_tend > 0) THEN 194 ! print some diagnostics if xxx_seri have changed 195 call cmp_seri(flag_inhib_tend,text) 196 ENDIF 197 RETURN ! on n ajoute pas les tendance 198 ENDIF 199 200 IF (abortphy==1) RETURN ! on n ajoute pas les tendance si le modele a deja plante. 201 202 debug_level=10 203 IF (first) THEN 204 print *,"TestJLD rcpv, rcw, rcs",rcpv, rcw, rcs 205 first=.false. 206 ENDIF 210 207 ! 211 208 ! print *,'add_phys_tend: paprs ',paprs … … 227 224 ! layer air mass 228 225 zairm(:, k) = (paprs(:,k)-paprs(:,k+1))/rg 229 END 230 231 if (fl_ebil .GT. 0) then226 ENDDO 227 228 IF (fl_ebil .GT. 0) THEN 232 229 ! ------------------------------------------------ 233 230 ! Compute vertical sum for each atmospheric column … … 244 241 zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_qbs_col(:,n), zh_col(:,n)) 245 242 246 end if ! endif (fl_ebil .GT. 0)243 ENDIF ! endif (fl_ebil .GT. 0) 247 244 248 245 !====================================================================== … … 250 247 !====================================================================== 251 248 252 253 254 255 256 249 u_seri(:,:)=u_seri(:,:)+zdu(:,:) 250 v_seri(:,:)=v_seri(:,:)+zdv(:,:) 251 ql_seri(:,:)=ql_seri(:,:)+zdql(:,:) 252 qs_seri(:,:)=qs_seri(:,:)+zdqi(:,:) 253 qbs_seri(:,:)=qbs_seri(:,:)+zdqbs(:,:) 257 254 258 255 !====================================================================== … … 261 258 !====================================================================== 262 259 263 264 265 266 267 268 269 IF ( zt>370. .or. zt<130. .or. abs(zdt(i,k))>50. ) then270 271 272 273 274 IF ( zq<0. .or. zq>0.1 .or. abs(zdq(i,k))>1.e-2 ) then275 276 277 278 279 280 281 282 260 jbad=0 261 jqbad=0 262 DO k = 1, klev 263 DO i = 1, klon 264 zt=t_seri(i,k)+zdt(i,k) 265 zq=q_seri(i,k)+zdq(i,k) 266 IF ( zt>370. .or. zt<130. .or. abs(zdt(i,k))>50. ) THEN 267 jbad = jbad + 1 268 jadrs(jbad) = i 269 kadrs(jbad) = k 270 ENDIF 271 IF ( zq<0. .or. zq>0.1 .or. abs(zdq(i,k))>1.e-2 ) THEN 272 jqbad = jqbad + 1 273 jqadrs(jqbad) = i 274 kqadrs(jqbad) = k 275 ENDIF 276 t_seri(i,k)=zt 277 q_seri(i,k)=zq 278 ENDDO 279 ENDDO 283 280 284 281 !===================================================================================== … … 286 283 !===================================================================================== 287 284 288 IF (jbad .GT. 0) THEN289 290 291 if(prt_level.ge.debug_level) THEN285 IF (jbad .GT. 0) THEN 286 DO j = 1, jbad 287 i=jadrs(j) 288 IF (prt_level.ge.debug_level) THEN 292 289 print*,'PLANTAGE POUR LE POINT i lon lat =',& 293 290 i,longitude_deg(i),latitude_deg(i),text … … 297 294 ENDDO 298 295 call print_debug_phys(i,debug_level,text) 299 endif300 301 ENDIF296 ENDIF 297 ENDDO 298 ENDIF 302 299 ! 303 300 !===================================================================================== 304 301 ! Impression, warning et correction en cas de probleme moins important 305 302 !===================================================================================== 306 IF (jqbad .GT. 0) THEN303 IF (jqbad .GT. 0) THEN 307 304 done(:) = .false. !jyg 308 305 DO j = 1, jqbad … … 315 312 write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k) 316 313 ENDDO 317 endif314 ENDIF 318 315 IF (ok_conserv_q) THEN 319 316 !jyg<20140228 Corrections pour conservation de l'eau … … 328 325 zqp_int = zqp_int + zqp(k) *(paprs(i,k)-paprs(i,k+1))/Rg 329 326 ENDDO 330 if(prt_level.ge.debug_level) THEN327 IF (prt_level.ge.debug_level) THEN 331 328 print*,' cas q_seri<1.e-15 i k zq_int zqp_int zq_int/zqp_int :', & 332 329 i, kqadrs(j), zq_int, zqp_int, zq_int/zqp_int 333 endif330 ENDIF 334 331 DO k = 1, klev 335 332 zq_new = zqp(k)*zq_int/zqp_int … … 343 340 DO k = 1, klev 344 341 zq=q_seri(i,k)+zdq(i,k) 345 if (zq.lt.1.e-15) then346 if (q_seri(i,k).lt.1.e-15) then347 if(prt_level.ge.debug_level) THEN348 print*,' cas q_seri<1.e-15 i k q_seri zq zdq :',i,k,q_seri(i,k),zq,zdq(i,k)349 endif350 q_seri(i,k)=1.e-15351 zdq(i,k)=(1.e-15-q_seri(i,k))352 endif353 endif342 IF (zq.lt.1.e-15) THEN 343 IF (q_seri(i,k).lt.1.e-15) THEN 344 IF (prt_level.ge.debug_level) THEN 345 print*,' cas q_seri<1.e-15 i k q_seri zq zdq :',i,k,q_seri(i,k),zq,zdq(i,k) 346 ENDIF 347 q_seri(i,k)=1.e-15 348 zdq(i,k)=(1.e-15-q_seri(i,k)) 349 ENDIF 350 ENDIF 354 351 ! zq=q_seri(i,k)+zdq(i,k) 355 ! if (zq.lt.1.e-15) then352 ! if (zq.lt.1.e-15) THEN 356 353 ! zdq(i,k)=(1.e-15-q_seri(i,k)) 357 354 ! endif … … 361 358 !jyg> 362 359 ENDDO ! j = 1, jqbad 363 ENDIF360 ENDIF 364 361 ! 365 362 366 363 !IM ajout memes tests pour reverifier les jbad, jqbad beg 367 368 369 370 371 IF ( t_seri(i,k)>370. .or. t_seri(i,k)<130. .or. abs(zdt(i,k))>50. ) then372 373 374 ! 375 ! 376 ! 377 378 IF ( q_seri(i,k)<0. .or. q_seri(i,k)>0.1 .or. abs(zdq(i,k))>1.e-2 ) then379 380 381 382 ! 383 ! 384 ! 385 386 387 388 IF (jbad .GT. 0) THEN364 jbad=0 365 jqbad=0 366 DO k = 1, klev 367 DO i = 1, klon 368 IF ( t_seri(i,k)>370. .or. t_seri(i,k)<130. .or. abs(zdt(i,k))>50. ) THEN 369 jbad = jbad + 1 370 jadrs(jbad) = i 371 ! if(prt_level.ge.debug_level) THEN 372 ! print*,'cas2 i k t_seri zdt',i,k,t_seri(i,k),zdt(i,k) 373 ! endif 374 ENDIF 375 IF ( q_seri(i,k)<0. .or. q_seri(i,k)>0.1 .or. abs(zdq(i,k))>1.e-2 ) THEN 376 jqbad = jqbad + 1 377 jqadrs(jqbad) = i 378 kqadrs(jqbad) = k 379 ! if(prt_level.ge.debug_level) THEN 380 ! print*,'cas2 i k q_seri zdq',i,k,q_seri(i,k),zdq(i,k) 381 ! endif 382 ENDIF 383 ENDDO 384 ENDDO 385 IF (jbad .GT. 0) THEN 389 386 DO j = 1, jbad 390 387 i=jadrs(j) … … 400 397 ENDDO 401 398 call print_debug_phys(i,debug_level,text) 402 endif399 ENDIF 403 400 ENDDO 404 ENDIF405 ! 406 IF (jqbad .GT. 0) THEN401 ENDIF 402 ! 403 IF (jqbad .GT. 0) THEN 407 404 DO j = 1, jqbad 408 405 i=jqadrs(j) 409 406 k=kqadrs(j) 410 if(prt_level.ge.debug_level) THEN407 IF (prt_level.ge.debug_level) THEN 411 408 print*,'WARNING : EAU2 POUR LE POINT i itap lon lat txt jqbad zdq q zdql ql',& 412 409 i,itap,longitude_deg(i),latitude_deg(i),text,jqbad,& … … 418 415 ENDDO 419 416 call print_debug_phys(i,debug_level,text) 420 endif417 ENDIF 421 418 ENDDO 422 ENDIF419 ENDIF 423 420 424 421 !====================================================================== … … 429 426 !====================================================================== 430 427 431 432 433 Print*,'ERROR ABORT hgardfou dans ',text428 CALL hgardfou(t_seri,ftsol,text,abortphy) 429 IF (abortphy==1) THEN 430 print*,'ERROR ABORT hgardfou dans ',text 434 431 ! JLD pourquoi on ne modifie pas de meme t_seri et q_seri ? 435 436 437 438 439 440 432 u_seri(:,:)=u_seri(:,:)-zdu(:,:) 433 v_seri(:,:)=v_seri(:,:)-zdv(:,:) 434 ql_seri(:,:)=ql_seri(:,:)-zdql(:,:) 435 qs_seri(:,:)=qs_seri(:,:)-zdqi(:,:) 436 qbs_seri(:,:)=qbs_seri(:,:)-zdqbs(:,:) 437 ENDIF 441 438 442 439 !====================================================================== … … 444 441 !====================================================================== 445 442 446 if (fl_ebil .GT. 0) then443 IF (fl_ebil .GT. 0) THEn 447 444 448 445 ! ------------------------------------------------ … … 476 473 d_h_col = (zh_col(:,2)-zh_col(:,1))/phys_tstep 477 474 478 end if ! endif (fl_ebil .GT. 0)475 ENDIF ! endif (fl_ebil .GT. 0) 479 476 ! 480 477 ! When in diagnostic mode, restore "out" variables to initial values. … … 530 527 REAL, DIMENSION(nlon,nlev) :: uu_n, vv_n 531 528 REAL, DIMENSION(nlon,nlev) :: temp_n, qv_n, ql_n, qs_n, qbs_n 532 533 534 529 ! 535 530 INTEGER k, n 536 531 537 integerdebug_level538 logical, save:: first=.true.532 INTEGER debug_level 533 LOGICAL, SAVE :: first=.true. 539 534 !$OMP THREADPRIVATE(first) 540 535 ! … … 569 564 ! Initialisations 570 565 571 IF (prt_level >= 5) then572 573 callflush574 end if575 576 577 if (first) then578 579 580 endif566 IF (prt_level >= 5) THEN 567 write (*,*) "In diag_phys_tend, after ",text 568 CALL flush 569 ENDIF 570 571 debug_level=10 572 IF (first) THEN 573 print *,"TestJLD rcpv, rcw, rcs",rcpv, rcw, rcs 574 first=.false. 575 ENDIF 581 576 ! 582 577 ! print *,'add_phys_tend: paprs ',paprs … … 587 582 ! layer air mass 588 583 zairm(:, k) = (paprs(:,k)-paprs(:,k+1))/rg 589 END 590 591 if (fl_ebil .GT. 0) then584 ENDDO 585 586 IF (fl_ebil .GT. 0) THEN 592 587 ! ------------------------------------------------ 593 588 ! Compute vertical sum for each atmospheric column … … 600 595 zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_qbs_col(:,n), zh_col(:,n)) 601 596 602 end if ! endif (fl_ebil .GT. 0)597 ENDIF ! endif (fl_ebil .GT. 0) 603 598 604 599 !====================================================================== … … 606 601 !====================================================================== 607 602 608 uu_n(:,:)=uu(:,:)+zdu(:,:) 609 vv_n(:,:)=vv(:,:)+zdv(:,:) 610 qv_n(:,:)=qv(:,:)+zdq(:,:) 611 ql_n(:,:)=ql(:,:)+zdql(:,:) 612 qs_n(:,:)=qs(:,:)+zdqs(:,:) 613 qbs_n(:,:)=qbs(:,:)+zdqbs(:,:) 614 temp_n(:,:)=temp(:,:)+zdt(:,:) 615 616 603 uu_n(:,:)=uu(:,:)+zdu(:,:) 604 vv_n(:,:)=vv(:,:)+zdv(:,:) 605 qv_n(:,:)=qv(:,:)+zdq(:,:) 606 ql_n(:,:)=ql(:,:)+zdql(:,:) 607 qs_n(:,:)=qs(:,:)+zdqs(:,:) 608 qbs_n(:,:)=qbs(:,:)+zdqbs(:,:) 609 temp_n(:,:)=temp(:,:)+zdt(:,:) 617 610 618 611 !====================================================================== … … 620 613 !====================================================================== 621 614 622 if (fl_ebil .GT. 0) then615 IF (fl_ebil .GT. 0) THEN 623 616 624 617 ! ------------------------------------------------ … … 644 637 d_ek_col(:) = (zek_col(:,2)-zek_col(:,1))/phys_tstep 645 638 646 print *,'zdu ', zdu647 print *,'zdv ', zdv648 print *,'d_ek_col, zek_col(2), zek_col(1) ',d_ek_col(1), zek_col(1,2), zek_col(1,1)639 print *,'zdu ', zdu 640 print *,'zdv ', zdv 641 print *,'d_ek_col, zek_col(2), zek_col(1) ',d_ek_col(1), zek_col(1,2), zek_col(1,1) 649 642 650 643 d_h_dair_col(:) = (zh_dair_col(:,2)-zh_dair_col(:,1))/phys_tstep … … 656 649 d_h_col = (zh_col(:,2)-zh_col(:,1))/phys_tstep 657 650 658 end if ! endif (fl_ebil .GT. 0)651 ENDIF ! endif (fl_ebil .GT. 0) 659 652 ! 660 653 … … 668 661 669 662 IMPLICIT none 670 include"YOMCST.h"663 INCLUDE "YOMCST.h" 671 664 672 665 INTEGER, INTENT(IN) :: nlon,nlev … … 685 678 686 679 INTEGER :: i, k 687 688 680 689 681 ! Reset variables … … 720 712 zh_qs_col(i) = zh_qs_col(i) + (zcpvap*temp(i, k) - rlstt)*qs(i, k)*zairm(i, k) !jyg 721 713 zh_qbs_col(i) = zh_qbs_col(i) + (zcpvap*temp(i, k) - rlstt)*qbs(i, k)*zairm(i, k) !jyg 722 END 723 END 714 ENDDO 715 ENDDO 724 716 ! compute total air enthalpy 725 717 zh_col(:) = zh_dair_col(:) + zh_qw_col(:) + zh_ql_col(:) + zh_qs_col(:) + zh_qbs_col(:) … … 749 741 USE climb_hq_mod, ONLY : d_h_col_vdf, f_h_bnd 750 742 IMPLICIT none 751 include"YOMCST.h"743 INCLUDE "YOMCST.h" 752 744 753 745 ! Arguments : 754 746 !------------ 755 747 CHARACTER*(*) text ! text specifing the involved parametrization 756 integeritap ! time step number748 INTEGER itap ! time step number 757 749 ! local variables 758 750 ! --------------- 759 realbilq_seuil, bilh_seuil ! thresold on error in Q and H budget760 realbilq_error, bilh_error ! erros in Q and H budget761 realbilq_bnd, bilh_bnd ! Q and H budget due to exchange with boundaries762 integerbilq_ok, bilh_ok751 REAL bilq_seuil, bilh_seuil ! thresold on error in Q and H budget 752 REAL bilq_error, bilh_error ! erros in Q and H budget 753 REAL bilq_bnd, bilh_bnd ! Q and H budget due to exchange with boundaries 754 INTEGER bilq_ok, bilh_ok 763 755 CHARACTER*(12) status 764 756 … … 769 761 770 762 !!print *,'prt_level:',prt_level,' fl_ebil:',fl_ebil,' fl_cor_ebil:',fl_cor_ebil 771 if ( (fl_ebil .GT. 0) .and. (klon .EQ. 1)) then 763 IF ((fl_ebil .GT. 0) .AND. (klon .EQ. 1)) THEN 772 764 773 765 bilq_bnd = 0. … … 801 793 bilh_error = d_h_col(1) - bilh_bnd 802 794 ! are the errors too large? 803 if ( abs(bilq_error) .gt. bilq_seuil) bilq_ok=1804 if ( abs(bilh_error) .gt. bilh_seuil) bilh_ok=1795 IF (abs(bilq_error) .GT. bilq_seuil) bilq_ok=1 796 IF (abs(bilh_error) .GT. bilh_seuil) bilh_ok=1 805 797 ! 806 798 ! Print diagnostics 807 799 ! ================= 808 if ( (bilq_ok .eq. 0).and.(bilh_ok .eq. 0) ) then800 IF ( (bilq_ok .eq. 0).AND.(bilh_ok .eq. 0) ) THEN 809 801 status="enerbil-OK" 810 else802 ELSE 811 803 status="enerbil-PB" 812 end if813 814 if ( prt_level .GE. 3) then804 ENDIF 805 806 IF (prt_level .GE. 3) THEN 815 807 write(*,9010) text,status," itap:",itap,"enerbilERROR: Q", bilq_error," H", bilh_error 816 808 9010 format (1x,A8,2x,A12,A6,I4,A18,E15.6,A5,E15.6) 817 end if818 if ( prt_level .GE. 3) then809 ENDIF 810 IF (prt_level .GE. 3) THEN 819 811 write(*,9000) text,"enerbil: Q,H,KE budget", d_qt_col(1), d_h_col(1),d_ek_col(1) 820 end if821 if ( prt_level .GE. 5) then812 ENDIF 813 IF (prt_level .GE. 5) THEN 822 814 write(*,9000) text,"enerbil at boundaries: Q, H",bilq_bnd, bilh_bnd 823 815 write(*,9000) text,"enerbil: water budget",d_qt_col(1),d_qw_col(1),d_ql_col(1),d_qs_col(1), d_qbs_col(1) 824 816 write(*,9000) text,"enerbil: enthalpy budget",d_h_col(1),d_h_dair_col(1),d_h_qw_col(1),d_h_ql_col(1),d_h_qs_col(1),d_h_qbs_col(1) 825 end if817 ENDIF 826 818 827 819 specific_diag: SELECT CASE (text) 828 820 CASE("vdf") specific_diag 829 if ( prt_level .GE. 5) then821 IF (prt_level .GE. 5) THEN 830 822 write(*,9000) text,"enerbil: d_h, bilh, sens,t_seri", d_h_col(1), bilh_bnd, sens(1), t_seri(1,1) 831 823 write(*,9000) text,"enerbil: d_h_col_vdf, f_h, diff",d_h_col_vdf, f_h_bnd, bilh_bnd-sens(1) 832 end if824 ENDIF 833 825 CASE("lsc") specific_diag 834 if ( prt_level .GE. 5) then826 IF (prt_level .GE. 5) THEN 835 827 write(*,9000) text,"enerbil: rain, bil_lat, bil_sens", rain_lsc(1), rlvtt * rain_lsc(1), -(rcw-rcpd)*t_seri(1,1) * rain_lsc(1) 836 828 write(*,9000) text,"enerbil: snow, bil_lat, bil_sens", snow_lsc(1), rlstt * snow_lsc(1), -(rcs-rcpd)*t_seri(1,1) * snow_lsc(1) 837 end if829 ENDIF 838 830 CASE("convection") specific_diag 839 if ( prt_level .GE. 5) then831 IF (prt_level .GE. 5) THEN 840 832 write(*,9000) text,"enerbil: rain, bil_lat, bil_sens", rain_con(1), rlvtt * rain_con(1), -(rcw-rcpd)*t_seri(1,1) * rain_con(1) 841 833 write(*,9000) text,"enerbil: snow, bil_lat, bil_sens", snow_con(1), rlstt * snow_con(1), -(rcs-rcpd)*t_seri(1,1) * snow_con(1) 842 end if834 ENDIF 843 835 END SELECT specific_diag 844 836 845 9000 format(1x,A8,2x,A35,10E15.6)846 847 end if ! endif (fl_ebil .GT. 0)837 9000 FORMAT(1x,A8,2x,A35,10E15.6) 838 839 ENDIF ! endif (fl_ebil .GT. 0) 848 840 849 841 END SUBROUTINE prt_enerbil
Note: See TracChangeset
for help on using the changeset viewer.