Changeset 109 for trunk/libf
- Timestamp:
- Apr 14, 2011, 11:47:04 AM (14 years ago)
- Location:
- trunk/libf
- Files:
-
- 2 added
- 1 deleted
- 13 edited
- 4 copied
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
trunk/libf/dyn3d/comvert.h
r1 r109 6 6 7 7 COMMON/comvert/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm), & 8 & pa,preff,nivsigs(llm),nivsig(llm+1) 8 & pa,preff,nivsigs(llm),nivsig(llm+1), & 9 & aps(llm),bps(llm) 9 10 10 REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig 11 REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig,aps,bps 11 12 12 13 !----------------------------------------------------------------------- -
trunk/libf/dyn3d/disvert_terre.F90
r107 r109 1 1 ! $Id: disvert.F90 1480 2011-01-31 21:29:58Z jghattas $ 2 2 3 SUBROUTINE disvert (pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)3 SUBROUTINE disvert_terre(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig) 4 4 5 5 ! Auteur : P. Le Van -
trunk/libf/dyn3d/exner_milieu.F
r107 r109 1 SUBROUTINE exner_ hyb( ngrid, ps, p,beta, pks, pk, pkf )1 SUBROUTINE exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf ) 2 2 c 3 3 c Auteurs : F. Forget , Y. Wanherdrick … … 17 17 c 18 18 c WARNING : CECI est une version speciale de exner_hyb originale 19 c Utilis dans la version martienne pour pouvoir20 c tourner avec des coordonn es verticales complexe21 c => Il ne verifie PAS la condition la proportionalit en22 c nergie totale/ interne / potentielle (F.Forget 2001)19 c Utilise dans la version martienne pour pouvoir 20 c tourner avec des coordonnees verticales complexe 21 c => Il ne verifie PAS la condition la proportionalite en 22 c energie totale/ interne / potentielle (F.Forget 2001) 23 23 c ( voir note de Fr.Hourdin ) , 24 24 c -
trunk/libf/dyn3d/iniconst.F
r1 r109 53 53 c----------------------------------------------------------------------- 54 54 55 CALL disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig) 56 c 55 if (planet_type.eq."earth") then 56 CALL disvert_terre(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig) 57 else 58 CALL disvert_noterre 59 endif 57 60 c 58 61 RETURN -
trunk/libf/dyn3d/leapfrog.F
r108 r109 238 238 dq(:,:,:)=0. 239 239 CALL pression ( ip1jmp1, ap, bp, ps, p ) 240 CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) 240 if (planet_type.eq."earth") then 241 CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) 242 else 243 CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf ) 244 endif 245 241 246 c------------------ 242 247 c TEST PK MONOTONE … … 404 409 405 410 CALL pression ( ip1jmp1, ap, bp, ps, p ) 406 CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf ) 411 if (planet_type.eq."earth") then 412 CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) 413 else 414 CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf ) 415 endif 407 416 408 417 ! rdaym_ini = itau * dtvr / daysec … … 519 528 520 529 CALL pression ( ip1jmp1, ap, bp, ps, p ) 521 CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) 530 if (planet_type.eq."earth") then 531 CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) 532 else 533 CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf ) 534 endif 522 535 523 536 -
trunk/libf/dyn3d/limy.F
r1 r109 40 40 REAL qbyv(ip1jm,llm) 41 41 42 REAL qpns,qpsn,ap n,aps,dyn1,dys1,dyn2,dys242 REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2 43 43 Logical extremum,first 44 44 save first … … 117 117 118 118 c print*,dyqv(iip1+1) 119 c ap n=abs(dyq(1)/dyqv(iip1+1))119 c appn=abs(dyq(1)/dyqv(iip1+1)) 120 120 c print*,dyq(ip1jm+1) 121 121 c print*,dyqv(ip1jm-iip1+1) 122 c ap s=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))122 c apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) 123 123 c do ij=2,iim 124 c ap n=amax1(abs(dyq(ij)/dyqv(ij)),apn)125 c ap s=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)124 c appn=amax1(abs(dyq(ij)/dyqv(ij)),appn) 125 c apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps) 126 126 c enddo 127 c ap n=min(pente_max/apn,1.)128 c ap s=min(pente_max/aps,1.)127 c appn=min(pente_max/appn,1.) 128 c apps=min(pente_max/apps,1.) 129 129 130 130 … … 132 132 133 133 c if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.) 134 c & ap n=0.134 c & appn=0. 135 135 c if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* 136 136 c & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.) 137 c & ap s=0.137 c & apps=0. 138 138 139 139 c limitation des pentes aux poles 140 140 c do ij=1,iip1 141 c dyq(ij)=ap n*dyq(ij)142 c dyq(ip1jm+ij)=ap s*dyq(ip1jm+ij)141 c dyq(ij)=appn*dyq(ij) 142 c dyq(ip1jm+ij)=apps*dyq(ip1jm+ij) 143 143 c enddo 144 144 -
trunk/libf/dyn3d/vlsplt.F
r1 r109 478 478 REAL qbyv(ip1jm,llm) 479 479 480 REAL qpns,qpsn,ap n,aps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs480 REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs 481 481 c REAL newq,oldmasse 482 482 Logical extremum,first,testcpu … … 602 602 C PRINT*,dyq(1) 603 603 C PRINT*,dyqv(iip1+1) 604 C ap n=abs(dyq(1)/dyqv(iip1+1))604 C appn=abs(dyq(1)/dyqv(iip1+1)) 605 605 C PRINT*,dyq(ip1jm+1) 606 606 C PRINT*,dyqv(ip1jm-iip1+1) 607 C ap s=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))607 C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) 608 608 C DO ij=2,iim 609 C ap n=amax1(abs(dyq(ij)/dyqv(ij)),apn)610 C ap s=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)609 C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn) 610 C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps) 611 611 C ENDDO 612 C ap n=min(pente_max/apn,1.)613 C ap s=min(pente_max/aps,1.)612 C appn=min(pente_max/appn,1.) 613 C apps=min(pente_max/apps,1.) 614 614 C 615 615 C … … 617 617 C 618 618 C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.) 619 C & ap n=0.619 C & appn=0. 620 620 C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* 621 621 C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.) 622 C & ap s=0.622 C & apps=0. 623 623 C 624 624 C limitation des pentes aux poles 625 625 C DO ij=1,iip1 626 C dyq(ij)=ap n*dyq(ij)627 C dyq(ip1jm+ij)=ap s*dyq(ip1jm+ij)626 C dyq(ij)=appn*dyq(ij) 627 C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij) 628 628 C ENDDO 629 629 C -
trunk/libf/dyn3d/vlspltqs.F
r5 r109 635 635 C PRINT*,dyq(1) 636 636 C PRINT*,dyqv(iip1+1) 637 C ap n=abs(dyq(1)/dyqv(iip1+1))637 C appn=abs(dyq(1)/dyqv(iip1+1)) 638 638 C PRINT*,dyq(ip1jm+1) 639 639 C PRINT*,dyqv(ip1jm-iip1+1) 640 C ap s=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))640 C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) 641 641 C DO ij=2,iim 642 C ap n=amax1(abs(dyq(ij)/dyqv(ij)),apn)643 C ap s=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)642 C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn) 643 C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps) 644 644 C ENDDO 645 C ap n=min(pente_max/apn,1.)646 C ap s=min(pente_max/aps,1.)645 C appn=min(pente_max/appn,1.) 646 C apps=min(pente_max/apps,1.) 647 647 C 648 648 C … … 650 650 C 651 651 C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.) 652 C & ap n=0.652 C & appn=0. 653 653 C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* 654 654 C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.) 655 C & ap s=0.655 C & apps=0. 656 656 C 657 657 C limitation des pentes aux poles 658 658 C DO ij=1,iip1 659 C dyq(ij)=ap n*dyq(ij)660 C dyq(ip1jm+ij)=ap s*dyq(ip1jm+ij)659 C dyq(ij)=appn*dyq(ij) 660 C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij) 661 661 C ENDDO 662 662 C -
trunk/libf/dyn3dpar/comvert.h
r1 r109 1 1 ! 2 ! $ Header$2 ! $Id: comvert.h 1279 2009-12-10 09:02:56Z fairhead $ 3 3 ! 4 4 !----------------------------------------------------------------------- 5 5 ! INCLUDE 'comvert.h' 6 6 7 COMMON/comvert/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm), & 8 & pa,preff,nivsigs(llm),nivsig(llm+1) 7 COMMON/comvert/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm), & 8 & pa,preff,nivsigs(llm),nivsig(llm+1), & 9 & aps(llm),bps(llm) 9 10 10 REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig 11 REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig,aps,bps 11 12 12 !-----------------------------------------------------------------------13 !----------------------------------------------------------------------- -
trunk/libf/dyn3dpar/disvert_terre.F90
r107 r109 1 1 ! $Id: disvert.F90 1480 2011-01-31 21:29:58Z jghattas $ 2 2 3 SUBROUTINE disvert (pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)3 SUBROUTINE disvert_terre(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig) 4 4 5 5 ! Auteur : P. Le Van -
trunk/libf/dyn3dpar/exner_milieu.F
r107 r109 1 SUBROUTINE exner_ hyb( ngrid, ps, p,beta, pks, pk, pkf )1 SUBROUTINE exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf ) 2 2 c 3 3 c Auteurs : F. Forget , Y. Wanherdrick … … 17 17 c 18 18 c WARNING : CECI est une version speciale de exner_hyb originale 19 c Utilis dans la version martienne pour pouvoir20 c tourner avec des coordonn es verticales complexe21 c => Il ne verifie PAS la condition la proportionalit en22 c nergie totale/ interne / potentielle (F.Forget 2001)19 c Utilise dans la version martienne pour pouvoir 20 c tourner avec des coordonnees verticales complexe 21 c => Il ne verifie PAS la condition la proportionalite en 22 c energie totale/ interne / potentielle (F.Forget 2001) 23 23 c ( voir note de Fr.Hourdin ) , 24 24 c -
trunk/libf/dyn3dpar/exner_milieu_p.F
r107 r109 1 SUBROUTINE exner_ hyb( ngrid, ps, p,beta, pks, pk, pkf )1 SUBROUTINE exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf ) 2 2 c 3 3 c Auteurs : F. Forget , Y. Wanherdrick … … 17 17 c 18 18 c WARNING : CECI est une version speciale de exner_hyb originale 19 c Utilis dans la version martienne pour pouvoir20 c tourner avec des coordonn es verticales complexe21 c => Il ne verifie PAS la condition la proportionalit en22 c nergie totale/ interne / potentielle (F.Forget 2001)19 c Utilise dans la version martienne pour pouvoir 20 c tourner avec des coordonnees verticales complexe 21 c => Il ne verifie PAS la condition la proportionalite en 22 c energie totale/ interne / potentielle (F.Forget 2001) 23 23 c ( voir note de Fr.Hourdin ) , 24 24 c 25 USE parallel 25 26 IMPLICIT NONE 26 27 c … … 44 45 REAL xpn, xps 45 46 REAL SSUM 46 EXTERNAL filtreg, SSUM 47 EXTERNAL SSUM 48 INTEGER ije,ijb,jje,jjb 47 49 50 c$OMP BARRIER 51 48 52 c ------------- 49 53 c Calcul de pks 50 54 c ------------- 51 55 52 DO ij = 1, ngrid 56 ijb=ij_begin 57 ije=ij_end 58 59 c$OMP DO SCHEDULE(STATIC) 60 DO ij = ijb, ije 53 61 pks(ij) = cpp * ( ps(ij)/preff ) ** kappa 54 62 ENDDO 63 c$OMP ENDDO 64 c Synchro OPENMP ici 55 65 56 DO ij = 1, iim 57 ppn(ij) = aire( ij ) * pks( ij ) 58 pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm ) 59 ENDDO 60 xpn = SSUM(iim,ppn,1) /apoln 61 xps = SSUM(iim,pps,1) /apols 62 63 DO ij = 1, iip1 64 pks( ij ) = xpn 65 pks( ij+ip1jm ) = xps 66 ENDDO 66 c$OMP MASTER 67 if (pole_nord) then 68 DO ij = 1, iim 69 ppn(ij) = aire( ij ) * pks( ij ) 70 ENDDO 71 xpn = SSUM(iim,ppn,1) /apoln 72 73 DO ij = 1, iip1 74 pks( ij ) = xpn 75 ENDDO 76 endif 77 78 if (pole_sud) then 79 DO ij = 1, iim 80 pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm ) 81 ENDDO 82 xps = SSUM(iim,pps,1) /apols 83 84 DO ij = 1, iip1 85 pks( ij+ip1jm ) = xps 86 ENDDO 87 endif 88 c$OMP END MASTER 67 89 c 68 90 c … … 72 94 dum1 = cpp * (2*preff)**(-kappa) 73 95 DO l = 1, llm-1 74 DO ij = 1, ngrid 96 c$OMP DO SCHEDULE(STATIC) 97 DO ij = ijb, ije 75 98 pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa 76 99 ENDDO 100 c$OMP ENDDO NOWAIT 77 101 ENDDO 78 102 … … 81 105 c et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2) 82 106 83 DO ij = 1, ngrid 107 c$OMP DO SCHEDULE(STATIC) 108 DO ij = ijb, ije 84 109 pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2) 85 110 ENDDO 111 c$OMP ENDDO NOWAIT 86 112 87 113 88 114 c calcul de pkf 89 115 c ------------- 90 CALL SCOPY ( ngrid * llm, pk, 1, pkf, 1 ) 91 CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 116 c CALL SCOPY ( ngrid * llm, pk, 1, pkf, 1 ) 117 DO l = 1, llm 118 c$OMP DO SCHEDULE(STATIC) 119 DO ij = ijb, ije 120 pkf(ij,l)=pk(ij,l) 121 ENDDO 122 c$OMP ENDDO NOWAIT 123 ENDDO 124 125 c$OMP BARRIER 126 127 jjb=jj_begin 128 jje=jj_end 129 CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 ) 92 130 93 131 c EST-CE UTILE ?? : calcul de beta 94 132 c -------------------------------- 95 133 DO l = 2, llm 96 DO ij = 1, ngrid 134 c$OMP DO SCHEDULE(STATIC) 135 DO ij = ijb, ije 97 136 beta(ij,l) = pk(ij,l) / pk(ij,l-1) 98 137 ENDDO 138 c$OMP ENDDO NOWAIT 99 139 ENDDO 100 140 -
trunk/libf/dyn3dpar/iniconst.F
r1 r109 53 53 c----------------------------------------------------------------------- 54 54 55 CALL disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig) 56 c 55 if (planet_type.eq."earth") then 56 CALL disvert_terre(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig) 57 else 58 CALL disvert_noterre 59 endif 57 60 c 58 61 RETURN -
trunk/libf/dyn3dpar/inidissip.F
r1 r109 31 31 INTEGER l,ij,idum,ii 32 32 REAL tetamin 33 REAL pseudoz33 REAL Pup 34 34 35 35 REAL ran1 … … 177 177 c -------------------------------------------------- 178 178 179 if (ok_strato .and. llm==39) then 180 do l=1,llm 181 pseudoz=8.*log(preff/presnivs(l)) 182 zvert(l)=1+ 183 s (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2. 184 s *(dissip_factz-1.) 185 enddo 186 else 187 DO l=1,llm 188 zvert(l)=1. 189 ENDDO 190 fact=2. 191 DO l = 1, llm 192 zz = 1. - preff/presnivs(l) 193 zvert(l)= fact -( fact-1.)/( 1.+zz*zz ) 194 ENDDO 179 c First step: going from 1 to dissip_fac_mid (in gcm.def) 180 c============ 181 DO l=1,llm 182 zz = 1. - preff/presnivs(l) 183 zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz ) 184 ENDDO 185 186 write(*,*) 'Dissipation : ' 187 write(*,*) 'Multiplication de la dissipation en altitude :', 188 write(*,*) ' dissip_fac_mid =', dissip_fac_mid 189 190 c Second step if ok_strato: from dissip_fac_mid to dissip_fac_up (in gcm.def) 191 c========================== 192 c Utilisation de la fonction tangente hyperbolique pour augmenter 193 c arbitrairement la dissipation et donc la stabilite du modele a 194 c partir d'une certaine altitude. 195 196 c Le facteur multiplicatif de basse atmosphere etant deja pris 197 c en compte, il faut diviser le facteur multiplicatif de haute 198 c atmosphere par celui-ci. 199 200 if (ok_strato) then 201 202 Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta) 203 do l=1,llm 204 zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1) 205 & *(1-(0.5*(1+tanh(-6./dissip_deltaz* 206 & (-dissip_hdelta*log(presnivs(l)/Pup)) )))) )) 207 enddo 208 209 write(*,*) ' dissip_fac_up =', dissip_fac_up 210 write(*,*) 'Transition mid /up: Pupstart,delta =', 211 & dissip_pupstart,'Pa', dissip_deltaz , '(km)' 212 195 213 endif 196 214 -
trunk/libf/dyn3dpar/leapfrog_p.F
r108 r109 271 271 272 272 CALL pression ( ip1jmp1, ap, bp, ps, p ) 273 CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) 273 if (planet_type.eq."earth") then 274 CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) 275 else 276 CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf ) 277 endif 274 278 c$OMP END MASTER 275 279 c----------------------------------------------------------------------- … … 746 750 747 751 c$OMP BARRIER 748 CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf ) 752 if (planet_type.eq."earth") then 753 CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) 754 else 755 CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf ) 756 endif 749 757 c$OMP BARRIER 750 758 jD_cur = jD_ref + day_ini - day_ref … … 1094 1102 CALL pression_p ( ip1jmp1, ap, bp, ps, p ) 1095 1103 c$OMP BARRIER 1096 CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) 1104 if (planet_type.eq."earth") then 1105 CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) 1106 else 1107 CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf ) 1108 endif 1097 1109 c$OMP BARRIER 1098 1110 -
trunk/libf/dyn3dpar/limy.F
r1 r109 40 40 REAL qbyv(ip1jm,llm) 41 41 42 REAL qpns,qpsn,ap n,aps,dyn1,dys1,dyn2,dys242 REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2 43 43 Logical extremum,first 44 44 save first … … 116 116 117 117 c print*,dyqv(iip1+1) 118 c ap n=abs(dyq(1)/dyqv(iip1+1))118 c appn=abs(dyq(1)/dyqv(iip1+1)) 119 119 c print*,dyq(ip1jm+1) 120 120 c print*,dyqv(ip1jm-iip1+1) 121 c ap s=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))121 c apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) 122 122 c do ij=2,iim 123 c ap n=amax1(abs(dyq(ij)/dyqv(ij)),apn)124 c ap s=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)123 c appn=amax1(abs(dyq(ij)/dyqv(ij)),appn) 124 c apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps) 125 125 c enddo 126 c ap n=min(pente_max/apn,1.)127 c ap s=min(pente_max/aps,1.)126 c appn=min(pente_max/appn,1.) 127 c apps=min(pente_max/apps,1.) 128 128 129 129 … … 131 131 132 132 c if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.) 133 c & ap n=0.133 c & appn=0. 134 134 c if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* 135 135 c & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.) 136 c & ap s=0.136 c & apps=0. 137 137 138 138 c limitation des pentes aux poles 139 139 c do ij=1,iip1 140 c dyq(ij)=ap n*dyq(ij)141 c dyq(ip1jm+ij)=ap s*dyq(ip1jm+ij)140 c dyq(ij)=appn*dyq(ij) 141 c dyq(ip1jm+ij)=apps*dyq(ip1jm+ij) 142 142 c enddo 143 143 -
trunk/libf/dyn3dpar/vlsplt_p.F
r1 r109 561 561 REAL qbyv(ip1jm,llm) 562 562 563 REAL qpns,qpsn,ap n,aps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs563 REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs 564 564 c REAL newq,oldmasse 565 565 Logical extremum,first,testcpu … … 732 732 C PRINT*,dyq(1) 733 733 C PRINT*,dyqv(iip1+1) 734 C ap n=abs(dyq(1)/dyqv(iip1+1))734 C appn=abs(dyq(1)/dyqv(iip1+1)) 735 735 C PRINT*,dyq(ip1jm+1) 736 736 C PRINT*,dyqv(ip1jm-iip1+1) 737 C ap s=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))737 C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) 738 738 C DO ij=2,iim 739 C ap n=amax1(abs(dyq(ij)/dyqv(ij)),apn)740 C ap s=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)739 C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn) 740 C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps) 741 741 C ENDDO 742 C ap n=min(pente_max/apn,1.)743 C ap s=min(pente_max/aps,1.)742 C appn=min(pente_max/appn,1.) 743 C apps=min(pente_max/apps,1.) 744 744 C 745 745 C … … 747 747 C 748 748 C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.) 749 C & ap n=0.749 C & appn=0. 750 750 C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* 751 751 C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.) 752 C & ap s=0.752 C & apps=0. 753 753 C 754 754 C limitation des pentes aux poles 755 755 C DO ij=1,iip1 756 C dyq(ij)=ap n*dyq(ij)757 C dyq(ip1jm+ij)=ap s*dyq(ip1jm+ij)756 C dyq(ij)=appn*dyq(ij) 757 C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij) 758 758 C ENDDO 759 759 C -
trunk/libf/dyn3dpar/vlspltqs_p.F
r8 r109 780 780 C PRINT*,dyq(1) 781 781 C PRINT*,dyqv(iip1+1) 782 C ap n=abs(dyq(1)/dyqv(iip1+1))782 C appn=abs(dyq(1)/dyqv(iip1+1)) 783 783 C PRINT*,dyq(ip1jm+1) 784 784 C PRINT*,dyqv(ip1jm-iip1+1) 785 C ap s=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))785 C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) 786 786 C DO ij=2,iim 787 C ap n=amax1(abs(dyq(ij)/dyqv(ij)),apn)788 C ap s=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)787 C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn) 788 C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps) 789 789 C ENDDO 790 C ap n=min(pente_max/apn,1.)791 C ap s=min(pente_max/aps,1.)790 C appn=min(pente_max/appn,1.) 791 C apps=min(pente_max/apps,1.) 792 792 C 793 793 C … … 795 795 C 796 796 C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.) 797 C & ap n=0.797 C & appn=0. 798 798 C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* 799 799 C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.) 800 C & ap s=0.800 C & apps=0. 801 801 C 802 802 C limitation des pentes aux poles 803 803 C DO ij=1,iip1 804 C dyq(ij)=ap n*dyq(ij)805 C dyq(ip1jm+ij)=ap s*dyq(ip1jm+ij)804 C dyq(ij)=appn*dyq(ij) 805 C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij) 806 806 C ENDDO 807 807 C
Note: See TracChangeset
for help on using the changeset viewer.