Changeset 108 for trunk/libf/dyn3d
- Timestamp:
- Apr 12, 2011, 11:16:02 AM (14 years ago)
- Location:
- trunk/libf/dyn3d
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/libf/dyn3d/calfis.F
r97 r108 587 587 588 588 ! ADAPTATION GCM POUR CP(T) 589 ztfi=ztfi+zdtfi*dtphys590 589 call t2tpot(ngridmx*llm,ztfi,zteta,zpk) 591 590 -
trunk/libf/dyn3d/comconst.h
r53 r108 6 6 7 7 COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl, & 8 & iflag_top_bound 8 & iflag_top_bound,mode_top_bound 9 9 COMMON/comconstr/dtvr,daysec, & 10 10 & pi,dtphys,dtdiss,rad,r,cpp,kappa,cotot,unsim,g,omeg & 11 & ,dissip_factz,dissip_deltaz,dissip_zref&12 & 11 & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta & 12 & ,dissip_pupstart ,tau_top_bound, & 13 13 & daylen,year_day,molmass, ihf 14 14 COMMON/cpdetvenus/nu_venus,t0_venus … … 28 28 REAL g ! (m/s2) gravity 29 29 REAL omeg ! (rad/s) rotation rate of the planet 30 REAL dissip_factz,dissip_deltaz,dissip_zref 31 INTEGER iflag_top_bound 30 REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta 31 REAL dissip_pupstart 32 INTEGER iflag_top_bound,mode_top_bound 32 33 REAL tau_top_bound 33 34 REAL daylen ! length of solar day, in 'standard' day length -
trunk/libf/dyn3d/conf_gcm.F
r97 r108 303 303 ! Parametres controlant la variation sur la verticale des constantes de 304 304 ! dissipation. 305 ! Pour le moment actifs uniquement dans la version a 39 niveaux 306 ! avec ok_strato=y 307 308 dissip_factz=4. 309 dissip_deltaz=10. 310 dissip_zref=30. 311 CALL getin('dissip_factz',dissip_factz ) 305 ! Actifs uniquement avec ok_strato=y 306 307 dissip_fac_mid=2. 308 dissip_fac_up=10. 309 dissip_deltaz=10.! Intervalle (km) pour le changement mid / up 310 dissip_hdelta=5. ! scale height (km) dans la zone de la transition(m) 311 dissip_pupstart=1.e3 ! pression (Pa) au bas la transition mid / up 312 CALL getin('dissip_fac_mid',dissip_fac_mid ) 313 CALL getin('dissip_fac_up',dissip_fac_up ) 312 314 CALL getin('dissip_deltaz',dissip_deltaz ) 313 CALL getin('dissip_zref',dissip_zref ) 314 315 CALL getin('dissip_hdelta',dissip_hdelta ) 316 CALL getin('dissip_pupstart',dissip_pupstart ) 317 318 ! Parametres controlant la couche eponge 319 ! Actifs uniquement avec ok_strato=y 320 ! mode = 0 : pas de sponge 321 ! mode = 1 : u et v -> 0 322 ! mode = 2 : u et v -> moyenne zonale 323 ! mode = 3 : u, v et h -> moyenne zonale 315 324 iflag_top_bound=1 325 mode_top_bound=3 316 326 tau_top_bound=1.e-5 317 327 CALL getin('iflag_top_bound',iflag_top_bound) 328 CALL getin('mode_top_bound',mode_top_bound) 318 329 CALL getin('tau_top_bound',tau_top_bound) 319 330 -
trunk/libf/dyn3d/inidissip.F
r1 r108 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/dyn3d/leapfrog.F
r101 r108 111 111 REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm) 112 112 REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1) 113 114 real :: duspg(ip1jmp1,llm) ! for bilan_dyn115 113 116 114 c variables pour le fichier histoire … … 463 461 c ------------------- 464 462 IF (ok_strato) THEN 465 CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi) 463 CALL top_bound( vcov,ucov,teta,masse,dutop,dvtop,dtetatop) 464 c dqtop=0, dptop=0 465 CALL addfi( dtphys, leapf, forward , 466 $ ucov, vcov, teta , q ,ps , 467 $ dutop, dvtop, dtetatop , dqtop ,dptop ) 466 468 ENDIF 467 c dqtop=0, dptop=0 468 469 c ajout des tendances physiques: 470 c ------------------------------ 471 CALL addfi( dtphys, leapf, forward , 472 $ ucov, vcov, teta , q ,ps , 473 $ dufi, dvfi, dtetafi , dqfi ,dpfi ) 474 c 469 475 470 c Diagnostique de conservation de l'énergie : difference 476 471 IF (ip_ebil_dyn.ge.1 ) THEN … … 512 507 ! Sponge layer (if any) 513 508 IF (ok_strato) THEN 514 dutop(:,:)=0.515 dvtop(:,:)=0.516 dtetatop(:,:)=0.517 dqtop(:,:,:)=0.518 dptop(:)=0.519 509 CALL top_bound(vcov,ucov,teta,masse,dutop,dvtop,dtetatop) 510 c dqtop=0, dptop=0 520 511 CALL addfi( dtvr, leapf, forward , 521 512 $ ucov, vcov, teta , q ,ps , … … 669 660 CALL bilan_dyn(dtvr*iperiod,dtvr*day_step*periodav, 670 661 & ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov, 671 & du,dudis,du spg,dufi)662 & du,dudis,dutop,dufi) 672 663 #endif 673 664 END IF … … 801 792 CALL bilan_dyn(dtvr*iperiod,dtvr*day_step*periodav, 802 793 & ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov, 803 & du,dudis,du spg,dufi)794 & du,dudis,dutop,dufi) 804 795 #endif 805 796 ENDIF -
trunk/libf/dyn3d/top_bound.F
r1 r108 52 52 LOGICAL,SAVE :: first=.true. 53 53 54 REAL zkm 54 55 INTEGER j,l 55 56 … … 68 69 rdamp(llm-3)=tau_top_bound/8. 69 70 else if (iflag_top_bound.eq.2) then 70 ! couc e eponge dans toutes les couches de pression plus faible que71 ! couche eponge dans toutes les couches de pression plus faible que 71 72 ! 100 fois la pression de la derniere couche 72 73 rdamp(:)=tau_top_bound … … 74 75 endif 75 76 first=.false. 76 print*,'TOP_BOUND rdamp=',rdamp 77 print*,'TOP_BOUND mode',mode_top_bound 78 print*,'Coeffs pour la couche eponge a l equateur' 79 print*,'p (Pa) z(km) tau (s) dt*rdamp' 80 do l=1,llm 81 if (rdamp(l).ne.0.) then 82 zkm = phi(iip1/2,jjp1/2,l)/(1000*g) 83 print*,presnivs(l),zkm, 84 . 1./rdamp(l), 85 . dt*rdamp(l) 86 endif 87 enddo 77 88 endif 78 89 79 90 CALL massbar(masse,massebx,masseby) 80 91 81 do l=1,llm 92 c mode = 0 : pas de sponge 93 c mode = 1 : u et v -> 0 94 c mode = 2 : u et v -> moyenne zonale 95 c mode = 3 : u, v et h -> moyenne zonale 96 97 if (mode_top_bound.ge.2) then 98 do l=1,llm 82 99 do j=1,jjm 83 100 vzon(j,l)=0. … … 92 109 vzon(j,l)=vzon(j,l)/zm 93 110 enddo 94 enddo111 enddo 95 112 96 do l=1,llm 97 do i=1,iip1 98 do j=1,jjm 99 dv(i,j,l)=dv(i,j,l)-rdamp(l)*(vcov(i,j,l)-vzon(j,l)) 100 enddo 101 enddo 102 enddo 103 104 do l=1,llm 113 do l=1,llm 105 114 do j=2,jjm 106 115 uzon(j,l)=0. … … 112 121 uzon(j,l)=uzon(j,l)/zm 113 122 enddo 114 enddo 123 enddo 124 else 125 vzon(:,:)=0. 126 uzon(:,:)=0. 127 endif 115 128 116 do l=1,llm 129 if (mode_top_bound.ge.3) then 130 do l=1,llm 117 131 do j=2,jjm 118 132 zm=0. … … 124 138 tzon(j,l)=tzon(j,l)/zm 125 139 enddo 126 enddo 140 enddo 141 endif 127 142 128 143 C AMORTISSEMENTS LINEAIRES: 129 144 130 do l=1,llm 145 if (mode_top_bound.ge.1) then 146 do l=1,llm 147 do i=1,iip1 148 do j=1,jjm 149 dv(i,j,l)= -rdamp(l)*(vcov(i,j,l)-vzon(j,l)) 150 enddo 151 enddo 152 enddo 153 154 do l=1,llm 131 155 do i=1,iip1 132 156 do j=2,jjm 133 du(i,j,l)=du(i,j,l) 134 s -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l)) 135 dh(i,j,l)=dh(i,j,l)-rdamp(l)*(teta(i,j,l)-tzon(j,l)) 157 du(i,j,l)= -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l)) 136 158 enddo 137 159 enddo 138 enddo 160 enddo 161 endif 162 163 if (mode_top_bound.ge.3) then 164 do l=1,llm 165 do i=1,iip1 166 do j=2,jjm 167 dh(i,j,l)= -rdamp(l)*(teta(i,j,l)-tzon(j,l)) 168 enddo 169 enddo 170 enddo 171 endif 139 172 140 141 173 RETURN 142 174 END
Note: See TracChangeset
for help on using the changeset viewer.