Changeset 5116 for LMDZ6/branches/Amaury_dev/libf/dyn3d_common
- Timestamp:
- Jul 24, 2024, 2:54:37 PM (4 months ago)
- Location:
- LMDZ6/branches/Amaury_dev/libf/dyn3d_common
- Files:
-
- 40 edited
-
academic.h (modified) (1 diff)
-
adaptdt.f90 (modified) (1 diff)
-
advn.f90 (modified) (42 diffs)
-
advy.f90 (modified) (1 diff)
-
com_io_dyn_mod.F90 (modified) (1 diff)
-
coordij.f90 (modified) (1 diff)
-
diagedyn.f90 (modified) (3 diffs)
-
disvert.F90 (modified) (8 diffs)
-
disvert_noterre.f90 (modified) (17 diffs)
-
exner_hyb_m.F90 (modified) (5 diffs)
-
exner_milieu_m.F90 (modified) (5 diffs)
-
fxhyp_m.F90 (modified) (5 diffs)
-
fyhyp_m.F90 (modified) (2 diffs)
-
gr_int_dyn.f90 (modified) (2 diffs)
-
grilles_gcm_netcdf_sub.F90 (modified) (1 diff)
-
iniconst.F90 (modified) (3 diffs)
-
inidissip.F90 (modified) (8 diffs)
-
inigeom.f90 (modified) (1 diff)
-
inigrads.f90 (modified) (4 diffs)
-
initfluxsto.f90 (modified) (3 diffs)
-
inithist.F90 (modified) (1 diff)
-
inter_barxy_m.F90 (modified) (21 diffs)
-
interpost.f90 (modified) (2 diffs)
-
interpre.f90 (modified) (2 diffs)
-
invert_zoom_x_m.F90 (modified) (2 diffs)
-
iso_verif_dyn.f90 (modified) (4 diffs)
-
limx.f90 (modified) (4 diffs)
-
limy.f90 (modified) (8 diffs)
-
limz.f90 (modified) (4 diffs)
-
pentes_ini.f90 (modified) (15 diffs)
-
ppm3d.f90 (modified) (78 diffs)
-
prather.f90 (modified) (3 diffs)
-
principal_cshift_m.F90 (modified) (2 diffs)
-
sortvarc.f90 (modified) (4 diffs)
-
test_period.f90 (modified) (2 diffs)
-
traceurpole.f90 (modified) (2 diffs)
-
ugeostr.F90 (modified) (2 diffs)
-
write_grads_dyn.h (modified) (2 diffs)
-
writedynav.F90 (modified) (1 diff)
-
writehist.f90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/academic.h
r5099 r5116 3 3 4 4 common/academic/tetarappel,knewt_t,kfrict,knewt_g,clat4 5 real:: tetarappel(ip1jmp1,llm)6 real:: knewt_t(llm)7 real:: kfrict(llm)8 real:: knewt_g9 real:: clat4(ip1jmp1)5 REAL :: tetarappel(ip1jmp1,llm) 6 REAL :: knewt_t(llm) 7 REAL :: kfrict(llm) 8 REAL :: knewt_g 9 REAL :: clat4(ip1jmp1) -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/adaptdt.f90
r5114 r5116 43 43 dtbon=dtvr/n 44 44 45 return46 end subroutineadaptdt45 RETURN 46 END SUBROUTINE adaptdt 47 47 48 48 -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/advn.f90
r5105 r5116 25 25 ! Arguments: 26 26 ! ---------- 27 integer:: mode28 real:: masse(ip1jmp1,llm)27 INTEGER :: mode 28 REAL :: masse(ip1jmp1,llm) 29 29 REAL :: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm) 30 30 REAL :: q(ip1jmp1,llm) … … 35 35 ! 36 36 INTEGER :: i,ij,l,j,ii 37 integer:: ijlqmin,iqmin,jqmin,lqmin38 integer:: ismin39 ! 40 real:: zm(ip1jmp1,llm),newmasse41 real:: mu(ip1jmp1,llm)42 real:: mv(ip1jm,llm)43 real:: mw(ip1jmp1,llm+1)44 real:: zq(ip1jmp1,llm),zz,qpn,qps45 real:: zqg(ip1jmp1,llm),zqd(ip1jmp1,llm)46 real:: zqs(ip1jmp1,llm),zqn(ip1jmp1,llm)47 real:: zqh(ip1jmp1,llm),zqb(ip1jmp1,llm)48 real:: temps0,temps1,temps2,temps349 real:: ztemps1,ztemps2,ssum37 INTEGER :: ijlqmin,iqmin,jqmin,lqmin 38 INTEGER :: ismin 39 ! 40 REAL :: zm(ip1jmp1,llm),newmasse 41 REAL :: mu(ip1jmp1,llm) 42 REAL :: mv(ip1jm,llm) 43 REAL :: mw(ip1jmp1,llm+1) 44 REAL :: zq(ip1jmp1,llm),zz,qpn,qps 45 REAL :: zqg(ip1jmp1,llm),zqd(ip1jmp1,llm) 46 REAL :: zqs(ip1jmp1,llm),zqn(ip1jmp1,llm) 47 REAL :: zqh(ip1jmp1,llm),zqb(ip1jmp1,llm) 48 REAL :: temps0,temps1,temps2,temps3 49 REAL :: ztemps1,ztemps2,ssum 50 50 save temps1,temps2,temps3 51 real:: zzpbar,zzw52 53 real:: qmin,qmax51 REAL :: zzpbar,zzw 52 53 REAL :: qmin,qmax 54 54 data qmin,qmax/0.,1./ 55 55 data temps1,temps2,temps3/0.,0.,0./ … … 139 139 ! Arguments: 140 140 ! ---------- 141 real:: q(ip1jmp1,llm),qg(ip1jmp1,llm),qd(ip1jmp1,llm)141 REAL :: q(ip1jmp1,llm),qg(ip1jmp1,llm),qd(ip1jmp1,llm) 142 142 ! 143 143 ! Local … … 146 146 INTEGER :: ij,l 147 147 ! 148 real:: dxqu(ip1jmp1),zqu(ip1jmp1)149 real:: zqmax(ip1jmp1),zqmin(ip1jmp1)148 REAL :: dxqu(ip1jmp1),zqu(ip1jmp1) 149 REAL :: zqmax(ip1jmp1),zqmin(ip1jmp1) 150 150 logical :: extremum(ip1jmp1) 151 151 152 integer:: mode152 INTEGER :: mode 153 153 save mode 154 154 data mode/1/ … … 156 156 ! calcul des pentes en u: 157 157 ! ----------------------- 158 if (mode==0) then158 if (mode==0) THEN 159 159 do l=1,llm 160 160 do ij=1,ip1jm … … 206 206 enddo 207 207 do ij=iip2+1,ip1jm 208 if(extremum(ij)) then208 IF(extremum(ij)) THEN 209 209 qg(ij,l)=q(ij,l) 210 210 qd(ij,l)=q(ij,l) … … 222 222 223 223 do ij=iip2+1,ip1jm 224 if(extremum(ij).and..not.extremum(ij-1)) &224 IF(extremum(ij).and..not.extremum(ij-1)) & 225 225 qd(ij-1,l)=q(ij,l) 226 226 enddo … … 256 256 ! Arguments: 257 257 ! ---------- 258 real:: q(ip1jmp1,llm),qs(ip1jmp1,llm),qn(ip1jmp1,llm)258 REAL :: q(ip1jmp1,llm),qs(ip1jmp1,llm),qn(ip1jmp1,llm) 259 259 ! 260 260 ! Local … … 263 263 INTEGER :: ij,l 264 264 ! 265 real:: dyqv(ip1jm),zqv(ip1jm,llm)266 real:: zqmax(ip1jm),zqmin(ip1jm)265 REAL :: dyqv(ip1jm),zqv(ip1jm,llm) 266 REAL :: zqmax(ip1jm),zqmin(ip1jm) 267 267 logical :: extremum(ip1jmp1) 268 268 269 integer:: mode269 INTEGER :: mode 270 270 save mode 271 271 data mode/1/ 272 272 273 if (mode==0) then273 if (mode==0) THEN 274 274 do l=1,llm 275 275 do ij=1,ip1jmp1 … … 315 315 316 316 do ij=iip2,ip1jm 317 if(extremum(ij)) then317 IF(extremum(ij)) THEN 318 318 qs(ij,l)=q(ij,l) 319 319 qn(ij,l)=q(ij,l) … … 352 352 ! Arguments: 353 353 ! ---------- 354 real:: q(ip1jmp1,llm),qh(ip1jmp1,llm),qb(ip1jmp1,llm)354 REAL :: q(ip1jmp1,llm),qh(ip1jmp1,llm),qb(ip1jmp1,llm) 355 355 ! 356 356 ! Local … … 359 359 INTEGER :: ij,l 360 360 ! 361 real:: dzqw(ip1jmp1,llm+1),zqw(ip1jmp1,llm+1)362 real:: zqmax(ip1jmp1,llm),zqmin(ip1jmp1,llm)361 REAL :: dzqw(ip1jmp1,llm+1),zqw(ip1jmp1,llm+1) 362 REAL :: zqmax(ip1jmp1,llm),zqmin(ip1jmp1,llm) 363 363 logical :: extremum(ip1jmp1,llm) 364 364 365 integer:: mode365 INTEGER :: mode 366 366 save mode 367 367 … … 371 371 ! ----------------------- 372 372 373 if (mode==0) then373 if (mode==0) THEN 374 374 do l=1,llm 375 375 do ij=1,ip1jmp1 … … 424 424 do l=2,llm-1 425 425 do ij=1,ip1jmp1 426 if(extremum(ij,l)) then426 IF(extremum(ij,l)) THEN 427 427 qh(ij,l)=q(ij,l) 428 428 qb(ij,l)=q(ij,l) … … 435 435 ! do l=2,llm-1 436 436 ! do ij=1,ip1jmp1 437 ! if(extremum(ij,l)) then437 ! IF(extremum(ij,l)) THEN 438 438 ! if (.not.extremum(ij,l-1)) qh(ij,l-1)=q(ij,l) 439 439 ! if (.not.extremum(ij,l+1)) qb(ij,l+1)=q(ij,l) … … 474 474 ! Arguments: 475 475 ! ---------- 476 integer:: mode477 real:: masse(ip1jmp1,llm)478 real:: u_m( ip1jmp1,llm )479 real:: q(ip1jmp1,llm),qd(ip1jmp1,llm),qg(ip1jmp1,llm)476 INTEGER :: mode 477 REAL :: masse(ip1jmp1,llm) 478 REAL :: u_m( ip1jmp1,llm ) 479 REAL :: q(ip1jmp1,llm),qd(ip1jmp1,llm),qg(ip1jmp1,llm) 480 480 ! 481 481 ! Local … … 483 483 ! 484 484 INTEGER :: i,j,ij,l,indu(ip1jmp1),niju,iju,ijq 485 integer:: n0,nl(llm)486 ! 487 real:: new_m,zu_m,zdq,zz488 real:: zsigg(ip1jmp1,llm),zsigd(ip1jmp1,llm),zsig489 real:: u_mq(ip1jmp1,llm)490 491 real:: zm,zq,zsigm,zsigp,zqm,zqp,zu485 INTEGER :: n0,nl(llm) 486 ! 487 REAL :: new_m,zu_m,zdq,zz 488 REAL :: zsigg(ip1jmp1,llm),zsigd(ip1jmp1,llm),zsig 489 REAL :: u_mq(ip1jmp1,llm) 490 491 REAL :: zm,zq,zsigm,zsigp,zqm,zqp,zu 492 492 493 493 logical :: ladvplus(ip1jmp1,llm) 494 494 495 real:: prec495 REAL :: prec 496 496 save prec 497 497 … … 501 501 do ij=iip2,ip1jm 502 502 zdq=qd(ij,l)-qg(ij,l) 503 ! if((qd(ij,l)-q(ij,l))*(q(ij,l)-qg(ij,l)).lt.0.) then503 ! if((qd(ij,l)-q(ij,l))*(q(ij,l)-qg(ij,l)).lt.0.) THEN 504 504 ! PRINT*,'probleme au point ij=',ij,' l=',l 505 505 ! PRINT*,qd(ij,l),q(ij,l),qg(ij,l) … … 507 507 ! qg(ij,l)=q(ij,l) 508 508 ! endif 509 if(abs(zdq)>prec) then509 IF(abs(zdq)>prec) THEN 510 510 zsigd(ij,l)=(q(ij,l)-qg(ij,l))/zdq 511 511 zsigg(ij,l)=1.-zsigd(ij,l) 512 ! if(.not.(zsigd(ij,l).ge.0..and.zsigd(ij,l).le.1. .and.513 ! s zsigg(ij,l).ge.0..or.zsigg(ij,l).le.1.) ) then512 ! IF(.not.(zsigd(ij,l).ge.0..and.zsigd(ij,l).le.1. .and. 513 ! s zsigg(ij,l).ge.0..or.zsigg(ij,l).le.1.) ) THEN 514 514 ! PRINT*,'probleme au point ij=',ij,' l=',l 515 515 ! PRINT*,'sigg=',zsigg(ij,l),' sigd=',zsigd(ij,l) … … 530 530 do l=1,llm 531 531 do ij=iip2,ip1jm-1 532 if (u_m(ij,l)>=0.) then532 if (u_m(ij,l)>=0.) THEN 533 533 zsigp=zsigd(ij,l) 534 534 zsigm=zsigg(ij,l) … … 548 548 ladvplus(ij,l)=zu>zm 549 549 zsig=zu/zm 550 if(zsig==0.) zsigp=0.1551 if (mode==1) then552 if (zsig<=zsigp) then550 IF(zsig==0.) zsigp=0.1 551 if (mode==1) THEN 552 if (zsig<=zsigp) THEN 553 553 u_mq(ij,l)=u_m(ij,l)*zqp 554 else if (mode==1) then554 else if (mode==1) THEN 555 555 u_mq(ij,l)= & 556 556 sign(zm,u_m(ij,l))*(zsigp*zqp+(zsig-zsigp)*zqm) 557 557 endif 558 558 else 559 if (zsig<=zsigp) then559 if (zsig<=zsigp) THEN 560 560 u_mq(ij,l)=u_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq)) 561 561 else … … 565 565 endif 566 566 endif 567 ! if(zsig.lt.0.) then567 ! IF(zsig.lt.0.) THEN 568 568 ! PRINT*,'au point ij=',ij,' l=',l,' sig=',zsig 569 569 ! stop … … 587 587 nl(l)=0 588 588 do ij=iip2,ip1jm 589 if(ladvplus(ij,l)) then589 IF(ladvplus(ij,l)) THEN 590 590 nl(l)=nl(l)+1 591 591 u_mq(ij,l)=0. … … 595 595 enddo 596 596 597 if(n0>1) then597 IF(n0>1) THEN 598 598 IF (prt_level > 9) WRITE(lunout,*) & 599 599 'Nombre de points pour lesquels on advect plus que le' & … … 601 601 602 602 do l=1,llm 603 if(nl(l)>0) then603 IF(nl(l)>0) THEN 604 604 iju=0 605 605 ! indicage des mailles concernees par le traitement special 606 606 do ij=iip2,ip1jm 607 if(ladvplus(ij,l).and.mod(ij,iip1)/=0) then607 IF(ladvplus(ij,l).and.mod(ij,iip1)/=0) THEN 608 608 iju=iju+1 609 609 indu(iju)=ij … … 619 619 zu_m=u_m(ij,l) 620 620 u_mq(ij,l)=0. 621 if(zu_m>0.) then621 IF(zu_m>0.) THEN 622 622 ijq=ij 623 623 i=ijq-(j-1)*iip1 … … 632 632 ! ajout de la maille non completement advectee 633 633 zsig=zu_m/masse(ijq,l) 634 if(zsig<=zsigd(ijq,l)) then634 IF(zsig<=zsigd(ijq,l)) THEN 635 635 u_mq(ij,l)=u_mq(ij,l)+zu_m*(qd(ijq,l) & 636 636 -0.5*zsig/zsigd(ijq,l)*(qd(ijq,l)-q(ijq,l))) … … 639 639 ! goto 8888 640 640 zz=0.5*(zsig-zsigd(ijq,l))/zsigg(ijq,l) 641 if(.not.(zz>0..and.zz<=0.5)) then641 IF(.not.(zz>0..and.zz<=0.5)) THEN 642 642 WRITE(lunout,*)'probleme2 au point ij=',ij, & 643 643 ' l=',l … … 662 662 ! 2eme MODIF SPECIFIQUE 663 663 zsig=-zu_m/masse(ij+1,l) 664 if(zsig<=zsigg(ijq,l)) then664 IF(zsig<=zsigg(ijq,l)) THEN 665 665 u_mq(ij,l)=u_mq(ij,l)+zu_m*(qg(ijq,l) & 666 666 -0.5*zsig/zsigg(ijq,l)*(qg(ijq,l)-q(ijq,l))) … … 669 669 ! goto 9999 670 670 zz=0.5*(zsig-zsigg(ijq,l))/zsigd(ijq,l) 671 if(.not.(zz>0..and.zz<=0.5)) then671 IF(.not.(zz>0..and.zz<=0.5)) THEN 672 672 WRITE(lunout,*)'probleme22 au point ij=',ij & 673 673 ,' l=',l … … 736 736 ! Arguments: 737 737 ! ---------- 738 real:: masse(ip1jmp1,llm)739 real:: v_m( ip1jm,llm )740 real:: q(ip1jmp1,llm),qn(ip1jmp1,llm),qs(ip1jmp1,llm)738 REAL :: masse(ip1jmp1,llm) 739 REAL :: v_m( ip1jm,llm ) 740 REAL :: q(ip1jmp1,llm),qn(ip1jmp1,llm),qs(ip1jmp1,llm) 741 741 ! 742 742 ! Local … … 745 745 INTEGER :: ij,l 746 746 ! 747 real:: new_m,zdq,zz748 real:: zsigs(ip1jmp1),zsign(ip1jmp1),zsig749 real:: v_mq(ip1jm,llm)750 real:: convpn,convps,convmpn,convmps,massen,masses751 real:: zm,zq,zsigm,zsigp,zqm,zqp752 real:: ssum753 real:: prec747 REAL :: new_m,zdq,zz 748 REAL :: zsigs(ip1jmp1),zsign(ip1jmp1),zsig 749 REAL :: v_mq(ip1jm,llm) 750 REAL :: convpn,convps,convmpn,convmps,massen,masses 751 REAL :: zm,zq,zsigm,zsigp,zqm,zqp 752 REAL :: ssum 753 REAL :: prec 754 754 save prec 755 755 … … 758 758 do ij=1,ip1jmp1 759 759 zdq=qn(ij,l)-qs(ij,l) 760 ! if((qn(ij,l)-q(ij,l))*(q(ij,l)-qs(ij,l)).lt.0.) then760 ! if((qn(ij,l)-q(ij,l))*(q(ij,l)-qs(ij,l)).lt.0.) THEN 761 761 ! PRINT*,'probleme au point ij=',ij,' l=',l,' advnqx' 762 762 ! PRINT*,qn(ij,l),q(ij,l),qs(ij,l) … … 764 764 ! qs(ij,l)=q(ij,l) 765 765 ! endif 766 if(abs(zdq)>prec) then766 IF(abs(zdq)>prec) THEN 767 767 zsign(ij)=(q(ij,l)-qs(ij,l))/zdq 768 768 zsigs(ij)=1.-zsign(ij) 769 ! if(.not.(zsign(ij).ge.0..and.zsign(ij).le.1. .and.770 ! s zsigs(ij).ge.0..or.zsigs(ij).le.1.) ) then769 ! IF(.not.(zsign(ij).ge.0..and.zsign(ij).le.1. .and. 770 ! s zsigs(ij).ge.0..or.zsigs(ij).le.1.) ) THEN 771 771 ! PRINT*,'probleme au point ij=',ij,' l=',l 772 772 ! PRINT*,'sigs=',zsigs(ij),' sign=',zsign(ij) … … 782 782 783 783 do ij=1,ip1jm 784 if (v_m(ij,l)>=0.) then784 if (v_m(ij,l)>=0.) THEN 785 785 zsigp=zsign(ij+iip1) 786 786 zsigm=zsigs(ij+iip1) … … 798 798 endif 799 799 zsig=abs(v_m(ij,l))/zm 800 if(zsig==0.) zsigp=0.1801 if (zsig<=zsigp) then800 IF(zsig==0.) zsigp=0.1 801 if (zsig<=zsigp) THEN 802 802 v_mq(ij,l)=v_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq)) 803 803 else … … 863 863 ! Arguments: 864 864 ! ---------- 865 real:: masse(ip1jmp1,llm)866 real:: w_m( ip1jmp1,llm+1)867 real:: q(ip1jmp1,llm),qb(ip1jmp1,llm),qh(ip1jmp1,llm)865 REAL :: masse(ip1jmp1,llm) 866 REAL :: w_m( ip1jmp1,llm+1) 867 REAL :: q(ip1jmp1,llm),qb(ip1jmp1,llm),qh(ip1jmp1,llm) 868 868 869 869 ! … … 873 873 INTEGER :: ij,l 874 874 ! 875 real:: new_m,zdq,zz876 real:: zsigh(ip1jmp1,llm),zsigb(ip1jmp1,llm),zsig877 real:: w_mq(ip1jmp1,llm+1)878 real:: zm,zq,zsigm,zsigp,zqm,zqp879 real:: prec875 REAL :: new_m,zdq,zz 876 REAL :: zsigh(ip1jmp1,llm),zsigb(ip1jmp1,llm),zsig 877 REAL :: w_mq(ip1jmp1,llm+1) 878 REAL :: zm,zq,zsigm,zsigp,zqm,zqp 879 REAL :: prec 880 880 save prec 881 881 … … 885 885 do ij=1,ip1jmp1 886 886 zdq=qb(ij,l)-qh(ij,l) 887 ! if((qh(ij,l)-q(ij,l))*(q(ij,l)-qb(ij,l)).lt.0.) then887 ! if((qh(ij,l)-q(ij,l))*(q(ij,l)-qb(ij,l)).lt.0.) THEN 888 888 ! PRINT*,'probleme au point ij=',ij,' l=',l 889 889 ! PRINT*,qh(ij,l),q(ij,l),qb(ij,l) … … 892 892 ! endif 893 893 894 if(abs(zdq)>prec) then894 IF(abs(zdq)>prec) THEN 895 895 zsigb(ij,l)=(q(ij,l)-qh(ij,l))/zdq 896 896 zsigh(ij,l)=1.-zsigb(ij,l) … … 907 907 do l=2,llm 908 908 do ij=1,ip1jmp1 909 if (w_m(ij,l)>=0.) then909 if (w_m(ij,l)>=0.) THEN 910 910 zsigp=zsigb(ij,l) 911 911 zsigm=zsigh(ij,l) … … 923 923 endif 924 924 zsig=abs(w_m(ij,l))/zm 925 if(zsig==0.) zsigp=0.1926 if (zsig<=zsigp) then925 IF(zsig==0.) zsigp=0.1 926 if (zsig<=zsigp) THEN 927 927 w_mq(ij,l)=w_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq)) 928 928 else -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/advy.f90
r5105 r5116 91 91 REAL :: cy1(llm,ntra), cyLAT(llm,ntra) 92 92 REAL :: z1(iim), zcos(iim), zsin(iim) 93 real:: smpn,smps,s0pn,s0ps93 REAL :: smpn,smps,s0pn,s0ps 94 94 REAL :: SSUM 95 95 EXTERNAL SSUM -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/com_io_dyn_mod.F90
r5113 r5116 8 8 ! Names of various files for outputs (in the dynamics) 9 9 ! to store instantaneous values: 10 character(len=18),parameter :: dynhist_file="dyn_hist.nc" ! on scalar grid11 character(len=18),parameter :: dynhistv_file="dyn_histv.nc" ! on v grid12 character(len=18),parameter :: dynhistu_file="dyn_histu.nc" ! on u grid10 CHARACTER(LEN=18),parameter :: dynhist_file="dyn_hist.nc" ! on scalar grid 11 CHARACTER(LEN=18),parameter :: dynhistv_file="dyn_histv.nc" ! on v grid 12 CHARACTER(LEN=18),parameter :: dynhistu_file="dyn_histu.nc" ! on u grid 13 13 14 14 ! to store averaged values: 15 character(len=18),parameter :: dynhistave_file="dyn_hist_ave.nc"16 character(len=18),parameter :: dynhistvave_file="dyn_histv_ave.nc"17 character(len=18),parameter :: dynhistuave_file="dyn_histu_ave.nc"15 CHARACTER(LEN=18),parameter :: dynhistave_file="dyn_hist_ave.nc" 16 CHARACTER(LEN=18),parameter :: dynhistvave_file="dyn_histv_ave.nc" 17 CHARACTER(LEN=18),parameter :: dynhistuave_file="dyn_histu_ave.nc" 18 18 19 19 ! Ids of various files for outputs (in the dynamics) 20 20 21 21 ! instantaneous (these are set by inithist.F) 22 integer:: histid23 integer:: histvid24 integer:: histuid22 INTEGER :: histid 23 INTEGER :: histvid 24 INTEGER :: histuid 25 25 26 26 ! averages (these are set by initdynav.F) 27 integer:: histaveid28 integer:: histvaveid29 integer:: histuaveid27 INTEGER :: histaveid 28 INTEGER :: histvaveid 29 INTEGER :: histuaveid 30 30 31 31 end module com_io_dyn_mod -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/coordij.f90
r5105 r5116 22 22 include "comgeom.h" 23 23 24 real:: zlon,zlat24 REAL :: zlon,zlat 25 25 26 26 zlon=lon*pi/180. -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/diagedyn.f90
r5105 r5116 138 138 INTEGER :: ndiag ! max number of diagnostic in parallel 139 139 PARAMETER (ndiag=10) 140 integer:: pas(ndiag)140 INTEGER :: pas(ndiag) 141 141 save pas 142 142 data pas/ndiag*0/ … … 160 160 ! 161 161 PRINT*,'MAIS POURQUOI DONC DIAGEDYN NE MARCHE PAS ?' 162 return162 RETURN 163 163 ! On ne garde les donnees que dans les colonnes i=1,iim 164 164 DO jj = 1,jjp1 … … 325 325 ! 326 326 ELSE 327 write(lunout,*)'diagedyn: set to function with Earth parameters'327 WRITE(lunout,*)'diagedyn: set to function with Earth parameters' 328 328 ENDIF ! of if (planet_type=="earth") 329 329 RETURN -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/disvert.F90
r5113 r5116 3 3 SUBROUTINE disvert() 4 4 5 use ioipsl, only: getin6 use new_unit_m, only: new_unit7 use lmdz_assert, only: assert5 use ioipsl, ONLY: getin 6 use new_unit_m, ONLY: new_unit 7 use lmdz_assert, ONLY: assert 8 8 USE comvert_mod, ONLY: ap, bp, aps, bps, nivsigs, nivsig, dpres, presnivs, & 9 9 pseudoalt, pa, preff, scaleheight, presinter … … 47 47 REAL alpha, beta, deltaz 48 48 REAL x 49 character(len=*),parameter :: modname="disvert"50 51 character(len=24):: vert_sampling49 CHARACTER(LEN=*),parameter :: modname="disvert" 50 51 CHARACTER(LEN=24):: vert_sampling 52 52 ! (allowed values are "param", "tropo", "strato" and "read") 53 53 … … 62 62 CALL getin('vert_sampling', vert_sampling) 63 63 WRITE(lunout,*) TRIM(modname)//' vert_sampling = ' // vert_sampling 64 if (llm==39 .and. vert_sampling=="strato") then64 if (llm==39 .and. vert_sampling=="strato") THEN 65 65 dsigmin=0.3 ! Vieille option par défaut pour CMIP5 66 66 else … … 82 82 CLOSE(99) 83 83 alpha=deltaz/(llm*scaleheight) 84 write(lunout, *)trim(modname),':scaleheight, alpha, k0, k1, beta', &84 WRITE(lunout, *)trim(modname),':scaleheight, alpha, k0, k1, beta', & 85 85 scaleheight, alpha, k0, k1, beta 86 86 … … 96 96 dzk1=alpha*tanh(l/k0) 97 97 dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta) 98 write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk298 WRITE(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2 99 99 zkm1=zk 100 100 enddo … … 332 332 ENDDO 333 333 334 write(lunout, *) trim(modname),': BP '335 write(lunout, *) bp336 write(lunout, *) trim(modname),': AP '337 write(lunout, *) ap338 339 write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'340 write(lunout, *)'couches calcules pour une pression de surface =', preff341 write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de '342 write(lunout, *) scaleheight,' km'334 WRITE(lunout, *) trim(modname),': BP ' 335 WRITE(lunout, *) bp 336 WRITE(lunout, *) trim(modname),': AP ' 337 WRITE(lunout, *) ap 338 339 WRITE(lunout, *) 'Niveaux de pressions approximatifs aux centres des' 340 WRITE(lunout, *)'couches calcules pour une pression de surface =', preff 341 WRITE(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de ' 342 WRITE(lunout, *) scaleheight,' km' 343 343 DO l = 1, llm 344 344 dpres(l) = bp(l) - bp(l+1) … … 347 347 presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) 348 348 pseudoalt(l) = log(preff/presnivs(l))*scaleheight 349 write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &349 WRITE(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', & 350 350 pseudoalt(l) & 351 351 , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ & … … 354 354 DO l=1, llmp1 355 355 presinter(l)= ( ap(l)+bp(l)*preff) 356 write(lunout, *)'PRESINTER(', l, ')=', presinter(l)356 WRITE(lunout, *)'PRESINTER(', l, ')=', presinter(l) 357 357 ENDDO 358 358 359 write(lunout, *) trim(modname),': PRESNIVS '360 write(lunout, *) presnivs359 WRITE(lunout, *) trim(modname),': PRESNIVS ' 360 WRITE(lunout, *) presnivs 361 361 362 362 CONTAINS -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/disvert_noterre.f90
r5113 r5116 30 30 REAL :: snorm 31 31 REAL :: alpha,beta,gama,delta,deltaz 32 real:: quoi,quand32 REAL :: quoi,quand 33 33 REAL :: zsig(llm),sig(llm+1) 34 34 INTEGER :: np,ierr 35 integer:: ierr1,ierr2,ierr3,ierr435 INTEGER :: ierr1,ierr2,ierr3,ierr4 36 36 REAL :: x 37 37 38 38 REAL :: SSUM 39 39 EXTERNAL SSUM 40 real:: newsig40 REAL :: newsig 41 41 REAL :: dz0,dz1,nhaut,sig1,esig,csig,zz 42 real:: tt,rr,gg, prevz43 real:: s(llm),dsig(llm)44 45 integer:: iz46 real:: z, ps,p47 character(len=*),parameter :: modname="disvert_noterre"42 REAL :: tt,rr,gg, prevz 43 REAL :: s(llm),dsig(llm) 44 45 INTEGER :: iz 46 REAL :: z, ps,p 47 CHARACTER(LEN=*),parameter :: modname="disvert_noterre" 48 48 49 49 ! … … 55 55 hybrid=.TRUE. ! default value for hybrid (ie: use hybrid coordinates) 56 56 CALL getin('hybrid',hybrid) 57 write(lunout,*) trim(modname),': hybrid=',hybrid57 WRITE(lunout,*) trim(modname),': hybrid=',hybrid 58 58 59 59 ! Ouverture possible de fichiers typiquement E.T. … … 61 61 open(99,file="esasig.def",status='old',form='formatted', & 62 62 iostat=ierr2) 63 if(ierr2/=0) then63 IF(ierr2/=0) THEN 64 64 close(99) 65 65 open(99,file="z2sig.def",status='old',form='formatted', & … … 71 71 ! ---------------------------------------- 72 72 73 IF(ierr2==0) then 74 73 IF(ierr2==0) THEN 75 74 ! Lecture de esasig.def : 76 75 ! Systeme peu souple, mais qui respecte en theorie … … 78 77 ! <-> energie cinetique, d'apres la note de Frederic Hourdin... 79 78 80 write(lunout,*)'*****************************'81 write(lunout,*)'WARNING reading esasig.def'82 write(lunout,*)'*****************************'79 WRITE(lunout,*)'*****************************' 80 WRITE(lunout,*)'WARNING reading esasig.def' 81 WRITE(lunout,*)'*****************************' 83 82 READ(99,*) scaleheight 84 83 READ(99,*) dz0 … … 126 125 ! ---------------------------------------- 127 126 128 ELSE IF(ierr4==0) then129 write(lunout,*)'****************************'130 write(lunout,*)'Reading z2sig.def'131 write(lunout,*)'****************************'127 ELSE IF(ierr4==0) THEN 128 WRITE(lunout,*)'****************************' 129 WRITE(lunout,*)'Reading z2sig.def' 130 WRITE(lunout,*)'****************************' 132 131 133 132 READ(99,*) scaleheight … … 146 145 !----------------------------------------------------------------------- 147 146 ELSE 148 write(lunout,*) 'didn t you forget something ??? '149 write(lunout,*) 'We need file z2sig.def ! (OR esasig.def)'147 WRITE(lunout,*) 'didn t you forget something ??? ' 148 WRITE(lunout,*) 'We need file z2sig.def ! (OR esasig.def)' 150 149 stop 151 150 ENDIF … … 170 169 171 170 if (hybrid) then ! use hybrid coordinates 172 write(lunout,*) "*********************************"173 write(lunout,*) "Using hybrid vertical coordinates"174 write(lunout,*)171 WRITE(lunout,*) "*********************************" 172 WRITE(lunout,*) "Using hybrid vertical coordinates" 173 WRITE(lunout,*) 175 174 ! Coordonnees hybrides avec mod 176 175 DO l = 1, llm … … 183 182 ap(llmp1) = 0. 184 183 else ! use sigma coordinates 185 write(lunout,*) "********************************"186 write(lunout,*) "Using sigma vertical coordinates"187 write(lunout,*)184 WRITE(lunout,*) "********************************" 185 WRITE(lunout,*) "Using sigma vertical coordinates" 186 WRITE(lunout,*) 188 187 ! Pour ne pas passer en coordonnees hybrides 189 188 DO l = 1, llm … … 196 195 bp(llmp1) = 0. 197 196 198 write(lunout,*) trim(modname),': BP '199 write(lunout,*) bp200 write(lunout,*) trim(modname),': AP '201 write(lunout,*) ap197 WRITE(lunout,*) trim(modname),': BP ' 198 WRITE(lunout,*) bp 199 WRITE(lunout,*) trim(modname),': AP ' 200 WRITE(lunout,*) ap 202 201 203 202 ! Calcul au milieu des couches : … … 214 213 ENDDO 215 214 216 if (hybrid) then215 if (hybrid) THEN 217 216 aps(llm) = aps(llm-1)**2 / aps(llm-2) 218 217 bps(llm) = 0.5*(bp(llm) + bp(llm+1)) … … 222 221 end if 223 222 224 write(lunout,*) trim(modname),': BPs '225 write(lunout,*) bps226 write(lunout,*) trim(modname),': APs'227 write(lunout,*) aps223 WRITE(lunout,*) trim(modname),': BPs ' 224 WRITE(lunout,*) bps 225 WRITE(lunout,*) trim(modname),': APs' 226 WRITE(lunout,*) aps 228 227 229 228 DO l = 1, llm … … 232 231 ENDDO 233 232 234 write(lunout,*)trim(modname),' : PRESNIVS'235 write(lunout,*)presnivs236 write(lunout,*)'Pseudo altitude of Presnivs : (for a scale ', &233 WRITE(lunout,*)trim(modname),' : PRESNIVS' 234 WRITE(lunout,*)presnivs 235 WRITE(lunout,*)'Pseudo altitude of Presnivs : (for a scale ', & 237 236 'height of ',scaleheight,' km)' 238 write(lunout,*)pseudoalt237 WRITE(lunout,*)pseudoalt 239 238 240 239 ! -------------------------------------------------- … … 252 251 ! do l=2,llm 253 252 ! approximation of scale height for Venus 254 ! if (zsig(l-1).le.55.) then253 ! if (zsig(l-1).le.55.) THEN 255 254 ! scaleheight = 15.5 - zsig(l-1)/55.*10. 256 255 ! else … … 260 259 ! . log((aps(l) + bps(l)*ps)/(aps(l-1) + bps(l-1)*ps)) 261 260 ! END DO 262 ! write(53,'(I3,50F10.5)') iz, zsig261 ! WRITE(53,'(I3,50F10.5)') iz, zsig 263 262 ! END DO 264 263 ! close(53) … … 296 295 297 296 IMPLICIT NONE 298 real:: x1, x2, sig,pa,preff, newsig, F299 integer:: j297 REAL :: x1, x2, sig,pa,preff, newsig, F 298 INTEGER :: j 300 299 301 300 newsig = sig 302 301 x1=0 303 302 x2=1 304 if (sig>=1) then303 if (sig>=1) THEN 305 304 newsig= sig 306 else if (sig*preff/pa>=0.25) then305 else if (sig*preff/pa>=0.25) THEN 307 306 DO J=1,9999 ! nombre d''iteration max 308 307 F=((1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig)/sig 309 ! write(0,*) J, ' newsig =', newsig, ' F= ', F310 if (F>1) then308 ! WRITE(0,*) J, ' newsig =', newsig, ' F= ', F 309 if (F>1) THEN 311 310 X2 = newsig 312 311 newsig=(X1+newsig)*0.5 … … 319 318 IF(abs(10.*log(F))<1.E-5) goto 999 320 319 END DO 321 else ! if (sig*preff/pa.le.0.25) then320 else ! if (sig*preff/pa.le.0.25) THEN 322 321 newsig= sig*preff/pa 323 322 end if -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/exner_hyb_m.F90
r5106 r5116 54 54 55 55 logical,save :: firstcall=.TRUE. 56 character(len=*),parameter :: modname="exner_hyb"56 CHARACTER(LEN=*),parameter :: modname="exner_hyb" 57 57 58 58 ! Sanity check 59 if (firstcall) then59 if (firstcall) THEN 60 60 ! sanity checks for Shallow Water case (1 vertical layer) 61 if (llm==1) then62 if (kappa/=1) then61 if (llm==1) THEN 62 if (kappa/=1) THEN 63 63 CALL abort_gcm(modname, & 64 64 "kappa!=1 , but running in Shallow Water mode!!",42) 65 65 endif 66 if (cpp/=r) then66 if (cpp/=r) THEN 67 67 CALL abort_gcm(modname, & 68 68 "cpp!=r , but running in Shallow Water mode!!",42) … … 74 74 75 75 ! Specific behaviour for Shallow Water (1 vertical layer) case: 76 if (llm==1) then 77 76 if (llm==1) THEN 78 77 ! Compute pks(:),pk(:),pkf(:) 79 78 … … 83 82 ENDDO 84 83 85 if (present(pkf)) then84 if (present(pkf)) THEN 86 85 pkf = pk 87 86 CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) … … 89 88 90 89 ! our work is done, exit routine 91 return90 RETURN 92 91 endif ! of if (llm.eq.1) 93 92 … … 138 137 ENDDO 139 138 140 if (present(pkf)) then139 if (present(pkf)) THEN 141 140 ! calcul de pkf 142 141 pkf = pk -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/exner_milieu_m.F90
r5106 r5116 51 51 52 52 logical,save :: firstcall=.TRUE. 53 character(len=*),parameter :: modname="exner_milieu"53 CHARACTER(LEN=*),parameter :: modname="exner_milieu" 54 54 55 55 ! Sanity check 56 if (firstcall) then56 if (firstcall) THEN 57 57 ! sanity checks for Shallow Water case (1 vertical layer) 58 if (llm==1) then59 if (kappa/=1) then58 if (llm==1) THEN 59 if (kappa/=1) THEN 60 60 CALL abort_gcm(modname, & 61 61 "kappa!=1 , but running in Shallow Water mode!!",42) 62 62 endif 63 if (cpp/=r) then63 if (cpp/=r) THEN 64 64 CALL abort_gcm(modname, & 65 65 "cpp!=r , but running in Shallow Water mode!!",42) … … 71 71 72 72 ! Specific behaviour for Shallow Water (1 vertical layer) case: 73 if (llm==1) then 74 73 if (llm==1) THEN 75 74 ! Compute pks(:),pk(:),pkf(:) 76 75 … … 80 79 ENDDO 81 80 82 if (present(pkf)) then81 if (present(pkf)) THEN 83 82 pkf = pk 84 83 CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) … … 86 85 87 86 ! our work is done, exit routine 88 return87 RETURN 89 88 endif ! of if (llm.eq.1) 90 89 … … 117 116 ENDDO 118 117 119 if (present(pkf)) then118 if (present(pkf)) THEN 120 119 ! calcul de pkf 121 120 pkf = pk -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/fxhyp_m.F90
r5113 r5116 18 18 ! 1., taux=0., clon=0.) est à - 180 degrés. 19 19 20 use lmdz_arth, only: arth21 use invert_zoom_x_m, only: invert_zoom_x, nmax22 use nrtype, only: pi, pi_d, twopi, twopi_d, k823 use principal_cshift_m, only: principal_cshift24 use serre_mod, only: clon, grossismx, dzoomx, taux20 use lmdz_arth, ONLY: arth 21 use invert_zoom_x_m, ONLY: invert_zoom_x, nmax 22 use nrtype, ONLY: pi, pi_d, twopi, twopi_d, k8 23 use principal_cshift_m, ONLY: principal_cshift 24 use serre_mod, ONLY: clon, grossismx, dzoomx, taux 25 25 26 26 include "dimensions.h" … … 45 45 print *, "Call sequence information: fxhyp" 46 46 47 test_iim: if (iim==1) then47 test_iim: if (iim==1) THEN 48 48 rlonv(1)=0. 49 49 rlonu(1)=pi … … 56 56 xprimp025(:)=1. 57 57 else test_iim 58 test_grossismx: if (grossismx == 1.) then58 test_grossismx: if (grossismx == 1.) THEN 59 59 step = twopi / iim 60 60 … … 193 193 END DO 194 194 195 if (rlonm025(is2) < - pi) then195 if (rlonm025(is2) < - pi) THEN 196 196 print *, 'Rlonm025 plus petit que - pi !' 197 197 STOP 1 … … 204 204 END DO 205 205 206 if (rlonm025(is2) > pi) then206 if (rlonm025(is2) > pi) THEN 207 207 print *, 'Rlonm025 plus grand que pi !' 208 208 STOP 1 -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/fyhyp_m.F90
r5114 r5116 16 16 ! Il vaut mieux avoir : grossismy * dzoom < pi / 2 17 17 18 use lmdz_coefpoly, only: coefpoly19 use nrtype, only: k820 use serre_mod, only: clat, grossismy, dzoomy, tauy18 use lmdz_coefpoly, ONLY: coefpoly 19 use nrtype, ONLY: k8 20 use serre_mod, ONLY: clat, grossismy, dzoomy, tauy 21 21 22 22 include "dimensions.h" … … 207 207 yprimin = a1 + 2.*a2*yi + 3.*a3*yi2 208 208 END DO 209 if (abs(yi-yo1) > epsilon) then209 if (abs(yi-yo1) > epsilon) THEN 210 210 print *, 'Pas de solution.', j, ylon2 211 211 STOP 1 -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/gr_int_dyn.f90
r5113 r5116 12 12 13 13 INTEGER :: iim 14 integer:: ip1, jp114 INTEGER :: ip1, jp1 15 15 REAL :: champin(iim, jp1) 16 16 REAL :: champdyn(iim+1, jp1) 17 17 18 18 INTEGER :: i, j 19 real:: polenord, polesud19 REAL :: polenord, polesud 20 20 21 21 !----------------------------------------------------------------------- … … 34 34 do j = 1, jp1 35 35 do i = 1, iim 36 if (j == 1) then36 if (j == 1) THEN 37 37 champdyn(i, j) = polenord 38 else if (j == jp1) then38 else if (j == jp1) THEN 39 39 champdyn(i, j) = polesud 40 40 else -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/grilles_gcm_netcdf_sub.F90
r5103 r5116 226 226 ! CLOSE netcdf file 227 227 CALL ncclos(ncid_out,rcode_out) 228 write(*,*) "END grilles_gcm_netcdf_sub OK"228 WRITE(*,*) "END grilles_gcm_netcdf_sub OK" 229 229 230 230 END SUBROUTINE grilles_gcm_netcdf_sub -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/iniconst.F90
r5103 r5116 21 21 include "iniprint.h" 22 22 23 character(len=*),parameter :: modname="iniconst"24 character(len=80) :: abort_message23 CHARACTER(LEN=*),parameter :: modname="iniconst" 24 CHARACTER(LEN=80) :: abort_message 25 25 26 26 … … 48 48 r = cpp * kappa 49 49 50 write(lunout,*) trim(modname),': R CP Kappa ',r,cpp,kappa50 WRITE(lunout,*) trim(modname),': R CP Kappa ',r,cpp,kappa 51 51 52 52 !----------------------------------------------------------------------- 53 53 54 54 ! vertical discretization: default behavior depends on planet_type flag 55 if (planet_type=="earth") then55 if (planet_type=="earth") THEN 56 56 disvert_type=1 57 57 else … … 60 60 ! but user can also specify using one or the other in run.def: 61 61 CALL getin('disvert_type',disvert_type) 62 write(lunout,*) trim(modname),': disvert_type=',disvert_type62 WRITE(lunout,*) trim(modname),': disvert_type=',disvert_type 63 63 64 64 pressure_exner = disvert_type == 1 ! default value 65 65 CALL getin('pressure_exner', pressure_exner) 66 66 67 if (disvert_type==1) then67 if (disvert_type==1) THEN 68 68 ! standard case for Earth (automatic generation of levels) 69 69 CALL disvert() 70 else if (disvert_type==2) then70 else if (disvert_type==2) THEN 71 71 ! standard case for planets (levels generated using z2sig.def file) 72 72 CALL disvert_noterre 73 73 else 74 write(abort_message,*) "Wrong value for disvert_type: ", disvert_type74 WRITE(abort_message,*) "Wrong value for disvert_type: ", disvert_type 75 75 CALL abort_gcm(modname,abort_message,0) 76 76 endif -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inidissip.F90
r5106 r5116 1 2 1 ! $Id$ 3 2 4 SUBROUTINE inidissip( lstardis,nitergdiv,nitergrot,niterh, &5 tetagdiv,tetagrot,tetatemp, vert_prof_dissip)3 SUBROUTINE inidissip(lstardis, nitergdiv, nitergrot, niterh, & 4 tetagdiv, tetagrot, tetatemp, vert_prof_dissip) 6 5 !======================================================================= 7 6 ! initialisation de la dissipation horizontale … … 11 10 ! ------------- 12 11 13 USE control_mod, ONLY: dissip_period, iperiod12 USE control_mod, ONLY: dissip_period, iperiod 14 13 USE comconst_mod, ONLY: dissip_deltaz, dissip_factz, dissip_zref, & 15 14 dtdiss, dtvr, rad 16 15 USE comvert_mod, ONLY: preff, presnivs 17 16 USE lmdz_filtreg, ONLY: filtreg 17 USE lmdz_libmath, ONLY: minmax 18 18 19 19 IMPLICIT NONE … … 23 23 include "iniprint.h" 24 24 25 LOGICAL, INTENT(in) :: lstardis26 INTEGER, INTENT(in) :: nitergdiv,nitergrot,niterh27 REAL, INTENT(in) :: tetagdiv,tetagrot,tetatemp28 29 integer, INTENT(in) :: vert_prof_dissip25 LOGICAL, INTENT(in) :: lstardis 26 INTEGER, INTENT(in) :: nitergdiv, nitergrot, niterh 27 REAL, INTENT(in) :: tetagdiv, tetagrot, tetatemp 28 29 integer, INTENT(in) :: vert_prof_dissip 30 30 ! Vertical profile of horizontal dissipation 31 31 ! Allowed values: … … 33 33 ! 1: tanh of altitude 34 34 35 ! Local variables:36 REAL fact, zvert(llm),zz37 REAL zh(ip1jmp1), zu(ip1jmp1), gx(ip1jmp1), divgra(ip1jmp1)38 real zv(ip1jm), gy(ip1jm), deltap(ip1jmp1, llm)39 REAL ullm, vllm,umin,vmin,zhmin,zhmax35 ! Local variables: 36 REAL fact, zvert(llm), zz 37 REAL zh(ip1jmp1), zu(ip1jmp1), gx(ip1jmp1), divgra(ip1jmp1) 38 real zv(ip1jm), gy(ip1jm), deltap(ip1jmp1, llm) 39 REAL ullm, vllm, umin, vmin, zhmin, zhmax 40 40 REAL zllm 41 41 42 INTEGER l, ij,idum,ii42 INTEGER l, ij, idum, ii 43 43 REAL tetamin 44 44 REAL pseudoz 45 character (len =80) :: abort_message45 character (len = 80) :: abort_message 46 46 47 47 REAL ran1 … … 53 53 ! ----------------------------------------------------------------- 54 54 55 crot 56 cdivu 57 cdivh 55 crot = -1. 56 cdivu = -1. 57 cdivh = -1. 58 58 59 59 ! calcul de la valeur propre de divgrad: … … 61 61 idum = 0 62 62 DO l = 1, llm 63 64 deltap(ij,l) = 1.65 66 ENDDO 67 68 idum 69 zh(1) = RAN1(idum) -.570 idum 63 DO ij = 1, ip1jmp1 64 deltap(ij, l) = 1. 65 ENDDO 66 ENDDO 67 68 idum = -1 69 zh(1) = RAN1(idum) - .5 70 idum = 0 71 71 DO ij = 2, ip1jmp1 72 zh(ij) = RAN1(idum) -.573 ENDDO 74 75 CALL filtreg (zh, jjp1,1,2,1,.TRUE.,1)76 77 CALL minmax(iip1 *jjp1,zh,zhmin,zhmax)78 79 IF ( zhmin >= zhmax) THEN80 write(lunout,*)' Inidissip zh min max ',zhmin,zhmax81 abort_message='probleme generateur alleatoire dans inidissip'82 CALL abort_gcm('inidissip',abort_message,1)72 zh(ij) = RAN1(idum) - .5 73 ENDDO 74 75 CALL filtreg (zh, jjp1, 1, 2, 1, .TRUE., 1) 76 77 CALL minmax(iip1 * jjp1, zh, zhmin, zhmax) 78 79 IF (zhmin >= zhmax) THEN 80 WRITE(lunout, *)' Inidissip zh min max ', zhmin, zhmax 81 abort_message = 'probleme generateur alleatoire dans inidissip' 82 CALL abort_gcm('inidissip', abort_message, 1) 83 83 ENDIF 84 84 85 zllm = ABS( zhmax)86 DO l = 1, 5087 88 CALL divgrad2(1,zh,deltap,niterh,divgra)89 90 CALL divgrad (1,zh,niterh,divgra)91 92 93 zllm= ABS(maxval(divgra))94 85 zllm = ABS(zhmax) 86 DO l = 1, 50 87 IF(lstardis) THEN 88 CALL divgrad2(1, zh, deltap, niterh, divgra) 89 ELSE 90 CALL divgrad (1, zh, niterh, divgra) 91 ENDIF 92 93 zllm = ABS(maxval(divgra)) 94 zh = divgra / zllm 95 95 ENDDO 96 96 97 97 IF(lstardis) THEN 98 cdivh = 1./ zllm98 cdivh = 1. / zllm 99 99 ELSE 100 cdivh = zllm ** ( -1./niterh)100 cdivh = zllm ** (-1. / niterh) 101 101 ENDIF 102 102 103 103 ! calcul des valeurs propres de gradiv (ii =1) et nxgrarot(ii=2) 104 104 ! ----------------------------------------------------------------- 105 write(lunout,*)'inidissip: calcul des valeurs propres'105 WRITE(lunout, *)'inidissip: calcul des valeurs propres' 106 106 107 107 DO ii = 1, 2 108 108 109 DO ij = 1, ip1jmp1 110 zu(ij) = RAN1(idum) -.5 111 ENDDO 112 CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1) 113 DO ij = 1, ip1jm 114 zv(ij) = RAN1(idum) -.5 115 ENDDO 116 CALL filtreg (zv,jjm,1,2,1,.FALSE.,1) 117 118 CALL minmax(iip1*jjp1,zu,umin,ullm ) 119 CALL minmax(iip1*jjm, zv,vmin,vllm ) 120 121 ullm = ABS ( ullm ) 122 vllm = ABS ( vllm ) 123 124 DO l = 1, 50 125 IF(ii==1) THEN 126 !cccc CALL covcont( 1,zu,zv,zu,zv ) 127 IF(lstardis) THEN 128 CALL gradiv2( 1,zu,zv,nitergdiv,gx,gy ) 129 ELSE 130 CALL gradiv ( 1,zu,zv,nitergdiv,gx,gy ) 131 ENDIF 109 DO ij = 1, ip1jmp1 110 zu(ij) = RAN1(idum) - .5 111 ENDDO 112 CALL filtreg (zu, jjp1, 1, 2, 1, .TRUE., 1) 113 DO ij = 1, ip1jm 114 zv(ij) = RAN1(idum) - .5 115 ENDDO 116 CALL filtreg (zv, jjm, 1, 2, 1, .FALSE., 1) 117 118 CALL minmax(iip1 * jjp1, zu, umin, ullm) 119 CALL minmax(iip1 * jjm, zv, vmin, vllm) 120 121 ullm = ABS (ullm) 122 vllm = ABS (vllm) 123 124 DO l = 1, 50 125 IF(ii==1) THEN 126 !cccc CALL covcont( 1,zu,zv,zu,zv ) 127 IF(lstardis) THEN 128 CALL gradiv2(1, zu, zv, nitergdiv, gx, gy) 132 129 ELSE 133 IF(lstardis) THEN 134 CALL nxgraro2( 1,zu,zv,nitergrot,gx,gy ) 135 ELSE 136 CALL nxgrarot( 1,zu,zv,nitergrot,gx,gy ) 137 ENDIF 130 CALL gradiv (1, zu, zv, nitergdiv, gx, gy) 138 131 ENDIF 139 140 zllm = max(abs(maxval(gx)), abs(maxval(gy))) 141 zu = gx / zllm 142 zv = gy / zllm 143 end DO 144 145 IF ( ii==1 ) THEN 132 ELSE 146 133 IF(lstardis) THEN 147 cdivu = 1./zllm134 CALL nxgraro2(1, zu, zv, nitergrot, gx, gy) 148 135 ELSE 149 cdivu = zllm **( -1./nitergdiv)136 CALL nxgrarot(1, zu, zv, nitergrot, gx, gy) 150 137 ENDIF 151 ELSE 152 IF(lstardis) THEN 153 crot = 1./ zllm 154 ELSE 155 crot = zllm **( -1./nitergrot ) 156 ENDIF 157 ENDIF 138 ENDIF 139 140 zllm = max(abs(maxval(gx)), abs(maxval(gy))) 141 zu = gx / zllm 142 zv = gy / zllm 143 end DO 144 145 IF (ii==1) THEN 146 IF(lstardis) THEN 147 cdivu = 1. / zllm 148 ELSE 149 cdivu = zllm **(-1. / nitergdiv) 150 ENDIF 151 ELSE 152 IF(lstardis) THEN 153 crot = 1. / zllm 154 ELSE 155 crot = zllm **(-1. / nitergrot) 156 ENDIF 157 ENDIF 158 158 159 159 end DO … … 163 163 164 164 ! IF(.NOT.lstardis) THEN 165 fact = rad*24./REAL(jjm)166 fact = fact*fact167 write(lunout,*)'inidissip: coef u ', fact/cdivu, 1./cdivu168 write(lunout,*)'inidissip: coef r ', fact/crot , 1./crot169 write(lunout,*)'inidissip: coef h ', fact/cdivh, 1./cdivh165 fact = rad * 24. / REAL(jjm) 166 fact = fact * fact 167 WRITE(lunout, *)'inidissip: coef u ', fact / cdivu, 1. / cdivu 168 WRITE(lunout, *)'inidissip: coef r ', fact / crot, 1. / crot 169 WRITE(lunout, *)'inidissip: coef h ', fact / cdivh, 1. / cdivh 170 170 ! ENDIF 171 171 … … 174 174 ! -------------------------------------------------- 175 175 176 if (vert_prof_dissip == 1) then177 do l=1,llm178 pseudoz=8.*log(preff/presnivs(l))179 zvert(l)=1+ &180 (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2. &181 *(dissip_factz-1.)182 176 if (vert_prof_dissip == 1) THEN 177 do l = 1, llm 178 pseudoz = 8. * log(preff / presnivs(l)) 179 zvert(l) = 1 + & 180 (tanh((pseudoz - dissip_zref) / dissip_deltaz) + 1.) / 2. & 181 * (dissip_factz - 1.) 182 enddo 183 183 else 184 DO l=1,llm185 zvert(l)=1.186 187 fact=2.188 189 zz = 1. - preff/presnivs(l)190 zvert(l)= fact -( fact-1.)/( 1.+zz*zz)191 184 DO l = 1, llm 185 zvert(l) = 1. 186 ENDDO 187 fact = 2. 188 DO l = 1, llm 189 zz = 1. - preff / presnivs(l) 190 zvert(l) = fact - (fact - 1.) / (1. + zz * zz) 191 ENDDO 192 192 endif 193 193 194 195 write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale' 196 197 tetamin = 1.e+6 198 199 DO l=1,llm 200 tetaudiv(l) = zvert(l)/tetagdiv 201 tetaurot(l) = zvert(l)/tetagrot 202 tetah(l) = zvert(l)/tetatemp 203 204 IF( tetamin> (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l) 205 IF( tetamin> (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l) 206 IF( tetamin> (1./ tetah(l)) ) tetamin = 1./ tetah(l) 194 WRITE(lunout, *)'inidissip: Constantes de temps de la diffusion horizontale' 195 196 tetamin = 1.e+6 197 198 DO l = 1, llm 199 tetaudiv(l) = zvert(l) / tetagdiv 200 tetaurot(l) = zvert(l) / tetagrot 201 tetah(l) = zvert(l) / tetatemp 202 203 IF(tetamin> (1. / tetaudiv(l))) tetamin = 1. / tetaudiv(l) 204 IF(tetamin> (1. / tetaurot(l))) tetamin = 1. / tetaurot(l) 205 IF(tetamin> (1. / tetah(l))) tetamin = 1. / tetah(l) 207 206 ENDDO 208 207 209 208 ! If dissip_period=0 calculate value for dissipation period, else keep value read from gcm.def 210 209 IF (dissip_period == 0) THEN 211 dissip_period = INT( tetamin/( 2.*dtvr*iperiod)) * iperiod212 write(lunout,*)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ',tetamin,dtvr,iperiod,dissip_period213 dissip_period = MAX(iperiod,dissip_period)210 dissip_period = INT(tetamin / (2. * dtvr * iperiod)) * iperiod 211 WRITE(lunout, *)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ', tetamin, dtvr, iperiod, dissip_period 212 dissip_period = MAX(iperiod, dissip_period) 214 213 END IF 215 214 216 dtdiss 217 write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr218 219 DO l = 1, llm220 write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), &221 dtdiss*tetah(l)215 dtdiss = dissip_period * dtvr 216 WRITE(lunout, *)'inidissip: dissip_period=', dissip_period, ' dtdiss=', dtdiss, ' dtvr=', dtvr 217 218 DO l = 1, llm 219 WRITE(lunout, *)zvert(l), dtdiss * tetaudiv(l), dtdiss * tetaurot(l), & 220 dtdiss * tetah(l) 222 221 ENDDO 223 222 -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inigeom.f90
r5105 r5116 16 16 ! 17 17 ! 18 use fxhyp_m, only: fxhyp19 use fyhyp_m, only: fyhyp18 use fxhyp_m, ONLY: fxhyp 19 use fyhyp_m, ONLY: fyhyp 20 20 USE comconst_mod, ONLY: pi, g, omeg, rad 21 21 USE logic_mod, ONLY: fxyhypb, ysinus -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inigrads.f90
r5113 r5116 7 7 IMPLICIT NONE 8 8 9 integer:: if, im, jm, lm, i, j, l10 real:: x(im), y(jm), z(lm), fx, fy, fz, dt11 real:: xmin, xmax, ymin, ymax9 INTEGER :: if, im, jm, lm, i, j, l 10 REAL :: x(im), y(jm), z(lm), fx, fy, fz, dt 11 REAL :: xmin, xmax, ymin, ymax 12 12 13 character(len= *), intent(in) :: file14 character(len= *), intent(in) :: titlel13 CHARACTER(LEN = *), intent(in) :: file 14 CHARACTER(LEN = *), intent(in) :: titlel 15 15 16 16 INCLUDE "gradsdef.h" 17 17 18 18 ! data unit/66,32,34,36,38,40,42,44,46,48/ 19 integer:: nf19 INTEGER :: nf 20 20 save nf 21 21 data nf/0/ … … 49 49 do i = 1, im 50 50 xd(i, if) = x(i) * fx 51 if(xd(i, if)<xmin) iid(if) = i + 152 if(xd(i, if)<=xmax) ifd(if) = i51 IF(xd(i, if)<xmin) iid(if) = i + 1 52 IF(xd(i, if)<=xmax) ifd(if) = i 53 53 enddo 54 54 PRINT*, 'On stoke du point ', iid(if), ' a ', ifd(if), ' en x' … … 59 59 do j = 1, jm 60 60 yd(j, if) = y(j) * fy 61 if(yd(j, if)>ymax) jid(if) = j + 162 if(yd(j, if)>=ymin) jfd(if) = j61 IF(yd(j, if)>ymax) jid(if) = j + 1 62 IF(yd(j, if)>=ymin) jfd(if) = j 63 63 enddo 64 64 PRINT*, 'On stoke du point ', jid(if), ' a ', jfd(if), ' en y' … … 89 89 90 90 91 end subroutineinigrads91 END SUBROUTINE inigrads -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/initfluxsto.f90
r5114 r5116 47 47 ! Arguments 48 48 ! 49 character(len= *) :: infile50 real:: tstep, t_ops, t_wrt51 integer:: fileid, filevid, filedid49 CHARACTER(LEN = *) :: infile 50 REAL :: tstep, t_ops, t_wrt 51 INTEGER :: fileid, filevid, filedid 52 52 53 53 ! This routine needs IOIPSL to work 54 54 ! Variables locales 55 55 ! 56 real:: nivd(1)57 integer:: tau058 real:: zjulian59 character(len= 3) :: str60 character(len= 10) :: ctrac61 integer:: iq62 real:: rlong(iip1, jjp1), rlat(iip1, jjp1), rl(1, 1)63 integer:: uhoriid, vhoriid, thoriid, zvertiid, dhoriid, dvertiid64 integer:: ii, jj65 integer:: zan, idayref56 REAL :: nivd(1) 57 INTEGER :: tau0 58 REAL :: zjulian 59 CHARACTER(LEN = 3) :: str 60 CHARACTER(LEN = 10) :: ctrac 61 INTEGER :: iq 62 REAL :: rlong(iip1, jjp1), rlat(iip1, jjp1), rl(1, 1) 63 INTEGER :: uhoriid, vhoriid, thoriid, zvertiid, dhoriid, dvertiid 64 INTEGER :: ii, jj 65 INTEGER :: zan, idayref 66 66 logical :: ok_sync 67 67 ! … … 211 211 CALL histend(filevid) 212 212 CALL histend(filedid) 213 if (ok_sync) then213 if (ok_sync) THEN 214 214 CALL histsync(fileid) 215 215 CALL histsync(filevid) … … 217 217 endif 218 218 219 return220 end subroutineinitfluxsto219 RETURN 220 END SUBROUTINE initfluxsto -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inithist.F90
r5114 r5116 46 46 ! Arguments 47 47 ! 48 integer:: day0, anne049 real:: tstep, t_ops, t_wrt48 INTEGER :: day0, anne0 49 REAL :: tstep, t_ops, t_wrt 50 50 51 51 ! This routine needs IOIPSL to work 52 52 ! Variables locales 53 53 ! 54 integer:: tau055 real:: zjulian56 integer:: iq57 real:: rlong(iip1, jjp1), rlat(iip1, jjp1)58 integer:: uhoriid, vhoriid, thoriid, zvertiid59 integer:: ii, jj60 integer:: zan, dayref54 INTEGER :: tau0 55 REAL :: zjulian 56 INTEGER :: iq 57 REAL :: rlong(iip1, jjp1), rlat(iip1, jjp1) 58 INTEGER :: uhoriid, vhoriid, thoriid, zvertiid 59 INTEGER :: ii, jj 60 INTEGER :: zan, dayref 61 61 ! 62 62 ! Initialisations -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inter_barxy_m.F90
r5113 r5116 1 2 1 ! $Id$ 3 2 … … 15 14 SUBROUTINE inter_barxy(dlonid, dlatid, champ, rlonimod, rlatimod, champint) 16 15 17 use lmdz_assert_eq, only: assert_eq18 use lmdz_assert, only: assert16 use lmdz_assert_eq, ONLY: assert_eq 17 use lmdz_assert, ONLY: assert 19 18 20 19 include "dimensions.h" … … 27 26 ! (for "aire", "apoln", "apols") 28 27 29 REAL, intent(in) :: dlonid(:)28 REAL, intent(in) :: dlonid(:) 30 29 ! (longitude from input file, in rad, from -pi to pi) 31 30 32 REAL, intent(in) :: dlatid(:), champ(:, :), rlonimod(:)33 34 REAL, intent(in) :: rlatimod(:)31 REAL, intent(in) :: dlatid(:), champ(:, :), rlonimod(:) 32 33 REAL, intent(in) :: rlatimod(:) 35 34 ! (latitude angle, in degrees or rad, in strictly decreasing order) 36 35 37 real, intent(out) :: champint(:, :)36 real, intent(out) :: champint(:, :) 38 37 ! Si taille de la seconde dim = jjm + 1, on veut interpoler sur les 39 38 ! jjm+1 latitudes rlatu du modele (latitudes des scalaires et de U) … … 55 54 56 55 jnterfd = assert_eq(size(champ, 2) - 1, size(dlatid), & 57 "inter_barxy jnterfd")56 "inter_barxy jnterfd") 58 57 jmods = size(champint, 2) 59 58 CALL assert(size(champ, 1) == size(dlonid), "inter_barxy size(champ, 1)") 60 59 CALL assert((/size(rlonimod), size(champint, 1)/) == iim, & 61 "inter_barxy iim")60 "inter_barxy iim") 62 61 CALL assert(any(jmods == (/jjm, jjm + 1/)), 'inter_barxy jmods') 63 62 CALL assert(size(rlatimod) == jjm, "inter_barxy size(rlatimod)") … … 65 64 ! Check decreasing order for "rlatimod": 66 65 DO i = 2, jjm 67 IF (rlatimod(i) >= rlatimod(i-1)) stop &68 '"inter_barxy": "rlatimod" should be strictly decreasing'66 IF (rlatimod(i) >= rlatimod(i - 1)) stop & 67 '"inter_barxy": "rlatimod" should be strictly decreasing' 69 68 ENDDO 70 69 71 70 yjmod(:jjm) = ord_coordm(rlatimod) 72 71 IF (jmods == jjm + 1) THEN 73 74 '"inter_barxy": with jmods = jjm + 1, yjmod(jjm) should be < 90.'72 IF (90. - yjmod(jjm) < 0.01) stop & 73 '"inter_barxy": with jmods = jjm + 1, yjmod(jjm) should be < 90.' 75 74 ELSE 76 77 78 '"inter_barxy": with jmods = jjm, yjmod(jjm) should be 90.'75 ! jmods = jjm 76 IF (ABS(yjmod(jjm) - 90.) > 0.01) stop & 77 '"inter_barxy": with jmods = jjm, yjmod(jjm) should be 90.' 79 78 ENDIF 80 79 … … 82 81 83 82 DO j = 1, jnterfd + 1 84 85 ENDDO 86 87 CALL ord_coord(dlatid, yjdat, decrois) 83 champy(:, j) = inter_barx(dlonid, champ(:, j), rlonimod) 84 ENDDO 85 86 CALL ord_coord(dlatid, yjdat, decrois) 88 87 IF (decrois) champy(:, :) = champy(:, jnterfd + 1:1:-1) 89 88 DO i = 1, iim 90 89 champint(i, :) = inter_bary(yjdat, champy(i, :), yjmod) 91 90 ENDDO 92 91 champint(:, :) = champint(:, jmods:1:-1) 93 92 94 93 IF (jmods == jjm + 1) THEN 95 96 champint(:, 1) = SUM(aire(:iim,1) * champint(:, 1)) / apoln97 98 * champint(:, jjm + 1)) / apols94 ! Valeurs uniques aux poles 95 champint(:, 1) = SUM(aire(:iim, 1) * champint(:, 1)) / apoln 96 champint(:, jjm + 1) = SUM(aire(:iim, jjm + 1) & 97 * champint(:, jjm + 1)) / apols 99 98 ENDIF 100 99 … … 103 102 !****************************** 104 103 105 function inter_barx(dlonid, fdat, rlonimod) 104 function inter_barx(dlonid, fdat, rlonimod) 106 105 107 106 ! INTERPOLATION BARYCENTRIQUE BASEE SUR LES AIRES … … 117 116 ! ( Les abscisses sont exprimees en degres) 118 117 119 use lmdz_assert_eq, only: assert_eq 118 USE lmdz_assert_eq, ONLY: assert_eq 119 USE lmdz_libmath, ONLY: minmax 120 120 121 121 IMPLICIT NONE 122 122 123 REAL, intent(in) :: dlonid(:)124 real, intent(in) :: fdat(:)125 real, intent(in) :: rlonimod(:)123 REAL, intent(in) :: dlonid(:) 124 real, intent(in) :: fdat(:) 125 real, intent(in) :: rlonimod(:) 126 126 127 127 real inter_barx(size(rlonimod)) … … 130 130 131 131 INTEGER idatmax, imodmax 132 REAL xxid(size(dlonid) +1), xxd(size(dlonid)+1), fdd(size(dlonid)+1)133 REAL fxd(size(dlonid) +1), xchan(size(dlonid)+1), fdchan(size(dlonid)+1)132 REAL xxid(size(dlonid) + 1), xxd(size(dlonid) + 1), fdd(size(dlonid) + 1) 133 REAL fxd(size(dlonid) + 1), xchan(size(dlonid) + 1), fdchan(size(dlonid) + 1) 134 134 REAL xxim(size(rlonimod)) 135 135 … … 149 149 ! A L'INTERFACE OUEST DE LA PREMIERE MAILLE DU MODELE 150 150 DO imod = 1, imodmax 151 152 ENDDO 153 154 CALL minmax( 155 IF( chmax<6.50) THEN156 157 xxim(imod) = xxim(imod) * 180./pi158 151 xxim(imod) = rlonimod(imod) 152 ENDDO 153 154 CALL minmax(imodmax, xxim, chmin, chmax) 155 IF(chmax<6.50) THEN 156 DO imod = 1, imodmax 157 xxim(imod) = xxim(imod) * 180. / pi 158 ENDDO 159 159 ENDIF 160 160 … … 162 162 163 163 DO imod = 1, imodmax 164 165 ENDDO 166 167 idatmax1 = idatmax + 1164 xxim(imod) = xxim(imod) - xim0 165 ENDDO 166 167 idatmax1 = idatmax + 1 168 168 169 169 DO idat = 1, idatmax 170 171 ENDDO 172 173 CALL minmax( 174 IF( chmax<6.50) THEN175 176 xxd(idat) = xxd(idat) * 180./pi177 170 xxd(idat) = dlonid(idat) 171 ENDDO 172 173 CALL minmax(idatmax, xxd, chmin, chmax) 174 IF(chmax<6.50) THEN 175 DO idat = 1, idatmax 176 xxd(idat) = xxd(idat) * 180. / pi 177 ENDDO 178 178 ENDIF 179 179 180 180 DO idat = 1, idatmax 181 xxd(idat) = MOD( xxd(idat) - xim0, 360.)182 181 xxd(idat) = MOD(xxd(idat) - xim0, 360.) 182 fdd(idat) = fdat (idat) 183 183 ENDDO 184 184 185 185 i = 2 186 DO while (xxd(i) >= xxd(i -1) .and. i < idatmax)187 188 ENDDO 189 IF (xxd(i) < xxd(i -1)) THEN190 191 192 nid = idatmax - ichang +1193 194 xchan (i) = xxd(i+ichang -1)195 fdchan(i) = fdd(i+ichang -1)196 197 DO i=1, ichang -1198 xchan (i+ nid) = xxd(i)199 fdchan(i+nid) = fdd(i)200 201 DO i =1, idatmax202 203 204 186 DO while (xxd(i) >= xxd(i - 1) .and. i < idatmax) 187 i = i + 1 188 ENDDO 189 IF (xxd(i) < xxd(i - 1)) THEN 190 ichang = i 191 ! *** reorganisation des longitudes entre 0. et 360. degres **** 192 nid = idatmax - ichang + 1 193 DO i = 1, nid 194 xchan (i) = xxd(i + ichang - 1) 195 fdchan(i) = fdd(i + ichang - 1) 196 ENDDO 197 DO i = 1, ichang - 1 198 xchan (i + nid) = xxd(i) 199 fdchan(i + nid) = fdd(i) 200 ENDDO 201 DO i = 1, idatmax 202 xxd(i) = xchan(i) 203 fdd(i) = fdchan(i) 204 ENDDO 205 205 end IF 206 206 … … 213 213 214 214 DO idat = 1, idatmax 215 IF ( xxd( idatmax1- idat)<360.) exit216 215 IF (xxd(idatmax1 - idat)<360.) exit 216 id1 = id1 + 1 217 217 ENDDO 218 218 219 219 DO idat = 1, idatmax 220 221 220 IF (xxd(idat)>0.) exit 221 id0 = id0 + 1 222 222 END DO 223 223 224 IF( id1 /= 0 ) then225 226 227 fxd (idat) = fdd(idatmax - id1 + idat)228 229 230 231 232 224 IF(id1 /= 0) THEN 225 DO idat = 1, id1 226 xxid(idat) = xxd(idatmax - id1 + idat) - 360. 227 fxd (idat) = fdd(idatmax - id1 + idat) 228 END DO 229 DO idat = 1, idatmax - id1 230 xxid(idat + id1) = xxd(idat) 231 fxd (idat + id1) = fdd(idat) 232 END DO 233 233 end IF 234 234 235 IF(id0 /= 0) then236 237 238 239 240 241 242 xxid (idatmax - id0 + idat) =xxd(idat) + 360.243 fxd (idatmax - id0 + idat) = fdd(idat)244 245 else 246 247 xxid(idat)= xxd(idat)248 fxd (idat)= fdd(idat)249 235 IF(id0 /= 0) THEN 236 DO idat = 1, idatmax - id0 237 xxid(idat) = xxd(idat + id0) 238 fxd (idat) = fdd(idat + id0) 239 END DO 240 241 DO idat = 1, id0 242 xxid (idatmax - id0 + idat) = xxd(idat) + 360. 243 fxd (idatmax - id0 + idat) = fdd(idat) 244 END DO 245 else 246 DO idat = 1, idatmax 247 xxid(idat) = xxd(idat) 248 fxd (idat) = fdd(idat) 249 ENDDO 250 250 end IF 251 251 xxid(idatmax1) = xxid(1) + 360. … … 258 258 ! iteration 259 259 260 x0 261 dxm 260 x0 = xim0 261 dxm = 0. 262 262 imod = 1 263 263 idat = 1 264 264 265 265 do while (imod <= imodmax) 266 267 dx= xxid(idat) - x0268 dxm= dxm + dx269 270 x0= xxid(idat)271 272 273 274 dx= xxim(imod) - x0275 dxm= dxm + dx276 277 x0= xxim(imod)278 dxm= 0.279 280 281 dx= xxim(imod) - x0282 dxm= dxm + dx283 284 x0= xxim(imod)285 dxm= 0.286 287 288 266 do while (xxim(imod)>xxid(idat)) 267 dx = xxid(idat) - x0 268 dxm = dxm + dx 269 inter_barx(imod) = inter_barx(imod) + dx * fxd(idat) 270 x0 = xxid(idat) 271 idat = idat + 1 272 END DO 273 IF (xxim(imod)<xxid(idat)) THEN 274 dx = xxim(imod) - x0 275 dxm = dxm + dx 276 inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm 277 x0 = xxim(imod) 278 dxm = 0. 279 imod = imod + 1 280 ELSE 281 dx = xxim(imod) - x0 282 dxm = dxm + dx 283 inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm 284 x0 = xxim(imod) 285 dxm = 0. 286 imod = imod + 1 287 idat = idat + 1 288 END IF 289 289 END DO 290 290 … … 299 299 ! L'indice 1 correspond à l'interface maille 1 -- maille 2. 300 300 301 use lmdz_assert, only: assert301 use lmdz_assert, ONLY: assert 302 302 303 303 IMPLICIT NONE 304 304 305 REAL, intent(in) :: yjdat(:)305 REAL, intent(in) :: yjdat(:) 306 306 ! (angles, ordonnées des interfaces des mailles des données, in 307 307 ! degrees, in increasing order) 308 308 309 REAL, intent(in) :: fdat(:) ! champ de données310 311 REAL, intent(in) :: yjmod(:)309 REAL, intent(in) :: fdat(:) ! champ de données 310 311 REAL, intent(in) :: yjmod(:) 312 312 ! (ordonnées des interfaces des mailles du modèle) 313 313 ! (in degrees, in strictly increasing order) … … 317 317 ! Variables local to the procedure: 318 318 319 REAL y0, dy, dym 319 REAL y0, dy, dym 320 320 INTEGER jdat ! indice du champ de données 321 321 integer jmod ! indice du champ du modèle … … 327 327 ! Initialisation des variables 328 328 inter_bary(:) = 0. 329 y0 330 dym 331 jmod 332 jdat 329 y0 = -90. 330 dym = 0. 331 jmod = 1 332 jdat = 1 333 333 334 334 do while (jmod <= size(yjmod)) 335 336 dy= yjdat(jdat) - y0337 dym= dym + dy338 339 y0= yjdat(jdat)340 jdat= jdat + 1341 342 343 dy= yjmod(jmod) - y0344 dym= dym + dy345 346 y0= yjmod(jmod)347 dym= 0.348 jmod= jmod + 1349 350 351 dy= yjmod(jmod) - y0352 dym= dym + dy353 354 y0= yjmod(jmod)355 dym= 0.356 jmod= jmod + 1357 jdat= jdat + 1358 335 do while (yjmod(jmod) > yjdat(jdat)) 336 dy = yjdat(jdat) - y0 337 dym = dym + dy 338 inter_bary(jmod) = inter_bary(jmod) + dy * fdat(jdat) 339 y0 = yjdat(jdat) 340 jdat = jdat + 1 341 END DO 342 IF (yjmod(jmod) < yjdat(jdat)) THEN 343 dy = yjmod(jmod) - y0 344 dym = dym + dy 345 inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym 346 y0 = yjmod(jmod) 347 dym = 0. 348 jmod = jmod + 1 349 ELSE 350 ! {yjmod(jmod) == yjdat(jdat)} 351 dy = yjmod(jmod) - y0 352 dym = dym + dy 353 inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym 354 y0 = yjmod(jmod) 355 dym = 0. 356 jmod = jmod + 1 357 jdat = jdat + 1 358 END IF 359 359 END DO 360 360 ! Le test de fin suppose que l'interface 0 est commune aux deux … … 373 373 ! Finally, the procedure adds 90° as the last value of the array. 374 374 375 use lmdz_assert_eq, only: assert_eq376 use comconst_mod, only: pi375 use lmdz_assert_eq, ONLY: assert_eq 376 use comconst_mod, ONLY: pi 377 377 378 378 IMPLICIT NONE 379 379 380 REAL, intent(in) :: xi(:)380 REAL, intent(in) :: xi(:) 381 381 ! (latitude, in degrees or radians, in increasing or decreasing order) 382 382 ! ("xi" should contain latitudes from pole to pole. … … 385 385 ! So the extreme values should not be 90° and -90°.) 386 386 387 REAL, intent(out) :: xo(:) ! angles in degrees388 LOGICAL, intent(out) :: decrois387 REAL, intent(out) :: xo(:) ! angles in degrees 388 LOGICAL, intent(out) :: decrois 389 389 390 390 ! Variables local to the procedure: … … 398 398 decrois = xi(2) < xi(1) 399 399 DO i = 3, nmax 400 IF (decrois .neqv. xi(i) < xi(i-1)) stop &401 '"ord_coord": latitudes are not monotonic'402 ENDDO 403 404 IF (abs(xi(1)) < pi) then405 406 400 IF (decrois .neqv. xi(i) < xi(i - 1)) stop & 401 '"ord_coord": latitudes are not monotonic' 402 ENDDO 403 404 IF (abs(xi(1)) < pi) THEN 405 ! "xi" contains latitudes in radians 406 xo(:nmax) = xi(:) * 180. / pi 407 407 else 408 409 408 ! "xi" contains latitudes in degrees 409 xo(:nmax) = xi(:) 410 410 end IF 411 411 412 412 IF (ABS(abs(xo(1)) - 90) < 0.001 .or. ABS(abs(xo(nmax)) - 90) < 0.001) THEN 413 414 415 // 'grid cells, not the centers of grid cells.'416 413 print *, "ord_coord" 414 PRINT *, '"xi" should contain the latitudes of the boundaries of ' & 415 // 'grid cells, not the centers of grid cells.' 416 STOP 417 417 ENDIF 418 418 … … 429 429 ! order. 430 430 431 use comconst_mod, only: pi431 use comconst_mod, ONLY: pi 432 432 433 433 IMPLICIT NONE 434 434 435 REAL, intent(in) :: xi(:) ! angle, in rad or degrees435 REAL, intent(in) :: xi(:) ! angle, in rad or degrees 436 436 REAL ord_coordm(size(xi)) ! angle, in degrees 437 437 … … 439 439 440 440 IF (xi(1) < 6.5) THEN 441 442 441 ! "xi" is in rad 442 ord_coordm(:) = xi(size(xi):1:-1) * 180. / pi 443 443 else 444 445 444 ! "xi" is in degrees 445 ord_coordm(:) = xi(size(xi):1:-1) 446 446 ENDIF 447 447 -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/interpost.f90
r5113 r5116 12 12 13 13 ! Arguments 14 real:: q(iip1,jjp1,llm)15 real:: qppm(iim,jjp1,llm)14 REAL :: q(iip1,jjp1,llm) 15 REAL :: qppm(iim,jjp1,llm) 16 16 ! Local 17 integer:: l,i,j17 INTEGER :: l,i,j 18 18 19 19 ! RE-INVERSION DES NIVEAUX … … 41 41 return 42 42 43 end subroutineinterpost43 END SUBROUTINE interpost -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/interpre.f90
r5114 r5116 18 18 !--------------------------------------------------- 19 19 ! Arguments 20 real:: apppm(llm + 1), bpppm(llm + 1)21 real:: q(iip1, jjp1, llm), qppm(iim, jjp1, llm)20 REAL :: apppm(llm + 1), bpppm(llm + 1) 21 REAL :: q(iip1, jjp1, llm), qppm(iim, jjp1, llm) 22 22 !--------------------------------------------------- 23 real:: masse(iip1, jjp1, llm)24 real:: massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm)25 real:: w(iip1, jjp1, llm)26 real:: fluxwppm(iim, jjp1, llm)27 real:: pbaru(iip1, jjp1, llm)28 real:: pbarv(iip1, jjm, llm)29 real:: unatppm(iim, jjp1, llm)30 real:: vnatppm(iim, jjp1, llm)31 real:: psppm(iim, jjp1)23 REAL :: masse(iip1, jjp1, llm) 24 REAL :: massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm) 25 REAL :: w(iip1, jjp1, llm) 26 REAL :: fluxwppm(iim, jjp1, llm) 27 REAL :: pbaru(iip1, jjp1, llm) 28 REAL :: pbarv(iip1, jjm, llm) 29 REAL :: unatppm(iim, jjp1, llm) 30 REAL :: vnatppm(iim, jjp1, llm) 31 REAL :: psppm(iim, jjp1) 32 32 !--------------------------------------------------- 33 33 ! Local 34 real:: vnat(iip1, jjp1, llm)35 real:: unat(iip1, jjp1, llm)36 real:: fluxw(iip1, jjp1, llm)37 real:: smass(iip1, jjp1)34 REAL :: vnat(iip1, jjp1, llm) 35 REAL :: unat(iip1, jjp1, llm) 36 REAL :: fluxw(iip1, jjp1, llm) 37 REAL :: smass(iip1, jjp1) 38 38 !---------------------------------------------------- 39 integer:: l, ij, i, j39 INTEGER :: l, ij, i, j 40 40 41 41 ! CALCUL DE LA PRESSION DE SURFACE … … 118 118 enddo 119 119 120 return121 end subroutineinterpre120 RETURN 121 END SUBROUTINE interpre 122 122 123 123 -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/invert_zoom_x_m.F90
r5114 r5116 9 9 SUBROUTINE invert_zoom_x(xf, xtild, Xprimt, xlon, xprimm, xuv) 10 10 11 use lmdz_coefpoly, only: coefpoly12 use nrtype, only: pi, pi_d, twopi_d, k813 use serre_mod, only: clon11 use lmdz_coefpoly, ONLY: coefpoly 12 use nrtype, ONLY: pi, pi_d, twopi_d, k8 13 use serre_mod, ONLY: clon 14 14 15 15 include "dimensions.h" … … 64 64 end DO 65 65 66 if (ABS(xvrai(i) - xo1) > my_eps) then66 if (ABS(xvrai(i) - xo1) > my_eps) THEN 67 67 ! iter == 300 68 68 print *, 'Pas de solution.' -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/iso_verif_dyn.f90
r5113 r5116 5 5 6 6 ! input: 7 real:: x8 character(len=*) :: err_msg ! message d''erreur à afficher7 REAL :: x 8 CHARACTER(LEN=*) :: err_msg ! message d''erreur à afficher 9 9 10 10 ! output 11 real:: borne11 REAL :: borne 12 12 parameter (borne=1e19) 13 integer:: iso_verif_noNaN_nostop13 INTEGER :: iso_verif_noNaN_nostop 14 14 15 if ((x>-borne).and.(x<borne)) then15 if ((x>-borne).and.(x<borne)) THEN 16 16 iso_verif_noNAN_nostop=0 17 17 else 18 write(*,*) 'erreur detectee par iso_verif_nonNaN:'19 write(*,*) err_msg20 write(*,*) 'x=',x18 WRITE(*,*) 'erreur detectee par iso_verif_nonNaN:' 19 WRITE(*,*) err_msg 20 WRITE(*,*) 'x=',x 21 21 iso_verif_noNaN_nostop=1 22 22 endif 23 23 24 return25 end functioniso_verif_nonan_nostop24 RETURN 25 END FUNCTION iso_verif_nonan_nostop 26 26 27 27 function iso_verif_egalite_nostop & … … 33 33 34 34 ! input: 35 real:: a, b36 character(len=*) :: err_msg ! message d''erreur à afficher35 REAL :: a, b 36 CHARACTER(LEN=*) :: err_msg ! message d''erreur à afficher 37 37 38 38 ! locals 39 real:: errmax ! erreur maximale en absolu.40 real:: errmaxrel ! erreur maximale en relatif autorisée39 REAL :: errmax ! erreur maximale en absolu. 40 REAL :: errmaxrel ! erreur maximale en relatif autorisée 41 41 parameter (errmax=1e-8) 42 42 parameter (errmaxrel=1e-3) 43 43 44 44 ! output 45 integer:: iso_verif_egalite_nostop45 INTEGER :: iso_verif_egalite_nostop 46 46 47 47 iso_verif_egalite_nostop=0 48 48 49 if (abs(a-b)>errmax) then49 if (abs(a-b)>errmax) THEN 50 50 if (abs((a-b)/max(max(abs(b),abs(a)),1e-18)) & 51 >errmaxrel) then52 write(*,*) 'erreur detectee par iso_verif_egalite:'53 write(*,*) err_msg54 write(*,*) 'a=',a55 write(*,*) 'b=',b51 >errmaxrel) THEN 52 WRITE(*,*) 'erreur detectee par iso_verif_egalite:' 53 WRITE(*,*) err_msg 54 WRITE(*,*) 'a=',a 55 WRITE(*,*) 'b=',b 56 56 iso_verif_egalite_nostop=1 57 57 endif 58 58 endif 59 59 60 return61 end functioniso_verif_egalite_nostop60 RETURN 61 END FUNCTION iso_verif_egalite_nostop 62 62 63 63 … … 68 68 69 69 ! input: 70 real:: x,q71 integer:: iso ! 2=HDO, 1=O1872 character(len=*) :: err_msg ! message d''erreur à afficher70 REAL :: x,q 71 INTEGER :: iso ! 2=HDO, 1=O18 72 CHARACTER(LEN=*) :: err_msg ! message d''erreur à afficher 73 73 74 74 ! locals 75 real:: qmin,deltaD76 real:: deltaDmax,deltaDmin,tnat75 REAL :: qmin,deltaD 76 REAL :: deltaDmax,deltaDmin,tnat 77 77 parameter (qmin=1e-11) 78 78 parameter (deltaDmax=200.0,deltaDmin=-999.9) 79 79 80 80 ! output 81 integer:: iso_verif_aberrant_nostop81 INTEGER :: iso_verif_aberrant_nostop 82 82 83 83 iso_verif_aberrant_nostop=0 84 84 85 85 ! verifier que HDO est raisonable 86 if (q>qmin) then86 if (q>qmin) THEN 87 87 IF(getKey('tnat', tnat, isoName(iso))) THEN 88 88 err_msg = 'Missing isotopic parameter "tnat"' … … 91 91 END IF 92 92 deltaD=(x/q/tnat-1)*1000 93 if ((deltaD>deltaDmax).or.(deltaD<deltaDmin)) then94 write(*,*) 'erreur detectee par iso_verif_aberrant:'95 write(*,*) err_msg96 write(*,*) 'q=',q97 write(*,*) 'deltaD=',deltaD98 write(*,*) 'iso=',iso93 if ((deltaD>deltaDmax).or.(deltaD<deltaDmin)) THEN 94 WRITE(*,*) 'erreur detectee par iso_verif_aberrant:' 95 WRITE(*,*) err_msg 96 WRITE(*,*) 'q=',q 97 WRITE(*,*) 'deltaD=',deltaD 98 WRITE(*,*) 'iso=',iso 99 99 iso_verif_aberrant_nostop=1 100 endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then 101 endif !if (q(i,k,iq).gt.qmin) then 100 endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) THEN 101 endif !if (q(i,k,iq).gt.qmin) THEN 102 RETURN 103 END FUNCTION iso_verif_aberrant_nostop 102 104 103 104 return105 end function iso_verif_aberrant_nostop106 -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/limx.f90
r5105 r5116 22 22 ! Arguments: 23 23 ! ---------- 24 real:: pente_max24 REAL :: pente_max 25 25 REAL :: s0(ip1jmp1,llm),sm(ip1jmp1,llm) 26 real:: sx(ip1jmp1,llm)26 REAL :: sx(ip1jmp1,llm) 27 27 ! 28 28 ! Local … … 30 30 ! 31 31 INTEGER :: ij,l,j,i,iju,ijq,indu(ip1jmp1),niju 32 integer:: n0,iadvplus(ip1jmp1,llm),nl(llm)32 INTEGER :: n0,iadvplus(ip1jmp1,llm),nl(llm) 33 33 ! 34 34 REAL :: q(ip1jmp1,llm) 35 real:: dxq(ip1jmp1,llm)35 REAL :: dxq(ip1jmp1,llm) 36 36 37 37 38 38 REAL :: new_m,zm 39 real:: dxqu(ip1jmp1)40 real:: adxqu(ip1jmp1),dxqmax(ip1jmp1)39 REAL :: dxqu(ip1jmp1) 40 REAL :: adxqu(ip1jmp1),dxqmax(ip1jmp1) 41 41 42 42 Logical :: extremum,first … … 44 44 45 45 REAL :: SSUM 46 integer:: ismax,ismin46 INTEGER :: ismax,ismin 47 47 EXTERNAL SSUM, ismin,ismax 48 48 … … 84 84 85 85 do ij=iip2+1,ip1jm 86 if( dxqu(ij-1)*dxqu(ij)>0. &87 .and. dxq(ij,l)*dxqu(ij)>0.) then86 IF( dxqu(ij-1)*dxqu(ij)>0. & 87 .and. dxq(ij,l)*dxqu(ij)>0.) THEN 88 88 dxq(ij,l)= & 89 89 sign(min(abs(dxq(ij,l)),dxqmax(ij)),dxq(ij,l)) -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/limy.f90
r5105 r5116 15 15 ! -------------------------------------------------------------------- 16 16 USE comconst_mod, ONLY: pi 17 USE lmdz_libmath, ONLY: ismax, ismin 17 18 IMPLICIT NONE 18 19 ! … … 24 25 ! Arguments: 25 26 ! ---------- 26 real:: pente_max27 real:: s0(ip1jmp1,llm),sy(ip1jmp1,llm),sm(ip1jmp1,llm)27 REAL :: pente_max 28 REAL :: s0(ip1jmp1,llm),sy(ip1jmp1,llm),sm(ip1jmp1,llm) 28 29 ! 29 30 ! Local … … 34 35 REAL :: q(ip1jmp1,llm) 35 36 REAL :: airej2,airejjm,airescb(iim),airesch(iim) 36 real:: sigv,dyq(ip1jmp1),dyqv(ip1jm)37 real:: adyqv(ip1jm),dyqmax(ip1jmp1)37 REAL :: sigv,dyq(ip1jmp1),dyqv(ip1jm) 38 REAL :: adyqv(ip1jm),dyqmax(ip1jmp1) 38 39 REAL :: qbyv(ip1jm,llm) 39 40 … … 42 43 save first 43 44 44 real:: convpn,convps,convmpn,convmps45 real:: sinlon(iip1),sinlondlon(iip1)46 real:: coslon(iip1),coslondlon(iip1)45 REAL :: convpn,convps,convmpn,convmps 46 REAL :: sinlon(iip1),sinlondlon(iip1) 47 REAL :: coslon(iip1),coslondlon(iip1) 47 48 save sinlon,coslon,sinlondlon,coslondlon 48 49 ! 49 50 ! 50 51 REAL :: SSUM 51 integer :: ismax,ismin 52 EXTERNAL SSUM, convflu,ismin,ismax 53 EXTERNAL filtreg 52 EXTERNAL SSUM 54 53 55 54 data first/.TRUE./ 56 55 57 if(first) then56 IF(first) THEN 58 57 PRINT*,'SCHEMA AMONT NOUVEAU' 59 58 first=.FALSE. … … 129 128 ! cas ou on a un extremum au pole 130 129 131 ! if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)130 ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.) 132 131 ! & appn=0. 133 ! if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*132 ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* 134 133 ! & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.) 135 134 ! & apps=0. … … 150 149 ! enddo 151 150 152 if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1))<=0.) &153 then151 IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1))<=0.) & 152 THEN 154 153 do ij=1,iip1 155 154 dyqmax(ij)=0. … … 161 160 endif 162 161 163 if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* &162 IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* & 164 163 dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)<=0.) & 165 then164 THEN 166 165 do ij=ip1jm+1,ip1jmp1 167 166 dyqmax(ij)=0. … … 176 175 177 176 do ij=1,ip1jmp1 ! cf below: should it be ip1jm instead ? 178 if(dyqv(ij)*dyqv(ij-iip1)>0.) then ! /!\ causes Warning: iteration 1056 invokes undefined behavior [-Waggressive-loop-optimizations] in 32x32x39177 IF(dyqv(ij)*dyqv(ij-iip1)>0.) then ! /!\ causes Warning: iteration 1056 invokes undefined behavior [-Waggressive-loop-optimizations] in 32x32x39 179 178 dyq(ij)=sign(min(abs(dyq(ij)),dyqmax(ij)),dyq(ij)) 180 179 else -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/limz.f90
r5105 r5116 22 22 ! Arguments: 23 23 ! ---------- 24 real:: pente_max24 REAL :: pente_max 25 25 REAL :: s0(ip1jmp1,llm),sm(ip1jmp1,llm) 26 real:: sz(ip1jmp1,llm)26 REAL :: sz(ip1jmp1,llm) 27 27 ! 28 28 ! Local … … 30 30 ! 31 31 INTEGER :: ij,l,j,i,iju,ijq,indu(ip1jmp1),niju 32 integer:: n0,iadvplus(ip1jmp1,llm),nl(llm)32 INTEGER :: n0,iadvplus(ip1jmp1,llm),nl(llm) 33 33 ! 34 34 REAL :: q(ip1jmp1,llm) 35 real:: dzq(ip1jmp1,llm)35 REAL :: dzq(ip1jmp1,llm) 36 36 37 37 38 38 REAL :: new_m,zm 39 real:: dzqw(ip1jmp1)40 real:: adzqw(ip1jmp1),dzqmax(ip1jmp1)39 REAL :: dzqw(ip1jmp1) 40 REAL :: adzqw(ip1jmp1),dzqmax(ip1jmp1) 41 41 42 42 Logical :: extremum,first … … 44 44 45 45 REAL :: SSUM 46 integer:: ismax,ismin46 INTEGER :: ismax,ismin 47 47 EXTERNAL SSUM, ismin,ismax 48 48 … … 77 77 78 78 do l=2,llm-1 79 if( dzqw(l-1)*dzqw(l)>0. &80 .and. dzq(ij,l)*dzqw(l)>0.) then79 IF( dzqw(l-1)*dzqw(l)>0. & 80 .and. dzq(ij,l)*dzqw(l)>0.) THEN 81 81 dzq(ij,l)= & 82 82 sign(min(abs(dzq(ij,l)),dzqmax(l)),dzq(ij,l)) -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/pentes_ini.f90
r5106 r5116 31 31 ! Arguments: 32 32 ! ---------- 33 integer:: mode33 INTEGER :: mode 34 34 REAL :: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm ) 35 35 REAL :: q( iip1,jjp1,llm,0:3) … … 42 42 REAL :: s0( iip1,jjp1,llm ), sx( iip1,jjp1,llm ) 43 43 REAL :: sy( iip1,jjp1,llm ), sz( iip1,jjp1,llm ) 44 real:: masn,mass,zz44 REAL :: masn,mass,zz 45 45 INTEGER :: i,j,l,iq 46 46 47 47 ! modif Fred 24 03 96 48 48 49 real:: sinlon(iip1),sinlondlon(iip1)50 real:: coslon(iip1),coslondlon(iip1)49 REAL :: sinlon(iip1),sinlondlon(iip1) 50 REAL :: coslon(iip1),coslondlon(iip1) 51 51 save sinlon,coslon,sinlondlon,coslondlon 52 real:: dyn1,dyn2,dys1,dys253 real:: qpn,qps,dqzpn,dqzps54 real:: smn,sms,s0n,s0s,sxn(iip1),sxs(iip1)55 real:: qmin,zq,pente_max52 REAL :: dyn1,dyn2,dys1,dys2 53 REAL :: qpn,qps,dqzpn,dqzps 54 REAL :: smn,sms,s0n,s0s,sxn(iip1),sxs(iip1) 55 REAL :: qmin,zq,pente_max 56 56 ! 57 57 REAL :: SSUM 58 integer:: ismax,ismin,lati,latf58 INTEGER :: ismax,ismin,lati,latf 59 59 EXTERNAL SSUM, ismin,ismax 60 60 logical :: first … … 72 72 limit = .TRUE. 73 73 pente_max=2 74 ! if (mode.eq.1.or.mode.eq.3) then75 ! if (mode.eq.1) then76 if (mode>=1) then74 ! if (mode.eq.1.or.mode.eq.3) THEN 75 ! if (mode.eq.1) THEN 76 if (mode>=1) THEN 77 77 lati=2 78 78 latf=jjm … … 84 84 qmin=0.4995 85 85 qmin=0. 86 if(first) then86 IF(first) THEN 87 87 PRINT*,'SCHEMA AMONT NOUVEAU' 88 88 first=.FALSE. … … 181 181 ! do i=1,iip1 182 182 ! zq=s0(i,j,l)/sm(i,j,l) 183 ! if(zq.lt.qmin)183 ! IF(zq.lt.qmin) 184 184 ! , PRINT*,'avant advx1, s0(',i,',',j,',',l,')=',zq 185 185 ! enddo … … 187 187 ! enddo 188 188 !CC 189 if(mode==2) then189 IF(mode==2) THEN 190 190 do l=1,llm 191 191 s0s=0. … … 247 247 endif 248 248 249 if (mode==4) then249 if (mode==4) THEN 250 250 do l=1,llm 251 251 do i=1,iip1 … … 261 261 CALL advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf) 262 262 ! CALL minmaxq(zq,1.e33,-1.e33,'avant advy ') 263 if (mode==4) then263 if (mode==4) THEN 264 264 do l=1,llm 265 265 do i=1,iip1 … … 282 282 CALL limz(s0,sz,sm,pente_max) 283 283 CALL advz( limit,dtvr,w,sm,s0,sx,sy,sz ) 284 if (mode==4) then284 if (mode==4) THEN 285 285 do l=1,llm 286 286 do i=1,iip1 … … 306 306 307 307 ! CALL minmaxq(zq,1.e33,-1.e33,'avant advx ') 308 if (mode==4) then308 if (mode==4) THEN 309 309 do l=1,llm 310 310 do i=1,iip1 … … 323 323 ! do i=1,iip1 324 324 ! zq=s0(i,j,l)/sm(i,j,l) 325 ! if(zq.lt.qmin)325 ! IF(zq.lt.qmin) 326 326 ! , PRINT*,'apres advx2, s0(',i,',',j,',',l,')=',zq 327 327 ! enddo … … 346 346 ! Traitements specifiques au pole 347 347 348 if(mode>=1) then348 IF(mode>=1) THEN 349 349 DO l=1,llm 350 350 ! filtrages aux poles … … 361 361 q( i,jjp1,llm+1-l,0)=qps 362 362 enddo 363 if(mode==3) then363 IF(mode==3) THEN 364 364 dyn1=0. 365 365 dys1=0. … … 382 382 enddo 383 383 endif 384 if(mode==1) then384 IF(mode==1) THEN 385 385 ! on filtre les valeurs au bord de la "grande maille pole" 386 386 dyn1=0. … … 459 459 do j=1,jjp1 460 460 do i=1,iip1 461 if(q(i,j,l,0)<qmin) &461 IF(q(i,j,l,0)<qmin) & 462 462 PRINT*,'apres pentes, s0(',i,',',j,',',l,')=',q(i,j,l,0) 463 463 enddo -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/ppm3d.f90
r5113 r5116 276 276 ! 277 277 278 real:: Q(IMR,JNP,NLAY,NC),PS1(IMR,JNP),PS2(IMR,JNP), &278 REAL :: Q(IMR,JNP,NLAY,NC),PS1(IMR,JNP),PS2(IMR,JNP), & 279 279 U(IMR,JNP,NLAY),V(IMR,JNP,NLAY),AP(NLAY+1), & 280 280 BP(NLAY+1),W(IMR,JNP,NLAY),NDT,val(NLAY),Umax 281 integer:: IGD,IORD,JORD,KORD,NC,IMR,JNP,j1,NLAY,AE282 integer:: IMRD2283 real:: PT281 INTEGER :: IGD,IORD,JORD,KORD,NC,IMR,JNP,j1,NLAY,AE 282 INTEGER :: IMRD2 283 REAL :: PT 284 284 logical :: cross, fill, dum 285 285 ! 286 286 ! Local dynamic arrays 287 287 ! 288 real:: CRX(IMR,JNP),CRY(IMR,JNP),xmass(IMR,JNP),ymass(IMR,JNP), &288 REAL :: CRX(IMR,JNP),CRY(IMR,JNP),xmass(IMR,JNP),ymass(IMR,JNP), & 289 289 fx1(IMR+1),DPI(IMR,JNP,NLAY),delp1(IMR,JNP,NLAY), & 290 290 WK1(IMR,JNP,NLAY),PU(IMR,JNP),PV(IMR,JNP),DC2(IMR,JNP), & … … 294 294 ! Local static arrays 295 295 ! 296 real:: DTDX(Jmax), DTDX5(Jmax), acosp(Jmax), &296 REAL :: DTDX(Jmax), DTDX5(Jmax), acosp(Jmax), & 297 297 cosp(Jmax), cose(Jmax), DAP(kmax),DBK(Kmax) 298 298 data NDT0, NSTEP /0, 0/ … … 314 314 ! 315 315 ! *********** Initialization ********************** 316 if(NSTEP==1) then317 ! 318 write(6,*) '------------------------------------ '319 write(6,*) 'NASA/GSFC Transport Core Version 4.5'320 write(6,*) '------------------------------------ '316 IF(NSTEP==1) THEN 317 ! 318 WRITE(6,*) '------------------------------------ ' 319 WRITE(6,*) 'NASA/GSFC Transport Core Version 4.5' 320 WRITE(6,*) '------------------------------------ ' 321 321 ! 322 322 WRITE(6,*) 'IMR=',IMR,' JNP=',JNP,' NLAY=',NLAY,' j1=',j1 … … 324 324 ! 325 325 ! controles sur les parametres 326 if(NLAY<6) then327 write(6,*) 'NLAY must be >= 6'326 IF(NLAY<6) THEN 327 WRITE(6,*) 'NLAY must be >= 6' 328 328 stop 329 329 endif 330 if (JNP<NLAY) then331 write(6,*) 'JNP must be >= NLAY'330 if (JNP<NLAY) THEN 331 WRITE(6,*) 'JNP must be >= NLAY' 332 332 stop 333 333 endif 334 334 IMRD2=mod(IMR,2) 335 if (j1==2.and.IMRD2/=0) then336 write(6,*) 'if j1=2 IMR must be an even integer'335 if (j1==2.and.IMRD2/=0) THEN 336 WRITE(6,*) 'if j1=2 IMR must be an even integer' 337 337 stop 338 338 endif 339 339 340 340 ! 341 if(Jmax<JNP .or. Kmax<NLAY) then342 write(6,*) 'Jmax or Kmax is too small'341 IF(Jmax<JNP .or. Kmax<NLAY) THEN 342 WRITE(6,*) 'Jmax or Kmax is too small' 343 343 stop 344 344 endif … … 353 353 DP = PI / REAL(JMR) 354 354 ! 355 if(IGD==0) then355 IF(IGD==0) THEN 356 356 ! Compute analytic cosine at cell edges 357 357 CALL cosa(cosp,cose,JNP,PI,DP) … … 372 372 endif 373 373 ! 374 if(NDT0 /= NDT) then374 IF(NDT0 /= NDT) THEN 375 375 DT = NDT 376 376 NDT0 = NDT 377 377 378 if(Umax < 180.) then379 write(6,*) 'Umax may be too small!'378 IF(Umax < 180.) THEN 379 WRITE(6,*) 'Umax may be too small!' 380 380 endif 381 381 CR1 = abs(Umax*DT)/(DL*AE) 382 382 MaxDT = DP*AE / abs(Umax) + 0.5 383 write(6,*)'Largest time step for max(V)=',Umax,' is ',MaxDT384 if(MaxDT < abs(NDT)) then385 write(6,*) 'Warning!!! NDT maybe too large!'386 endif 387 ! 388 if(CR1>=0.95) then383 WRITE(6,*)'Largest time step for max(V)=',Umax,' is ',MaxDT 384 IF(MaxDT < abs(NDT)) THEN 385 WRITE(6,*) 'Warning!!! NDT maybe too large!' 386 endif 387 ! 388 IF(CR1>=0.95) THEN 389 389 JS0 = 0 390 390 JN0 = 0 … … 413 413 DTDY5 = 0.5*DTDY 414 414 ! 415 ! write(6,*) 'J1=',J1,' J2=', J2415 ! WRITE(6,*) 'J1=',J1,' J2=', J2 416 416 endif 417 417 ! … … 429 429 430 430 ! 431 if(j1/=2) then431 IF(j1/=2) THEN 432 432 DO IC=1,NC 433 433 DO L=1,NLAY … … 453 453 do k=1,NLAY 454 454 ! 455 if(IGD==0) then455 IF(IGD==0) THEN 456 456 ! Convert winds on A-Grid to Courant # on C-Grid. 457 457 CALL A2C(U(1,1,k),V(1,1,k),IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5) … … 480 480 do j=JS0,j1+1,-1 481 481 do i=1,IMR 482 if(abs(CRX(i,j))>1.) then482 IF(abs(CRX(i,j))>1.) THEN 483 483 JS = j 484 484 go to 2222 … … 490 490 do j=JN0,j2-1 491 491 do i=1,IMR 492 if(abs(CRX(i,j))>1.) then492 IF(abs(CRX(i,j))>1.) THEN 493 493 JN = j 494 494 go to 2233 … … 498 498 2233 continue 499 499 ! 500 if(j1/=2) then ! Enlarged polar cap.500 IF(j1/=2) then ! Enlarged polar cap. 501 501 do i=1,IMR 502 502 DPI(i, 2,k) = 0. … … 586 586 enddo 587 587 ! 588 if(j1==2) then588 IF(j1==2) THEN 589 589 IMH = IMR/2 590 590 do i=1,IMH … … 608 608 ! E-W advective cross term 609 609 do j=J1,J2 610 if(J>JS .and. J<JN) GO TO 250610 IF(J>JS .and. J<JN) GO TO 250 611 611 ! 612 612 do i=1,IMR … … 623 623 ru = UA(i,j) - iu 624 624 iiu = i-iu 625 if(UA(i,j)>=0.) then625 IF(UA(i,j)>=0.) THEN 626 626 wk1(i,j,1) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu)) 627 627 else … … 633 633 END DO 634 634 ! 635 if(JN/=0) then635 IF(JN/=0) THEN 636 636 do j=JS+1,JN-1 637 637 ! … … 661 661 enddo 662 662 ! 663 if(cross) then663 IF(cross) THEN 664 664 ! Add cross terms in the vertical direction. 665 if(IORD >= 2) then665 IF(IORD >= 2) THEN 666 666 iad = 2 667 667 else … … 669 669 endif 670 670 ! 671 if(JORD >= 2) then671 IF(JORD >= 2) THEN 672 672 jad = 2 673 673 else … … 750 750 ! 751 751 752 if(fill) CALL qckxyz(DQ(1,1,1,IC),DC2,IMR,JNP,NLAY,j1,j2, &752 IF(fill) CALL qckxyz(DQ(1,1,1,IC),DC2,IMR,JNP,NLAY,j1,j2, & 753 753 cosp,acosp,.FALSE.,IC,NSTEP) 754 754 ! … … 765 765 enddo 766 766 ! 767 if(j1/=2) then767 IF(j1/=2) THEN 768 768 DO k=1,NLAY 769 769 DO I=1,IMR … … 776 776 END DO 777 777 ! 778 if(j1/=2) then778 IF(j1/=2) THEN 779 779 DO k=1,NLAY 780 780 DO i=1,IMR … … 794 794 integer,parameter :: kmax = 150 795 795 real,parameter :: R23 = 2./3., R3 = 1./3. 796 integer:: IMR,JNP,NLAY,J1,KORD797 real:: WZ(IMR,JNP,NLAY),P(IMR,JNP,NLAY),DC(IMR,JNP,NLAY), &796 INTEGER :: IMR,JNP,NLAY,J1,KORD 797 REAL :: WZ(IMR,JNP,NLAY),P(IMR,JNP,NLAY),DC(IMR,JNP,NLAY), & 798 798 wk1(IMR,*),delp(IMR,JNP,NLAY),DQ(IMR,JNP,NLAY), & 799 799 DQDT(IMR,JNP,NLAY) 800 800 ! Assuming JNP >= NLAY 801 real:: AR(IMR,*),AL(IMR,*),A6(IMR,*),flux(IMR,*),wk2(IMR,*), &801 REAL :: AR(IMR,*),AL(IMR,*),A6(IMR,*),flux(IMR,*),wk2(IMR,*), & 802 802 wz2(IMR,*) 803 integer:: JMR,IMJM,NLAYM1,LMT,K,I,J804 real:: c0,c1,c2,tmp,qmax,qmin,a,b,fct,a1,a2,cm,cp803 INTEGER :: JMR,IMJM,NLAYM1,LMT,K,I,J 804 REAL :: c0,c1,c2,tmp,qmax,qmin,a,b,fct,a1,a2,cm,cp 805 805 ! 806 806 JMR = JNP - 1 … … 871 871 ! 872 872 ! Check if change sign 873 if(wk1(i,1)*AL(i,1)<=0.) then873 IF(wk1(i,1)*AL(i,1)<=0.) THEN 874 874 AL(i,1) = 0. 875 875 flux(i,1) = 0. … … 887 887 AR(i,NLAY) = wk1(i,NLAY) + fct 888 888 AL(i,NLAY) = wk1(i,NLAY) - (fct+fct) 889 if(wk1(i,NLAY)*AR(i,NLAY)<=0.) AR(i,NLAY) = 0.889 IF(wk1(i,NLAY)*AR(i,NLAY)<=0.) AR(i,NLAY) = 0. 890 890 flux(i,NLAY) = AR(i,NLAY) - wk1(i,NLAY) 891 891 END DO … … 926 926 ! 927 927 ! Interior depending on KORD 928 if(LMT<=2) &928 IF(LMT<=2) & 929 929 CALL lmtppm(flux(1,2),A6(1,2),AR(1,2),AL(1,2),wk1(1,2), & 930 930 IMR*(NLAY-2),LMT) … … 933 933 ! 934 934 DO i=1,IMR*NLAYM1 935 IF(wz2(i,1)>0.) then935 IF(wz2(i,1)>0.) THEN 936 936 CM = wz2(i,1) / wk2(i,1) 937 937 flux(i,2) = AR(i,1)+0.5*CM*(AL(i,1)-AR(i,1)+A6(i,1)*(1.-R23*CM)) … … 961 961 2000 continue 962 962 END DO 963 return964 end subroutinefzppm963 RETURN 964 END SUBROUTINE fzppm 965 965 ! 966 966 SUBROUTINE xtp(IMR,JNP,IML,j1,j2,JN,JS,PU,DQ,Q,UC, & 967 967 fx1,xmass,IORD) 968 968 IMPLICIT NONE 969 integer:: IMR,JNP,IML,j1,j2,JN,JS,IORD970 real:: PU,DQ,Q,UC,fx1,xmass971 real:: dc,qtmp972 integer:: ISAVE(IMR)969 INTEGER :: IMR,JNP,IML,j1,j2,JN,JS,IORD 970 REAL :: PU,DQ,Q,UC,fx1,xmass 971 REAL :: dc,qtmp 972 INTEGER :: ISAVE(IMR) 973 973 dimension UC(IMR,*),DC(-IML:IMR+IML+1),xmass(IMR,JNP) & 974 974 ,fx1(IMR+1),DQ(IMR,JNP),qtmp(-IML:IMR+1+IML) 975 975 dimension PU(IMR,JNP),Q(IMR,JNP) 976 integer:: jvan,j1vl,j2vl,j,i,iu,itmp,ist,imp977 real:: rut976 INTEGER :: jvan,j1vl,j2vl,j,i,iu,itmp,ist,imp 977 REAL :: rut 978 978 ! 979 979 IMP = IMR + 1 … … 990 990 enddo 991 991 ! 992 if(j>=JN .or. j<=JS) goto 2222992 IF(j>=JN .or. j<=JS) goto 2222 993 993 ! ************* Eulerian ********** 994 994 ! … … 1007 1007 DC(0) = DC(IMR) 1008 1008 ! 1009 if(IORD==2 .or. j<=j1vl .or. j>=j2vl) then1009 IF(IORD==2 .or. j<=j1vl .or. j>=j2vl) THEN 1010 1010 DO i=1,IMR 1011 1011 iu = REAL(i) - uc(i,j) … … 1058 1058 ! 1059 1059 do i=1,IMR 1060 IF(uc(i,j)>1.) then1060 IF(uc(i,j)>1.) THEN 1061 1061 !DIR$ NOVECTOR 1062 1062 do ist = ISAVE(i),i-1 1063 1063 fx1(i) = fx1(i) + qtmp(ist) 1064 1064 enddo 1065 elseIF(uc(i,j)<-1.) then1065 elseIF(uc(i,j)<-1.) THEN 1066 1066 do ist = i,ISAVE(i)-1 1067 1067 fx1(i) = fx1(i) - qtmp(ist) … … 1084 1084 ! 1085 1085 END DO 1086 return1087 end subroutinextp1086 RETURN 1087 END SUBROUTINE xtp 1088 1088 ! 1089 1089 SUBROUTINE fxppm(IMR,IML,UT,P,DC,flux,IORD) 1090 1090 IMPLICIT NONE 1091 integer:: IMR,IML,IORD1092 real:: UT,P,DC,flux1091 INTEGER :: IMR,IML,IORD 1092 REAL :: UT,P,DC,flux 1093 1093 real,parameter :: R3 = 1./3., R23 = 2./3. 1094 1094 DIMENSION UT(*),flux(*),P(-IML:IMR+IML+1),DC(-IML:IMR+IML+1) 1095 1095 REAL :: AR(0:IMR),AL(0:IMR),A6(0:IMR) 1096 integer:: LMT,IMP,JLVL,i1096 INTEGER :: LMT,IMP,JLVL,i 1097 1097 ! logical first 1098 1098 ! data first /.TRUE./ 1099 1099 ! SAVE LMT 1100 ! if(first) then1100 ! IF(first) THEN 1101 1101 ! 1102 1102 ! correction calcul de LMT a chaque passage pour pouvoir choisir 1103 1103 ! plusieurs schemas PPM pour differents traceurs 1104 ! IF (IORD.LE.0) then1105 ! if(IMR.GE.144) then1104 ! IF (IORD.LE.0) THEN 1105 ! IF(IMR.GE.144) THEN 1106 1106 ! LMT = 0 1107 ! elseif(IMR.GE.72) then1107 ! elseif(IMR.GE.72) THEN 1108 1108 ! LMT = 1 1109 1109 ! else … … 1115 1115 ! 1116 1116 LMT = IORD - 3 1117 ! write(6,*) 'PPM option in E-W direction = ', LMT1117 ! WRITE(6,*) 'PPM option in E-W direction = ', LMT 1118 1118 ! first = .FALSE. 1119 1119 ! endif … … 1132 1132 END DO 1133 1133 ! 1134 if(LMT<=2) CALL lmtppm(DC(1),A6(1),AR(1),AL(1),P(1),IMR,LMT)1134 IF(LMT<=2) CALL lmtppm(DC(1),A6(1),AR(1),AL(1),P(1),IMR,LMT) 1135 1135 ! 1136 1136 AL(0) = AL(IMR) … … 1139 1139 ! 1140 1140 DO i=1,IMR 1141 IF(UT(i)>0.) then1141 IF(UT(i)>0.) THEN 1142 1142 flux(i) = AR(i-1) + 0.5*UT(i)*(AL(i-1) - AR(i-1) + & 1143 1143 A6(i-1)*(1.-R23*UT(i)) ) … … 1147 1147 endif 1148 1148 enddo 1149 return1150 end subroutinefxppm1149 RETURN 1150 END SUBROUTINE fxppm 1151 1151 ! 1152 1152 SUBROUTINE xmist(IMR,IML,P,DC) 1153 1153 IMPLICIT NONE 1154 integer:: IMR,IML1154 INTEGER :: IMR,IML 1155 1155 real,parameter :: R24 = 1./24. 1156 real:: P(-IML:IMR+1+IML),DC(-IML:IMR+1+IML)1157 integer:: i1158 real:: tmp,pmax,pmin1156 REAL :: P(-IML:IMR+1+IML),DC(-IML:IMR+1+IML) 1157 INTEGER :: i 1158 REAL :: tmp,pmax,pmin 1159 1159 ! 1160 1160 do i=1,IMR … … 1164 1164 DC(i) = sign(min(abs(tmp),Pmax,Pmin), tmp) 1165 1165 END DO 1166 return1167 end subroutinexmist1166 RETURN 1167 END SUBROUTINE xmist 1168 1168 ! 1169 1169 SUBROUTINE ytp(IMR,JNP,j1,j2,acosp,RCAP,DQ,P,VC,DC2 & 1170 1170 ,ymass,fx,A6,AR,AL,JORD) 1171 1171 IMPLICIT NONE 1172 integer:: IMR,JNP,j1,j2,JORD1173 real:: acosp,RCAP,DQ,P,VC,DC2,ymass,fx,A6,AR,AL1172 INTEGER :: IMR,JNP,j1,j2,JORD 1173 REAL :: acosp,RCAP,DQ,P,VC,DC2,ymass,fx,A6,AR,AL 1174 1174 dimension P(IMR,JNP),VC(IMR,JNP),ymass(IMR,JNP) & 1175 1175 ,DC2(IMR,JNP),DQ(IMR,JNP),acosp(JNP) 1176 1176 ! Work array 1177 1177 DIMENSION fx(IMR,JNP),AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP) 1178 integer:: JMR,len,i,jt,j1179 real:: sum1,sum21178 INTEGER :: JMR,len,i,jt,j 1179 REAL :: sum1,sum2 1180 1180 ! 1181 1181 JMR = JNP - 1 1182 1182 len = IMR*(J2-J1+2) 1183 1183 ! 1184 if(JORD==1) then1184 IF(JORD==1) THEN 1185 1185 DO i=1,len 1186 1186 JT = REAL(J1) - VC(i,J1) … … 1191 1191 CALL ymist(IMR,JNP,j1,P,DC2,4) 1192 1192 ! 1193 if(JORD<=0 .or. JORD>=3) then 1194 1193 IF(JORD<=0 .or. JORD>=3) THEN 1195 1194 CALL fyppm(VC,P,DC2,fx,IMR,JNP,j1,j2,A6,AR,AL,JORD) 1196 1195 … … 1228 1227 enddo 1229 1228 ! 1230 if(j1/=2) then1229 IF(j1/=2) THEN 1231 1230 do i=1,IMR 1232 1231 DQ(i, 2) = sum1 … … 1235 1234 endif 1236 1235 ! 1237 return1238 end subroutineytp1236 RETURN 1237 END SUBROUTINE ytp 1239 1238 ! 1240 1239 subroutine ymist(IMR,JNP,j1,P,DC,ID) 1241 1240 IMPLICIT NONE 1242 integer:: IMR,JNP,j1,ID1241 INTEGER :: IMR,JNP,j1,ID 1243 1242 real,parameter :: R24 = 1./24. 1244 real:: P(IMR,JNP),DC(IMR,JNP)1245 integer:: iimh,jmr,ijm3,imh,i1246 real:: pmax,pmin,tmp1243 REAL :: P(IMR,JNP),DC(IMR,JNP) 1244 INTEGER :: iimh,jmr,ijm3,imh,i 1245 REAL :: pmax,pmin,tmp 1247 1246 ! 1248 1247 IMH = IMR / 2 … … 1291 1290 ENDIF 1292 1291 ! 1293 if(j1/=2) then1292 IF(j1/=2) THEN 1294 1293 do i=1,IMR 1295 1294 DC(i,1) = 0. … … 1317 1316 END DO 1318 1317 endif 1319 return1320 end subroutineymist1318 RETURN 1319 END SUBROUTINE ymist 1321 1320 ! 1322 1321 SUBROUTINE fyppm(VC,P,DC,flux,IMR,JNP,j1,j2,A6,AR,AL,JORD) 1323 1322 IMPLICIT NONE 1324 integer:: IMR,JNP,j1,j2,JORD1323 INTEGER :: IMR,JNP,j1,j2,JORD 1325 1324 real,parameter :: R3 = 1./3., R23 = 2./3. 1326 real:: VC(IMR,*),flux(IMR,*),P(IMR,*),DC(IMR,*)1325 REAL :: VC(IMR,*),flux(IMR,*),P(IMR,*),DC(IMR,*) 1327 1326 ! Local work arrays. 1328 real:: AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP)1329 integer:: LMT,i1330 integer:: IMH,JMR,j11,IMJM1,len1327 REAL :: AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP) 1328 INTEGER :: LMT,i 1329 INTEGER :: IMH,JMR,j11,IMJM1,len 1331 1330 ! logical first 1332 1331 ! data first /.TRUE./ … … 1338 1337 IMJM1 = IMR*(J2-J1+2) 1339 1338 len = IMR*(J2-J1+3) 1340 ! if(first) then1341 ! IF(JORD.LE.0) then1342 ! if(JMR.GE.90) then1339 ! IF(first) THEN 1340 ! IF(JORD.LE.0) THEN 1341 ! IF(JMR.GE.90) THEN 1343 1342 ! LMT = 0 1344 ! elseif(JMR.GE.45) then1343 ! elseif(JMR.GE.45) THEN 1345 1344 ! LMT = 1 1346 1345 ! else … … 1384 1383 END DO 1385 1384 ! 1386 if(LMT<=2) CALL lmtppm(DC(1,j11),A6(1,j11),AR(1,j11) &1385 IF(LMT<=2) CALL lmtppm(DC(1,j11),A6(1,j11),AR(1,j11) & 1387 1386 ,AL(1,j11),P(1,j11),len,LMT) 1388 1387 ! 1389 1388 1390 1389 DO i=1,IMJM1 1391 IF(VC(i,j1)>0.) then1390 IF(VC(i,j1)>0.) THEN 1392 1391 flux(i,j1) = AR(i,j11) + 0.5*VC(i,j1)*(AL(i,j11) - AR(i,j11) + & 1393 1392 A6(i,j11)*(1.-R23*VC(i,j1)) ) … … 1397 1396 endif 1398 1397 END DO 1399 return1400 end subroutinefyppm1398 RETURN 1399 END SUBROUTINE fyppm 1401 1400 ! 1402 1401 SUBROUTINE yadv(IMR,JNP,j1,j2,p,VA,ady,wk,IAD) 1403 1402 IMPLICIT NONE 1404 integer:: IMR,JNP,j1,j2,IAD1403 INTEGER :: IMR,JNP,j1,j2,IAD 1405 1404 REAL :: p(IMR,JNP),ady(IMR,JNP),VA(IMR,JNP) 1406 1405 REAL :: WK(IMR,-1:JNP+2) … … 1426 1425 wk(i+IMH,JNP+2) = p(i,JNP-2) 1427 1426 enddo 1428 ! write(*,*) 'toto 1'1427 ! WRITE(*,*) 'toto 1' 1429 1428 ! -------------------------------- 1430 IF(IAD==2) then1429 IF(IAD==2) THEN 1431 1430 do j=j1-1,j2+1 1432 1431 do i=1,IMR 1433 ! write(*,*) 'avt NINT','i=',i,'j=',j1432 ! WRITE(*,*) 'avt NINT','i=',i,'j=',j 1434 1433 JP = NINT(VA(i,j)) 1435 1434 rv = JP - VA(i,j) 1436 ! write(*,*) 'VA=',VA(i,j), 'JP1=',JP,'rv=',rv1435 ! WRITE(*,*) 'VA=',VA(i,j), 'JP1=',JP,'rv=',rv 1437 1436 JP = j - JP 1438 ! write(*,*) 'JP2=',JP1437 ! WRITE(*,*) 'JP2=',JP 1439 1438 a1 = 0.5*(wk(i,jp+1)+wk(i,jp-1)) - wk(i,jp) 1440 1439 b1 = 0.5*(wk(i,jp+1)-wk(i,jp-1)) 1441 ! write(*,*) 'a1=',a1,'b1=',b11440 ! WRITE(*,*) 'a1=',a1,'b1=',b1 1442 1441 ady(i,j) = wk(i,jp) + rv*(a1*rv + b1) - wk(i,j) 1443 1442 enddo 1444 1443 enddo 1445 ! write(*,*) 'toto 2'1446 ! 1447 ELSEIF(IAD==1) then1444 ! WRITE(*,*) 'toto 2' 1445 ! 1446 ELSEIF(IAD==1) THEN 1448 1447 do j=j1-1,j2+1 1449 1448 do i=1,imr … … 1454 1453 ENDIF 1455 1454 ! 1456 if(j1/=2) then1455 IF(j1/=2) THEN 1457 1456 sum1 = 0. 1458 1457 sum2 = 0. … … 1487 1486 endif 1488 1487 ! 1489 return1490 end subroutineyadv1488 RETURN 1489 END SUBROUTINE yadv 1491 1490 ! 1492 1491 SUBROUTINE xadv(IMR,JNP,j1,j2,p,UA,JS,JN,IML,adx,IAD) … … 1499 1498 JMR = JNP-1 1500 1499 do j=j1,j2 1501 if(J>JS .and. J<JN) GO TO 13091500 IF(J>JS .and. J<JN) GO TO 1309 1502 1501 ! 1503 1502 do i=1,IMR … … 1519 1518 adx(i,j) = qtmp(ip) + ru*(a1*ru + b1) 1520 1519 enddo 1521 ELSEIF(IAD==1) then1520 ELSEIF(IAD==1) THEN 1522 1521 DO i=1,IMR 1523 1522 iu = UA(i,j) 1524 1523 ru = UA(i,j) - iu 1525 1524 iiu = i-iu 1526 if(UA(i,j)>=0.) then1525 IF(UA(i,j)>=0.) THEN 1527 1526 adx(i,j) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu)) 1528 1527 else … … 1560 1559 adx(i,j) = qtmp(ip)- p(i,j) + ru*(a1*ru + b1) 1561 1560 enddo 1562 ELSEIF(IAD==1) then1561 ELSEIF(IAD==1) THEN 1563 1562 ! 1st order 1564 1563 DO i=1,IMR … … 1569 1568 enddo 1570 1569 ! 1571 if(j1/=2) then1570 IF(j1/=2) THEN 1572 1571 do i=1,IMR 1573 1572 adx(i, 2) = 0. … … 1580 1579 adx(i,JNP) = 0. 1581 1580 enddo 1582 return1583 end subroutinexadv1581 RETURN 1582 END SUBROUTINE xadv 1584 1583 ! 1585 1584 SUBROUTINE lmtppm(DC,A6,AR,AL,P,IM,LMT) … … 1600 1599 ! 1601 1600 real,parameter :: R12 = 1./12. 1602 real:: A6(IM),AR(IM),AL(IM),P(IM),DC(IM)1603 integer:: IM,LMT1601 REAL :: A6(IM),AR(IM),AL(IM),P(IM),DC(IM) 1602 INTEGER :: IM,LMT 1604 1603 INTEGER :: i 1605 1604 REAL :: da1,da2,a6da,fmin 1606 1605 ! 1607 if(LMT==0) then1606 IF(LMT==0) THEN 1608 1607 ! Full constraint 1609 1608 do i=1,IM 1610 if(DC(i)==0.) then1609 IF(DC(i)==0.) THEN 1611 1610 AR(i) = p(i) 1612 1611 AL(i) = p(i) … … 1616 1615 da2 = da1**2 1617 1616 A6DA = A6(i)*da1 1618 if(A6DA < -da2) then1617 IF(A6DA < -da2) THEN 1619 1618 A6(i) = 3.*(AL(i)-p(i)) 1620 1619 AR(i) = AL(i) - A6(i) 1621 elseif(A6DA > da2) then1620 elseif(A6DA > da2) THEN 1622 1621 A6(i) = 3.*(AR(i)-p(i)) 1623 1622 AL(i) = AR(i) - A6(i) … … 1625 1624 endif 1626 1625 END DO 1627 elseif(LMT==1) then1626 elseif(LMT==1) THEN 1628 1627 ! Semi-monotonic constraint 1629 1628 do i=1,IM 1630 if(abs(AR(i)-AL(i)) >= -A6(i)) go to 1501631 if(p(i)<AR(i) .and. p(i)<AL(i)) then1629 IF(abs(AR(i)-AL(i)) >= -A6(i)) go to 150 1630 IF(p(i)<AR(i) .and. p(i)<AL(i)) THEN 1632 1631 AR(i) = p(i) 1633 1632 AL(i) = p(i) 1634 1633 A6(i) = 0. 1635 elseif(AR(i) > AL(i)) then1634 elseif(AR(i) > AL(i)) THEN 1636 1635 A6(i) = 3.*(AL(i)-p(i)) 1637 1636 AR(i) = AL(i) - A6(i) … … 1642 1641 150 continue 1643 1642 END DO 1644 elseif(LMT==2) then1643 elseif(LMT==2) THEN 1645 1644 do i=1,IM 1646 if(abs(AR(i)-AL(i)) >= -A6(i)) go to 2501645 IF(abs(AR(i)-AL(i)) >= -A6(i)) go to 250 1647 1646 fmin = p(i) + 0.25*(AR(i)-AL(i))**2/A6(i) + A6(i)*R12 1648 if(fmin>=0.) go to 2501649 if(p(i)<AR(i) .and. p(i)<AL(i)) then1647 IF(fmin>=0.) go to 250 1648 IF(p(i)<AR(i) .and. p(i)<AL(i)) THEN 1650 1649 AR(i) = p(i) 1651 1650 AL(i) = p(i) 1652 1651 A6(i) = 0. 1653 elseif(AR(i) > AL(i)) then1652 elseif(AR(i) > AL(i)) THEN 1654 1653 A6(i) = 3.*(AL(i)-p(i)) 1655 1654 AR(i) = AL(i) - A6(i) … … 1661 1660 END DO 1662 1661 endif 1663 return1664 end subroutinelmtppm1662 RETURN 1663 END SUBROUTINE lmtppm 1665 1664 ! 1666 1665 SUBROUTINE A2C(U,V,IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5) 1667 1666 IMPLICIT NONE 1668 integer:: IMR,JMR,j1,j21669 real:: U(IMR,*),V(IMR,*),CRX(IMR,*),CRY(IMR,*),DTDX5(*),DTDY51670 integer:: i,j1667 INTEGER :: IMR,JMR,j1,j2 1668 REAL :: U(IMR,*),V(IMR,*),CRX(IMR,*),CRY(IMR,*),DTDX5(*),DTDY5 1669 INTEGER :: i,j 1671 1670 ! 1672 1671 do j=j1,j2 … … 1683 1682 CRY(i,2) = DTDY5*(V(i,2)+V(i,1)) 1684 1683 END DO 1685 return1686 end subroutinea2c1684 RETURN 1685 END SUBROUTINE a2c 1687 1686 ! 1688 1687 SUBROUTINE cosa(cosp,cose,JNP,PI,DP) 1689 1688 IMPLICIT NONE 1690 integer:: JNP1691 real:: cosp(*),cose(*),PI,DP1692 integer:: JMR,j,jeq1693 real:: ph51689 INTEGER :: JNP 1690 REAL :: cosp(*),cose(*),PI,DP 1691 INTEGER :: JMR,j,jeq 1692 REAL :: ph5 1694 1693 JMR = JNP-1 1695 1694 do j=2,JNP … … 1699 1698 ! 1700 1699 JEQ = (JNP+1) / 2 1701 if(JMR == 2*(JMR/2) ) then1700 IF(JMR == 2*(JMR/2) ) THEN 1702 1701 do j=JNP, JEQ+1, -1 1703 1702 cose(j) = cose(JNP+2-j) … … 1716 1715 cosp(1) = 0. 1717 1716 cosp(JNP) = 0. 1718 return1719 end subroutinecosa1717 RETURN 1718 END SUBROUTINE cosa 1720 1719 ! 1721 1720 SUBROUTINE cosc(cosp,cose,JNP,PI,DP) 1722 1721 IMPLICIT NONE 1723 integer:: JNP1724 real:: cosp(*),cose(*),PI,DP1725 real:: phi1726 integer:: j1722 INTEGER :: JNP 1723 REAL :: cosp(*),cose(*),PI,DP 1724 REAL :: phi 1725 INTEGER :: j 1727 1726 ! 1728 1727 phi = -0.5*PI … … 1741 1740 cosp(j) = 0.5*(cose(j)+cose(j+1)) 1742 1741 END DO 1743 return1744 end subroutinecosc1742 RETURN 1743 END SUBROUTINE cosc 1745 1744 ! 1746 1745 SUBROUTINE qckxyz(Q,qtmp,IMR,JNP,NLAY,j1,j2,cosp,acosp, & … … 1752 1751 logical :: cross 1753 1752 INTEGER :: NLAYM1,len,ip,L,icr,ipy,ipx,i 1754 real:: qup,qly,dup,sum1753 REAL :: qup,qly,dup,sum 1755 1754 ! 1756 1755 NLAYM1 = NLAY-1 … … 1762 1761 icr = 1 1763 1762 CALL filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny) 1764 if(ipy==0) goto 501763 IF(ipy==0) goto 50 1765 1764 CALL filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny) 1766 if(ipx==0) goto 501767 ! 1768 if(cross) then1765 IF(ipx==0) goto 50 1766 ! 1767 IF(cross) THEN 1769 1768 CALL filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny) 1770 1769 endif 1771 if(icr==0) goto 501770 IF(icr==0) goto 50 1772 1771 ! 1773 1772 ! Vertical filling... … … 1785 1784 ! 1786 1785 CALL filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny) 1787 if(ipy==0) goto 2251786 IF(ipy==0) goto 225 1788 1787 CALL filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny) 1789 if(ipx==0) go to 2251790 if(cross) then1788 IF(ipx==0) go to 225 1789 IF(cross) THEN 1791 1790 CALL filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny) 1792 1791 endif 1793 if(icr==0) goto 2251792 IF(icr==0) goto 225 1794 1793 ! 1795 1794 do i=1,len … … 1816 1815 ! 1817 1816 CALL filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny) 1818 if(ipy==0) goto 9111817 IF(ipy==0) goto 911 1819 1818 CALL filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny) 1820 if(ipx==0) goto 9111819 IF(ipx==0) goto 911 1821 1820 ! 1822 1821 CALL filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny) 1823 if(icr==0) goto 9111822 IF(icr==0) goto 911 1824 1823 ! 1825 1824 DO I=1,len … … 1841 1840 911 continue 1842 1841 ! 1843 if(ip>IMR) then1844 write(6,*) 'IC=',IC,' STEP=',NSTEP, &1842 IF(ip>IMR) THEN 1843 WRITE(6,*) 'IC=',IC,' STEP=',NSTEP, & 1845 1844 ' Vertical filling pts=',ip 1846 1845 endif 1847 1846 ! 1848 if(sum>1.e-25) then1849 write(6,*) IC,NSTEP,' Mass source from the ground=',sum1847 IF(sum>1.e-25) THEN 1848 WRITE(6,*) IC,NSTEP,' Mass source from the ground=',sum 1850 1849 endif 1851 1850 RETURN … … 1854 1853 SUBROUTINE filcr(q,IMR,JNP,j1,j2,cosp,acosp,icr,tiny) 1855 1854 IMPLICIT NONE 1856 integer:: IMR,JNP,j1,j2,icr1857 real:: q(IMR,*),cosp(*),acosp(*),tiny1858 integer:: i,j1859 real:: dq,dn,d0,d1,ds,d21855 INTEGER :: IMR,JNP,j1,j2,icr 1856 REAL :: q(IMR,*),cosp(*),acosp(*),tiny 1857 INTEGER :: i,j 1858 REAL :: dq,dn,d0,d1,ds,d2 1860 1859 icr = 0 1861 1860 do j=j1+1,j2-1 … … 1878 1877 endif 1879 1878 END DO 1880 if(icr==0 .and. q(IMR,j)>=0.) goto 651879 IF(icr==0 .and. q(IMR,j)>=0.) goto 65 1881 1880 DO i=2,IMR 1882 1881 IF(q(i,j)<0.) THEN … … 1940 1939 ! 1941 1940 do i=1,IMR 1942 if(q(i,j1)<0. .or. q(i,j2)<0.) then1941 IF(q(i,j1)<0. .or. q(i,j2)<0.) THEN 1943 1942 icr = 1 1944 1943 goto 80 … … 1948 1947 80 continue 1949 1948 ! 1950 if(q(1,1)<0. .or. q(1,jnp)<0.) then1949 IF(q(1,1)<0. .or. q(1,jnp)<0.) THEN 1951 1950 icr = 1 1952 1951 endif 1953 1952 ! 1954 return1955 end subroutinefilcr1953 RETURN 1954 END SUBROUTINE filcr 1956 1955 ! 1957 1956 SUBROUTINE filns(q,IMR,JNP,j1,j2,cosp,acosp,ipy,tiny) 1958 1957 IMPLICIT NONE 1959 integer:: IMR,JNP,j1,j2,ipy1960 real:: q(IMR,*),cosp(*),acosp(*),tiny1961 real:: DP,CAP1,dq,dn,d0,d1,ds,d21958 INTEGER :: IMR,JNP,j1,j2,ipy 1959 REAL :: q(IMR,*),cosp(*),acosp(*),tiny 1960 REAL :: DP,CAP1,dq,dn,d0,d1,ds,d2 1962 1961 INTEGER :: i,j 1963 1962 ! logical first … … 1965 1964 ! save cap1 1966 1965 ! 1967 ! if(first) then1966 ! IF(first) THEN 1968 1967 DP = 4.*ATAN(1.)/REAL(JNP-1) 1969 1968 CAP1 = IMR*(1.-COS((j1-1.5)*DP))/DP … … 2021 2020 ! 2022 2021 ! Check Poles. 2023 if(q(1,1)<0.) then2022 IF(q(1,1)<0.) THEN 2024 2023 dq = q(1,1)*cap1/REAL(IMR)*acosp(j1) 2025 2024 do i=1,imr 2026 2025 q(i,1) = 0. 2027 2026 q(i,j1) = q(i,j1) + dq 2028 if(q(i,j1)<0.) ipy = 12029 enddo 2030 endif 2031 ! 2032 if(q(1,JNP)<0.) then2027 IF(q(i,j1)<0.) ipy = 1 2028 enddo 2029 endif 2030 ! 2031 IF(q(1,JNP)<0.) THEN 2033 2032 dq = q(1,JNP)*cap1/REAL(IMR)*acosp(j2) 2034 2033 do i=1,imr 2035 2034 q(i,JNP) = 0. 2036 2035 q(i,j2) = q(i,j2) + dq 2037 if(q(i,j2)<0.) ipy = 12038 enddo 2039 endif 2040 ! 2041 return2042 end subroutinefilns2036 IF(q(i,j2)<0.) ipy = 1 2037 enddo 2038 endif 2039 ! 2040 RETURN 2041 END SUBROUTINE filns 2043 2042 ! 2044 2043 SUBROUTINE filew(q,qtmp,IMR,JNP,j1,j2,ipx,tiny) 2045 2044 IMPLICIT NONE 2046 integer:: IMR,JNP,j1,j2,ipx2047 real:: q(IMR,*),qtmp(JNP,IMR),tiny2048 integer:: i,j2049 real:: d0,d1,d22045 INTEGER :: IMR,JNP,j1,j2,ipx 2046 REAL :: q(IMR,*),qtmp(JNP,IMR),tiny 2047 INTEGER :: i,j 2048 REAL :: d0,d1,d2 2050 2049 ! 2051 2050 ipx = 0 … … 2059 2058 do i=2,imr-1 2060 2059 do j=j1,j2 2061 if(qtmp(j,i)<0.) then2060 IF(qtmp(j,i)<0.) THEN 2062 2061 ipx = 1 2063 2062 ! west … … 2077 2076 i=1 2078 2077 do j=j1,j2 2079 if(qtmp(j,i)<0.) then2078 IF(qtmp(j,i)<0.) THEN 2080 2079 ipx = 1 2081 2080 ! west … … 2094 2093 i=IMR 2095 2094 do j=j1,j2 2096 if(qtmp(j,i)<0.) then2095 IF(qtmp(j,i)<0.) THEN 2097 2096 ipx = 1 2098 2097 ! west … … 2110 2109 END DO 2111 2110 ! 2112 if(ipx/=0) then2111 IF(ipx/=0) THEN 2113 2112 do j=j1,j2 2114 2113 do i=1,imr … … 2119 2118 ! 2120 2119 ! Poles. 2121 if(q(1,1)<0 .or. q(1,JNP)<0.) ipx = 12122 endif 2123 return2124 end subroutinefilew2120 IF(q(1,1)<0 .or. q(1,JNP)<0.) ipx = 1 2121 endif 2122 RETURN 2123 END SUBROUTINE filew 2125 2124 ! 2126 2125 SUBROUTINE zflip(q,im,km,nc) 2127 2126 IMPLICIT NONE 2128 2127 ! This routine flip the array q (in the vertical). 2129 integer:: im,km,nc2130 real:: q(im,km,nc)2128 INTEGER :: im,km,nc 2129 REAL :: q(im,km,nc) 2131 2130 ! local dynamic array 2132 real:: qtmp(im,km)2133 integer:: IC,k,i2131 REAL :: qtmp(im,km) 2132 INTEGER :: IC,k,i 2134 2133 ! 2135 2134 do IC = 1, nc … … 2145 2144 END DO 2146 2145 END DO 2147 return2148 end subroutinezflip2146 RETURN 2147 END SUBROUTINE zflip -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/prather.f90
r5106 r5116 33 33 REAL :: q( iip1,jjp1,llm,0:9) 34 34 REAL :: w( ip1jmp1,llm ) 35 integer:: ordre,ilim35 INTEGER :: ordre,ilim 36 36 37 37 ! Local: 38 38 ! ------ 39 39 LOGICAL :: limit 40 real:: zq(iip1,jjp1,llm)40 REAL :: zq(iip1,jjp1,llm) 41 41 REAL :: sm ( iip1,jjp1, llm ) 42 42 REAL :: s0( iip1,jjp1,llm ), sx( iip1,jjp1,llm ) … … 49 49 REAL :: szz( iip1,jjp1,llm ),zz 50 50 INTEGER :: i,j,l,indice 51 real:: sxn(iip1),sxs(iip1)52 53 real:: sinlon(iip1),sinlondlon(iip1)54 real:: coslon(iip1),coslondlon(iip1)55 real:: qmin,qmax51 REAL :: sxn(iip1),sxs(iip1) 52 53 REAL :: sinlon(iip1),sinlondlon(iip1) 54 REAL :: coslon(iip1),coslondlon(iip1) 55 REAL :: qmin,qmax 56 56 save qmin,qmax 57 57 save sinlon,coslon,sinlondlon,coslondlon 58 real:: dyn1,dyn2,dys1,dys2,qpn,qps,dqzpn,dqzps59 real:: masn,mass58 REAL :: dyn1,dyn2,dys1,dys2,qpn,qps,dqzpn,dqzps 59 REAL :: masn,mass 60 60 ! 61 61 REAL :: SSUM 62 integer:: ismax,ismin62 INTEGER :: ismax,ismin 63 63 EXTERNAL SSUM, ismin,ismax 64 64 logical :: first … … 80 80 limit = .TRUE. 81 81 82 if(first) then82 IF(first) THEN 83 83 PRINT*,'SCHEMA PRATHER' 84 84 first=.FALSE. -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/principal_cshift_m.F90
r5113 r5116 11 11 ! xprimm. 12 12 13 use nrtype, only: twopi14 use serre_mod, only: clon13 use nrtype, ONLY: twopi 14 use serre_mod, ONLY: clon 15 15 16 16 include "dimensions.h" … … 22 22 !----------------------------------------------------- 23 23 24 if (is2 /= 0) then24 if (is2 /= 0) THEN 25 25 IF (clon <= 0.) THEN 26 26 IF (is2 /= 1) THEN -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/sortvarc.f90
r5113 r5116 70 70 ! Ehouarn: when no initialization fields from file, resetvarc should be 71 71 ! set to false 72 if (firstcal) then73 if (.not.read_start) then72 if (firstcal) THEN 73 if (.not.read_start) THEN 74 74 resetvarc=.TRUE. 75 75 endif … … 147 147 ang = SSUM( llm, angl, 1 ) 148 148 149 IF (firstcal.and.resetvarc) then149 IF (firstcal.and.resetvarc) THEN 150 150 WRITE(lunout,3500) itau, rjour, heure, time 151 151 WRITE(lunout,*) trim(modname), & … … 162 162 ! compute relative changes in etot,... (except if 'reference' values 163 163 ! are zero, which can happen when using iniacademic) 164 if (etot0/=0) then164 if (etot0/=0) THEN 165 165 etot= etot/etot0 166 166 else … … 168 168 endif 169 169 rmsv= SQRT(rmsv/ptot) 170 if (ptot0/=0) then170 if (ptot0/=0) THEN 171 171 ptot= ptot/ptot0 172 172 else 173 173 ptot=1. 174 174 endif 175 if (ztot0/=0) then175 if (ztot0/=0) THEN 176 176 ztot= ztot/ztot0 177 177 else 178 178 ztot=1. 179 179 endif 180 if (stot0/=0) then180 if (stot0/=0) THEN 181 181 stot= stot/stot0 182 182 else 183 183 stot=1. 184 184 endif 185 if (ang0/=0) then185 if (ang0/=0) THEN 186 186 ang = ang /ang0 187 187 else -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/test_period.f90
r5106 r5116 44 44 do ij=1,iim 45 45 if (teta(ij,l)/=teta(1,l) & 46 .or.teta(ip1jm+ij,l)/=teta(ip1jm+1,l) ) then46 .or.teta(ip1jm+ij,l)/=teta(ip1jm+1,l) ) THEN 47 47 PRINT *,'STOP dans test_period car --- TETA --- n est pas', & 48 48 ' constant aux poles ! ' … … 100 100 do ij=1,iim 101 101 if (p(ij,l)/=p(1,l) & 102 .or.p(ip1jm+ij,l)/=p(ip1jm+1,l) ) then102 .or.p(ip1jm+ij,l)/=p(ip1jm+1,l) ) THEN 103 103 PRINT *,'STOP dans test_period car --- P --- n est pas', & 104 104 ' constant aux poles ! ' -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/traceurpole.f90
r5114 r5116 12 12 13 13 ! Arguments 14 integer:: iq15 real:: masse(iip1, jjp1, llm)16 real:: q(iip1, jjp1, llm)14 INTEGER :: iq 15 REAL :: masse(iip1, jjp1, llm) 16 REAL :: q(iip1, jjp1, llm) 17 17 18 18 19 19 ! Locals 20 integer:: i, j, l21 real:: sommemassen(llm)22 real:: sommemqn(llm)23 real:: sommemasses(llm)24 real:: sommemqs(llm)25 real:: qpolen(llm), qpoles(llm)20 INTEGER :: i, j, l 21 REAL :: sommemassen(llm) 22 REAL :: sommemqn(llm) 23 REAL :: sommemasses(llm) 24 REAL :: sommemqs(llm) 25 REAL :: qpolen(llm), qpoles(llm) 26 26 27 27 … … 56 56 enddo 57 57 58 return59 end subroutinetraceurpole58 RETURN 59 END SUBROUTINE traceurpole -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/ugeostr.F90
r5113 r5116 11 11 ! levels are pressure levels. 12 12 13 use comconst_mod, only: omeg, rad13 use comconst_mod, ONLY: omeg, rad 14 14 15 15 IMPLICIT NONE … … 29 29 DO j=1,jjm 30 30 31 if (abs(sin(rlatv(j)))<1.e-4) then31 if (abs(sin(rlatv(j)))<1.e-4) THEN 32 32 zlat=1.e-4 33 33 else -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/write_grads_dyn.h
r5105 r5116 2 2 ! $Header$ 3 3 4 if (callinigrads) then 5 4 if (callinigrads) THEN 6 5 string10='dyn' 7 6 CALL inigrads(1,iip1 & … … 26 25 do iq=1,nqtot 27 26 string10='q' 28 write(string10(2:2),'(i1)') iq27 WRITE(string10(2:2),'(i1)') iq 29 28 CALL wrgrads(1,llm,q(:,:,iq),string10,string10) 30 29 enddo -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/writedynav.F90
r5114 r5116 120 120 ! CALL histwrite(histaveid, 'phis', itau_w, phis, iip1*jjp1, ndex2d) 121 121 122 if (ok_sync) then122 if (ok_sync) THEN 123 123 CALL histsync(histaveid) 124 124 CALL histsync(histvaveid) -
LMDZ6/branches/Amaury_dev/libf/dyn3d_common/writehist.f90
r5114 r5116 47 47 REAL :: phis(ip1jmp1) 48 48 REAL :: q(ip1jmp1, llm, nqtot) 49 integer:: time49 INTEGER :: time 50 50 51 51 … … 53 53 ! Variables locales 54 54 ! 55 integer:: iq, ii, ll56 integer:: ndexu(ip1jmp1 * llm), ndexv(ip1jm * llm), ndex2d(ip1jmp1)55 INTEGER :: iq, ii, ll 56 INTEGER :: ndexu(ip1jmp1 * llm), ndexv(ip1jm * llm), ndex2d(ip1jmp1) 57 57 logical :: ok_sync 58 integer:: itau_w58 INTEGER :: itau_w 59 59 REAL :: vnat(ip1jm, llm), unat(ip1jmp1, llm) 60 60 … … 114 114 ! Fin 115 115 ! 116 if (ok_sync) then116 if (ok_sync) THEN 117 117 CALL histsync(histid) 118 118 CALL histsync(histvid) 119 119 CALL histsync(histuid) 120 120 endif 121 return122 end subroutinewritehist121 RETURN 122 END SUBROUTINE writehist
Note: See TracChangeset
for help on using the changeset viewer.
