Changeset 2720 for LMDZ5/branches/testing/libf/phylmd/yamada_c.F90
- Timestamp:
- Nov 30, 2016, 1:28:41 PM (8 years ago)
- Location:
- LMDZ5/branches/testing
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/branches/testing
- Property svn:mergeinfo changed
/LMDZ5/trunk merged: 2665-2668,2670-2674,2677-2681,2683-2684,2686,2690-2719
- Property svn:mergeinfo changed
-
LMDZ5/branches/testing/libf/phylmd/yamada_c.F90
r2408 r2720 4 4 SUBROUTINE yamada_c(ngrid,timestep,plev,play & 5 5 & ,pu,pv,pt,d_u,d_v,d_t,cd,q2,km,kn,kq,d_t_diss,ustar & 6 & ,iflag_pbl ,okiophys)6 & ,iflag_pbl) 7 7 USE dimphy, ONLY: klon, klev 8 8 USE print_control_mod, ONLY: prt_level 9 USE ioipsl_getin_p_mod, ONLY : getin_p 10 9 11 IMPLICIT NONE 10 12 #include "YOMCST.h" … … 50 52 REAL, DIMENSION(klon,klev) :: pu,pv,pt 51 53 REAL, DIMENSION(klon,klev) :: d_t_diss 52 INTEGER okiophys53 54 54 55 REAL timestep … … 68 69 REAL unsdzdec(klon,klev+1) 69 70 70 REAL km(klon,klev +1)71 REAL km(klon,klev) 71 72 REAL kmpre(klon,klev+1),tmp2 72 73 REAL mpre(klon,klev+1) 73 REAL kn(klon,klev +1)74 REAL kq(klon,klev +1)74 REAL kn(klon,klev) 75 REAL kq(klon,klev) 75 76 real ff(klon,klev+1),delta(klon,klev+1) 76 77 real aa(klon,klev+1),aa0,aa1 … … 84 85 data first,ipas/.false.,0/ 85 86 !$OMP THREADPRIVATE( first,ipas) 87 INTEGER, SAVE :: iflag_tke_diff=0 88 !$OMP THREADPRIVATE(iflag_tke_diff) 89 86 90 87 91 integer ig,k … … 119 123 REAL, DIMENSION(klon,klev) :: exner,masse 120 124 REAL, DIMENSION(klon,klev+1) :: masseb,q2old,q2neg 125 LOGICAL okiophys 121 126 122 127 frif(ri)=0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156)) … … 128 133 129 134 135 okiophys=klon==1 130 136 if (firstcall) then 137 CALL getin_p('iflag_tke_diff',iflag_tke_diff) 131 138 allocate(l0(klon)) 132 #ifdef IOPHYS 133 call iophys_ini 139 #define IOPHYS 140 #ifdef IOPHYS 141 ! call iophys_ini 134 142 #endif 135 143 firstcall=.false. 136 144 endif 137 145 138 139 #ifdef IOPHYS 140 if (okiophys==1) then 146 IF (ngrid<=0) RETURN ! Bizarre : on n a pas ce probeleme pour coef_diff_turb 147 148 #ifdef IOPHYS 149 if (okiophys) then 141 150 call iophys_ecrit('q2i',klev,'q2 debut my','m2/s2',q2(:,1:klev)) 142 151 call iophys_ecrit('kmi',klev,'Kz debut my','m/s2',km(:,1:klev)) … … 146 155 nlay=klev 147 156 nlev=klev+1 157 148 158 149 159 !------------------------------------------------------------------------- … … 152 162 153 163 154 zu(:,:)=pu(:,:)+0.5*d_u(:,:) 155 zv(:,:)=pv(:,:)+0.5*d_v(:,:) 156 zt(:,:)=pt(:,:)+0.5*d_t(:,:) 164 zalpha=0.5 ! Anciennement 0.5. Essayer de voir pourquoi ? 165 zu(:,:)=pu(:,:)+zalpha*d_u(:,:) 166 zv(:,:)=pv(:,:)+zalpha*d_v(:,:) 167 zt(:,:)=pt(:,:)+zalpha*d_t(:,:) 168 157 169 do k=1,klev 158 170 exner(:,k)=(play(:,k)/plev(:,1))**RKAPPA 159 171 masse(:,k)=(plev(:,k)-plev(:,k+1))/RG 172 teta(:,k)=zt(:,k)/exner(:,k) 160 173 enddo 161 teta(:,:)=zt(:,:)/exner(:,:)162 174 163 175 ! Atmospheric mass at layer interfaces, where the TKE is computed … … 168 180 enddo 169 181 masseb(:,:)=0.5*masseb(:,:) 170 171 172 182 173 183 zlev(:,1)=0. … … 202 212 203 213 #ifdef IOPHYS 204 if (okiophys ==1) then214 if (okiophys) then 205 215 call iophys_ecrit('zlay',klev,'Geop','m',zlay) 206 216 call iophys_ecrit('teta',klev,'teta','K',teta) 207 217 call iophys_ecrit('temp',klev,'temp','K',zt) 208 218 call iophys_ecrit('pt',klev,'temp','K',pt) 219 call iophys_ecrit('pu',klev,'u','m/s',pu) 220 call iophys_ecrit('pv',klev,'v','m/s',pv) 209 221 call iophys_ecrit('d_u',klev,'d_u','m/s2',d_u) 210 222 call iophys_ecrit('d_v',klev,'d_v','m/s2',d_v) … … 213 225 call iophys_ecrit('masse',klev,'masse','',masse) 214 226 call iophys_ecrit('masseb',klev,'masseb','',masseb) 215 call iophys_ecrit('Cm2',klev,'m2 conserv','m/s',(dddu(:,1:klev)+dddv(:,1:klev))/(masseb(:,1:klev)*timestep))216 call iophys_ecrit('Cn2',klev,'m2 conserv','m/s',(rcpd*dddt(:,1:klev)/masseb(:,1:klev))/timestep)217 call iophys_ecrit('rifc',klev,'rif conservative','',rcpd*dddt(:,1:klev)/min(dddu(:,1:klev)+dddv(:,1:klev),-1.e-20))218 227 endif 219 228 #endif … … 273 282 rif(ig,k)=rifc 274 283 endif 275 if(rif(ig,k) .lt.0.16) then284 if(rif(ig,k)<0.16) then 276 285 alpha(ig,k)=falpha(rif(ig,k)) 277 286 sm(ig,k)=fsm(rif(ig,k)) … … 338 347 339 348 #ifdef IOPHYS 340 if (okiophys ==1) then349 if (okiophys) then 341 350 call iophys_ecrit('rif',klev,'Flux Richardson','m',rif(:,1:klev)) 342 351 call iophys_ecrit('m2',klev,'m2 ','m/s',m2(:,1:klev)) 343 call iophys_ecrit('Km2 ',klev,'m2 conserv','m/s',km(:,1:klev)*m2(:,1:klev))352 call iophys_ecrit('Km2app',klev,'m2 conserv','m/s',km(:,1:klev)*m2(:,1:klev)) 344 353 call iophys_ecrit('Km',klev,'Km','m2/s',km(:,1:klev)) 345 354 endif … … 357 366 ! Evolution of TKE under source terms K M2 and K N2 358 367 leff(:,:)=max(l(:,:),1.) 359 IF (iflag_pbl==29) THEN 360 km2(:,:)=km(:,:)*m2(:,:) 361 kn2(:,:)=kn2(:,:)*rif(:,:) 362 ELSEIF (iflag_pbl==25) THEN 363 DO k=1,klev 364 km2(:,k)=-0.5*(dddu(:,k)+dddv(:,k)+dddu(:,k+1)+dddv(:,k+1)) & 365 & /(masse(:,k)*timestep) 366 kn2(:,k)=rcpd*0.5*(dddt(:,k)+dddt(:,k+1))/(masse(:,k)*timestep) 367 leff(:,k)=0.5*(leff(:,k)+leff(:,k+1)) 368 ENDDO 369 km2(:,klev+1)=0. ; kn2(:,klev+1)=0. 370 ELSE 368 369 !################################################################## 370 !# IF (iflag_pbl==29) THEN 371 !# STOP'Ne pas utiliser iflag_pbl=29' 372 !# km2(:,:)=km(:,:)*m2(:,:) 373 !# kn2(:,:)=kn2(:,:)*rif(:,:) 374 !# ELSEIF (iflag_pbl==25) THEN 375 ! VERSION AVEC LA TKE EN MILIEU DE COUCHE 376 !# STOP'Ne pas utiliser iflag_pbl=25' 377 !# DO k=1,klev 378 !# km2(:,k)=-0.5*(dddu(:,k)+dddv(:,k)+dddu(:,k+1)+dddv(:,k+1)) & 379 !# & /(masse(:,k)*timestep) 380 !# kn2(:,k)=rcpd*0.5*(dddt(:,k)+dddt(:,k+1))/(masse(:,k)*timestep) 381 !# leff(:,k)=0.5*(leff(:,k)+leff(:,k+1)) 382 !# ENDDO 383 !# km2(:,klev+1)=0. ; kn2(:,klev+1)=0. 384 !# ELSE 385 !################################################################# 386 371 387 km2(:,:)=-(dddu(:,:)+dddv(:,:))/(masseb(:,:)*timestep) 372 388 kn2(:,:)=rcpd*dddt(:,:)/(masseb(:,:)*timestep) 373 ENDIF389 ! ENDIF 374 390 q2neg(:,:)=q2(:,:)+timestep*(km2(:,:)-kn2(:,:)) 375 391 q2(:,:)=min(max(q2neg(:,:),1.e-10),1.e4) 392 393 394 #ifdef IOPHYS 395 if (okiophys) then 396 call iophys_ecrit('km2',klev,'m2 conserv','m/s',km2(:,1:klev)) 397 call iophys_ecrit('kn2',klev,'n2 conserv','m/s',kn2(:,1:klev)) 398 endif 399 #endif 376 400 377 401 ! Dissipation of TKE … … 379 403 q2(:,:)=1./(1./sqrt(q2(:,:))+timestep/(2*leff(:,:)*b1)) 380 404 q2(:,:)=q2(:,:)*q2(:,:) 381 405 ! IF (iflag_pbl<=24) THEN 382 406 DO k=1,klev 383 407 d_t_diss(:,k)=(masseb(:,k)*(q2neg(:,k)-q2(:,k))+masseb(:,k+1)*(q2neg(:,k+1)-q2(:,k+1)))/(2.*rcpd*masse(:,k)) 384 408 ENDDO 385 ELSE IF (iflag_pbl<=27) THEN 386 DO k=1,klev 387 d_t_diss(:,k)=(q2neg(:,k)-q2(:,k))/rcpd 388 ENDDO 389 ENDIF 409 410 !################################################################### 411 ! ELSE IF (iflag_pbl<=27) THEN 412 ! DO k=1,klev 413 ! d_t_diss(:,k)=(q2neg(:,k)-q2(:,k))/rcpd 414 ! ENDDO 415 ! ENDIF 390 416 ! print*,'iflag_pbl ',d_t_diss 417 !################################################################### 391 418 392 419 393 420 ! Compuation of stability functions 394 IF (iflag_pbl/=29) THEN421 ! IF (iflag_pbl/=29) THEN 395 422 DO k=1,klev 396 423 DO ig=1,ngrid … … 409 436 ENDDO 410 437 ENDDO 411 ENDIF438 ! ENDIF 412 439 413 440 ! Computation of turbulent diffusivities 414 IF (25<=iflag_pbl.and.iflag_pbl<=28) THEN 415 DO k=2,klev 416 sqrtq(:,k)=sqrt(0.5*(q2(:,k)+q2(:,k-1))) 417 ENDDO 418 ELSE 419 DO k=2,klev 420 sqrtq(:,k)=sqrt(q2(:,k)) 421 ENDDO 422 ENDIF 423 DO k=2,klev 424 DO ig=1,ngrid 425 km(ig,k)=leff(ig,k)*sqrtq(ig,k)*sm(ig,k) 426 kn(ig,k)=km(ig,k)*alpha(ig,k) 427 kq(ig,k)=leff(ig,k)*zq*0.2 428 ! print*,q2(ig,k),zq,km(ig,k) 441 ! IF (25<=iflag_pbl.and.iflag_pbl<=28) THEN 442 ! DO k=2,klev 443 ! sqrtq(:,k)=sqrt(0.5*(q2(:,k)+q2(:,k-1))) 444 ! ENDDO 445 ! ELSE 446 kq(:,:)=0. 447 DO k=1,klev 448 ! Coefficient au milieu des couches pour diffuser la TKE 449 kq(:,k)=0.5*leff(:,k)*sqrt(q2(:,k))*0.2 429 450 ENDDO 451 452 #ifdef IOPHYS 453 if (okiophys) then 454 call iophys_ecrit('q2b',klev,'KTE inter','m2/s',q2(:,1:klev)) 455 endif 456 #endif 457 458 IF (iflag_tke_diff==1) THEN 459 CALL vdif_q2(timestep, RG, RD, ngrid, plev, pt, kq, q2) 460 ENDIF 461 462 km(:,:)=0. 463 kn(:,:)=0. 464 DO k=1,klev 465 km(:,k)=leff(:,k)*sqrt(q2(:,k))*sm(:,k) 466 kn(:,k)=km(:,k)*alpha(:,k) 430 467 ENDDO 431 468 432 469 433 434 #ifdef IOPHYS 435 if (okiophys==1) then 470 #ifdef IOPHYS 471 if (okiophys) then 436 472 call iophys_ecrit('mixingl',klev,'Mixing length','m',leff(:,1:klev)) 437 473 call iophys_ecrit('rife',klev,'Flux Richardson','m',rif(:,1:klev)) … … 447 483 #endif 448 484 485 449 486 ENDIF 450 487 … … 452 489 ! print*,'OK2' 453 490 RETURN 454 !==================================================================== 455 ! Yamada 2.0 456 !==================================================================== 457 if (iflag_pbl.eq.6) then 458 459 do k=2,klev 460 q2(:,k)=l(:,k)**2*zz(:,k) 461 enddo 462 463 464 else if (iflag_pbl.eq.7) then 465 !==================================================================== 466 ! Yamada 2.Fournier 467 !==================================================================== 468 469 ! Calcul de l, km, au pas precedent 470 do k=2,klev 471 do ig=1,ngrid 472 ! print*,'SMML=',sm(ig,k),l(ig,k) 473 delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k)) 474 kmpre(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k) 475 mpre(ig,k)=sqrt(m2(ig,k)) 476 ! print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k) 477 enddo 478 enddo 479 480 do k=2,klev-1 481 do ig=1,ngrid 482 m2cstat=max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1,1.e-12) 483 mcstat=sqrt(m2cstat) 484 485 ! print*,'M2 L=',k,mpre(ig,k),mcstat 486 ! 487 ! -----{puis on ecrit la valeur de q qui annule l'equation de m 488 ! supposee en q3} 489 ! 490 IF (k.eq.2) THEN 491 kmcstat=1.E+0 / mcstat & 492 & *( unsdz(ig,k)*kmpre(ig,k+1) & 493 & *mpre(ig,k+1) & 494 & +unsdz(ig,k-1) & 495 & *cd(ig) & 496 & *( sqrt(zu(ig,3)**2+zv(ig,3)**2) & 497 & -mcstat/unsdzdec(ig,k) & 498 & -mpre(ig,k+1)/unsdzdec(ig,k+1) )**2) & 499 & /( unsdz(ig,k)+unsdz(ig,k-1) ) 500 ELSE 501 kmcstat=1.E+0 / mcstat & 502 & *( unsdz(ig,k)*kmpre(ig,k+1) & 503 & *mpre(ig,k+1) & 504 & +unsdz(ig,k-1)*kmpre(ig,k-1) & 505 & *mpre(ig,k-1) ) & 506 & /( unsdz(ig,k)+unsdz(ig,k-1) ) 507 ENDIF 508 ! print*,'T2 L=',k,tmp2 509 tmp2=kmcstat & 510 & /( sm(ig,k)/q2(ig,k) ) & 511 & /l(ig,k) 512 q2(ig,k)=max(tmp2,1.e-12)**(2./3.) 513 ! print*,'Q2 L=',k,q2(ig,k) 514 ! 515 enddo 516 enddo 517 518 else if (iflag_pbl==8.or.iflag_pbl==9) then 519 !==================================================================== 520 ! Yamada 2.5 a la Didi 521 !==================================================================== 522 523 524 ! Calcul de l, km, au pas precedent 525 do k=2,klev 526 do ig=1,ngrid 527 ! print*,'SMML=',sm(ig,k),l(ig,k) 528 delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k)) 529 if (delta(ig,k).lt.1.e-20) then 530 ! print*,'ATTENTION L=',k,' Delta=',delta(ig,k) 531 delta(ig,k)=1.e-20 532 endif 533 km(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k) 534 aa0=(m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1) 535 aa1=(m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1) 536 ! abder print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20) 537 aa(ig,k)=aa1*timestep/(delta(ig,k)*l(ig,k)) 538 ! print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k) 539 qpre=sqrt(q2(ig,k)) 540 ! if (iflag_pbl.eq.8 ) then 541 if (aa(ig,k).gt.0.) then 542 q2(ig,k)=(qpre+aa(ig,k)*qpre*qpre)**2 543 else 544 q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2 545 endif 546 ! else ! iflag_pbl=9 547 ! if (aa(ig,k)*qpre.gt.0.9) then 548 ! q2(ig,k)=(qpre*10.)**2 549 ! else 550 ! q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2 551 ! endif 552 ! endif 553 q2(ig,k)=min(max(q2(ig,k),1.e-10),1.e4) 554 ! print*,'Q2 L=',k,q2(ig,k),qpre*qpre 555 enddo 556 enddo 557 558 else if (iflag_pbl>=10) then 559 560 ! print*,'Schema mixte D' 561 ! print*,'Longueur ',l(:,:) 562 do k=2,klev-1 563 l(:,k)=max(l(:,k),1.) 564 km(:,k)=l(:,k)*sqrt(q2(:,k))*sm(:,k) 565 q2(:,k)=q2(:,k)+timestep*km(:,k)*m2(:,k)*(1.-rif(:,k)) 566 q2(:,k)=min(max(q2(:,k),1.e-10),1.e4) 567 q2(:,k)=1./(1./sqrt(q2(:,k))+timestep/(2*l(:,k)*b1)) 568 q2(:,k)=q2(:,k)*q2(:,k) 569 enddo 570 571 572 else 573 CALL abort_physic(modname,'Cas nom prevu dans yamada4',1) 574 575 endif ! Fin du cas 8 576 577 ! print*,'OK8' 578 579 !==================================================================== 580 ! Calcul des coefficients de m�ange 581 !==================================================================== 582 do k=2,klev 583 ! print*,'k=',k 584 do ig=1,ngrid 585 !abde print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k) 586 zq=sqrt(q2(ig,k)) 587 km(ig,k)=l(ig,k)*zq*sm(ig,k) 588 kn(ig,k)=km(ig,k)*alpha(ig,k) 589 kq(ig,k)=l(ig,k)*zq*0.2 590 ! print*,'KML=',km(ig,k),kn(ig,k) 591 enddo 592 enddo 593 594 ! Transport diffusif vertical de la TKE. 595 if (iflag_pbl.ge.12) then 596 ! print*,'YAMADA VDIF' 597 q2(:,1)=q2(:,2) 598 call vdif_q2(timestep,RG,RD,ngrid,plev,zt,kq,q2) 599 endif 600 601 ! Traitement des cas noctrunes avec l'introduction d'une longueur 602 ! minilale. 603 604 !==================================================================== 605 ! Traitement particulier pour les cas tres stables. 606 ! D'apres Holtslag Boville. 607 608 if (prt_level>1) THEN 609 print*,'YAMADA4 0' 610 endif !(prt_level>1) THEN 611 do ig=1,ngrid 612 coriol(ig)=1.e-4 613 pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5) 614 enddo 615 616 ! print*,'pblhmin ',pblhmin 617 !Test a remettre 21 11 02 618 ! test abd 13 05 02 if(0.eq.1) then 619 if(1==1) then 620 if(iflag_pbl==8.or.iflag_pbl==10) then 621 622 do k=2,klev 623 do ig=1,ngrid 624 if (teta(ig,2).gt.teta(ig,1)) then 625 qmin=ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2 626 kmin=kap*zlev(ig,k)*qmin 627 else 628 kmin=-1. ! kmin n'est utilise que pour les SL stables. 629 endif 630 if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then 631 ! print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k) 632 ! s ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k) 633 kn(ig,k)=kmin 634 km(ig,k)=kmin 635 kq(ig,k)=kmin 636 ! la longueur de melange est suposee etre l= kap z 637 ! K=l q Sm d'ou q2=(K/l Sm)**2 638 q2(ig,k)=(qmin/sm(ig,k))**2 639 endif 640 enddo 641 enddo 642 643 else 644 645 do k=2,klev 646 do ig=1,ngrid 647 if (teta(ig,2).gt.teta(ig,1)) then 648 qmin=ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2 649 kmin=kap*zlev(ig,k)*qmin 650 else 651 kmin=-1. ! kmin n'est utilise que pour les SL stables. 652 endif 653 if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then 654 ! print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k) 655 ! s ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k) 656 kn(ig,k)=kmin 657 km(ig,k)=kmin 658 kq(ig,k)=kmin 659 ! la longueur de melange est suposee etre l= kap z 660 ! K=l q Sm d'ou q2=(K/l Sm)**2 661 sm(ig,k)=1. 662 alpha(ig,k)=1. 663 q2(ig,k)=min((qmin/sm(ig,k))**2,10.) 664 zq=sqrt(q2(ig,k)) 665 km(ig,k)=l(ig,k)*zq*sm(ig,k) 666 kn(ig,k)=km(ig,k)*alpha(ig,k) 667 kq(ig,k)=l(ig,k)*zq*0.2 668 endif 669 enddo 670 enddo 671 endif 672 673 endif 674 675 if (prt_level>1) THEN 676 print*,'YAMADA4 1' 677 endif !(prt_level>1) THEN 678 ! Diagnostique pour stokage 679 680 if(1.eq.0)then 681 rino=rif 682 smyam(1:ngrid,1)=0. 683 styam(1:ngrid,1)=0. 684 lyam(1:ngrid,1)=0. 685 knyam(1:ngrid,1)=0. 686 w2yam(1:ngrid,1)=0. 687 t2yam(1:ngrid,1)=0. 688 689 smyam(1:ngrid,2:klev)=sm(1:ngrid,2:klev) 690 styam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)*alpha(1:ngrid,2:klev) 691 lyam(1:ngrid,2:klev)=l(1:ngrid,2:klev) 692 knyam(1:ngrid,2:klev)=kn(1:ngrid,2:klev) 693 694 ! Estimations de w'2 et T'2 d'apres Abdela et McFarlane 695 696 w2yam(1:ngrid,2:klev)=q2(1:ngrid,2:klev)*0.24 & 697 & +lyam(1:ngrid,2:klev)*5.17*kn(1:ngrid,2:klev) & 698 & *n2(1:ngrid,2:klev)/sqrt(q2(1:ngrid,2:klev)) 699 700 t2yam(1:ngrid,2:klev)=9.1*kn(1:ngrid,2:klev) & 701 & *dtetadz(1:ngrid,2:klev)**2 & 702 & /sqrt(q2(1:ngrid,2:klev))*lyam(1:ngrid,2:klev) 703 endif 704 705 ! print*,'OKFIN' 706 first=.false. 707 return 708 end 491 END
Note: See TracChangeset
for help on using the changeset viewer.