Changeset 1322 for LMDZ4/branches/LMDZ4V5.0-dev/libf
- Timestamp:
- Mar 12, 2010, 11:54:11 AM (15 years ago)
- Location:
- LMDZ4/branches/LMDZ4V5.0-dev/libf
- Files:
-
- 14 edited
-
dyn3d/etat0_netcdf.F (modified) (2 diffs)
-
dyn3dpar/etat0_netcdf.F (modified) (2 diffs)
-
dyn3dpar/leapfrog_p.F (modified) (1 diff)
-
phylmd/calwake.F (modified) (6 diffs)
-
phylmd/conf_phys.F90 (modified) (5 diffs)
-
phylmd/cv3_cine.F (modified) (3 diffs)
-
phylmd/cv3_routines.F (modified) (8 diffs)
-
phylmd/cv3p1_closure.F (modified) (4 diffs)
-
phylmd/ini_wake.F (modified) (5 diffs)
-
phylmd/phyetat0.F (modified) (4 diffs)
-
phylmd/phyredem.F (modified) (2 diffs)
-
phylmd/phys_state_var_mod.F90 (modified) (5 diffs)
-
phylmd/physiq.F (modified) (13 diffs)
-
phylmd/wake.F (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3d/etat0_netcdf.F
r1299 r1322 144 144 REAL :: solarlong0 145 145 real :: seuil_inversion 146 real :: alp_offset 146 147 147 148 integer read_climoz ! read ozone climatology … … 184 185 & iflag_thermals,nsplit_thermals,tau_thermals, & 185 186 & iflag_thermals_ed,iflag_thermals_optflux, & 186 & iflag_coupl,iflag_clos,iflag_wake, read_climoz ) 187 & iflag_coupl,iflag_clos,iflag_wake, read_climoz, & 188 & alp_offset) 187 189 188 190 ! co2_ppm0 : initial value of atmospheric CO2 from .def file (co2_ppm value) -
LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/etat0_netcdf.F
r1299 r1322 144 144 REAL :: solarlong0 145 145 real :: seuil_inversion 146 REAL :: alp_offset 146 147 147 148 integer read_climoz ! read ozone climatology … … 184 185 & iflag_thermals,nsplit_thermals,tau_thermals, & 185 186 & iflag_thermals_ed,iflag_thermals_optflux, & 186 & iflag_coupl,iflag_clos,iflag_wake, read_climoz ) 187 & iflag_coupl,iflag_clos,iflag_wake, read_climoz, & 188 & alp_offset) 187 189 188 190 ! co2_ppm0 : initial value of atmospheric CO2 from .def file (co2_ppm value) -
LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/leapfrog_p.F
r1299 r1322 677 677 call suspend_timer(timer_caldyn) 678 678 679 if (prt_level >= 10) then 679 680 write(lunout,*) 680 681 & 'leapfrog_p: Entree dans la physique : Iteration No ',true_itau 682 endif 681 683 c$OMP END MASTER 682 684 -
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/calwake.F
r990 r1322 1 ! 2 ! $Id$ 3 ! 1 4 SUBROUTINE CALWAKE(paprs,pplay,dtime 2 5 : ,t,q,omgb … … 189 192 $ ,Cstar,d_deltat_gw 190 193 $ ,d_deltatw,d_deltaqw) 191 192 DO i=1,klon193 IF (ktopw(i) .GT. 0) THEN194 DO l=1,klev194 c 195 DO l=1,klev 196 DO i=1,klon 197 IF (ktopw(i) .GT. 0) THEN 195 198 wake_deltat(i,l)= dtw(i,l) 196 199 wake_deltaq(i,l)= dqw(i,l) … … 212 215 wake_ddeltat(i,l) = d_deltatw(i,l) 213 216 wake_ddeltaq(i,l) = d_deltaqw(i,l) 214 ENDDO 215 ELSE 216 DO l = 1,klev 217 ELSE 217 218 wake_deltat(i,l)= 0. 218 219 wake_deltaq(i,l)= 0. … … 222 223 wake_dtKE(i,l) = 0. 223 224 wake_dqKE(i,l) = 0. 225 wake_dtPBL(i,l) = 0. 226 wake_dqPBL(i,l) = 0. 224 227 wake_omg(i,l) = 0. 225 228 wake_dp_deltomg(i,l) = 0. … … 230 233 undi_t(i,l)=te(i,l) 231 234 undi_q(i,l)=qe(i,l) 232 ENDDO 233 ENDIF 234 235 wake_ddeltat(i,l) = 0. 236 wake_ddeltaq(i,l) = 0. 237 ENDIF 238 ENDDO 239 ENDDO 240 c 241 DO i=1,klon 235 242 wake_h(i)= hw(i) 236 243 wake_s(i)= sigmaw(i) … … 241 248 wake_Cstar(i) = Cstar(i) 242 249 wake_dens(i) = wdens(i) 243 c 244 cIM 290108 999 CONTINUE 245 c 246 ENDDO 250 ENDDO 251 c 247 252 RETURN 248 253 END 254 249 255 SUBROUTINE CALWAKE_scal(paprs,pplay,dtime 250 256 : ,t,q,omgb -
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/conf_phys.F90
r1299 r1322 22 22 iflag_thermals,nsplit_thermals,tau_thermals, & 23 23 iflag_thermals_ed,iflag_thermals_optflux, & 24 iflag_coupl,iflag_clos,iflag_wake, read_climoz) 24 iflag_coupl,iflag_clos,iflag_wake, read_climoz, & 25 alp_offset) 25 26 26 27 use IOIPSL … … 110 111 integer :: iflag_clos 111 112 integer :: iflag_wake 113 real :: alp_offset 114 REAL, SAVE :: alp_offset_omp 112 115 integer,SAVE :: iflag_coupl_omp,iflag_clos_omp,iflag_wake_omp 113 116 integer,SAVE :: iflag_cvl_sigd_omp … … 1032 1035 iflag_wake_omp = 0 1033 1036 call getin('iflag_wake',iflag_wake_omp) 1037 1038 !Config Key = alp_offset 1039 !Config Desc = 1040 !Config Def = 0 1041 !Config Help = 1042 ! 1043 alp_offset_omp = 0.01 1044 call getin('alp_offset',alp_offset_omp) 1034 1045 1035 1046 ! … … 1449 1460 iflag_clos = iflag_clos_omp 1450 1461 iflag_wake = iflag_wake_omp 1462 alp_offset = alp_offset_omp 1451 1463 iflag_cvl_sigd = iflag_cvl_sigd_omp 1452 1464 type_run = type_run_omp … … 1616 1628 write(numout,*)' Fmax = ', Fmax 1617 1629 write(numout,*)' alphas = ', alphas 1630 write(numout,*)' iflag_wake = ', iflag_wake 1631 write(numout,*)' alp_offset = ', alp_offset 1618 1632 1619 1633 write(numout,*)' lonmin lonmax latmin latmax bilKP_ins =',& -
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/cv3_cine.F
r1146 r1322 1 ! 2 ! $Id$ 3 ! 1 4 SUBROUTINE cv3_cine(nloc,ncum,nd,icb,inb 2 5 : ,pbase,plcl,p,ph,tv,tvp 3 : ,cina,cinb )6 : ,cina,cinb,plfc) 4 7 5 8 *************************************************************** … … 26 29 c 27 30 c output 28 real cina(nloc),cinb(nloc) 31 real cina(nloc),cinb(nloc),plfc(nloc) 29 32 c 30 33 c local variables … … 34 37 logical lswitch(nloc),lswitch1(nloc),lswitch2(nloc) 35 38 logical exist_lfc(nloc) 36 real plfc(nloc)37 39 real dpmax 38 40 real deltap,dcin -
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/cv3_routines.F
r1299 r1322 74 74 c dtcrit = -5.0 75 75 c tau = 3000. 76 cc tau = 1800. 77 c tau= 2800. 78 tau=8000. 76 tau = 1800. 77 cc tau=8000. 79 78 beta = 1.0 - delt/tau 80 79 alpha1 = 1.5e-3 … … 2094 2093 cc---end jyg--- 2095 2094 c 2096 c--------retour à la formulation originale d' 'Emanuel.2095 c--------retour à la formulation originale d'Emanuel. 2097 2096 b6=bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac 2098 2097 c6=water(il,i+1)+bfac*wdtrain(il) … … 2100 2099 if(c6.gt.0.0)then 2101 2100 revap=0.5*(-b6+sqrt(b6*b6+4.*c6)) 2102 water(il,i)=revap*revap !equation de conservation 2101 cjyg Dans sa formulation originale, Emanuel calcule l'evaporation par: 2102 cc evap(il,i)=sigt*afac*revap 2103 c ce qui n'est pas correct. Dans cv_routines, la formulation a été modifiee. 2104 c Ici,l'evaporation evap est simplement calculee par l'equation de 2105 c conservation. 2106 water(il,i)=revap*revap 2103 2107 else 2108 cjyg---- Correction : si c6 <= 0, water(il,i)=0. 2104 2109 water(il,i) = 0. 2105 2110 endif … … 2345 2350 real esum(nloc), fsum(nloc), gsum(nloc), hsum(nloc) 2346 2351 real th_wake(nloc,nd) 2347 real alpha_qpos(nloc) 2352 real alpha_qpos(nloc),alpha_qpos1(nloc) 2348 2353 real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd) ! cld 2349 2354 real siga(nloc,nd), sax(nloc,nd), mac(nloc,nd) ! cld … … 3050 3055 do il=1,ncum 3051 3056 IF (iflag(il) .le. 1) THEN 3057 IF (cvflag_grav) then 3058 ex=0.01*grav*ment(il,inb(il),inb(il)) 3059 : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j)) 3060 : /(ph(i l,inb(il))-ph(il,inb(il)+1)) 3061 ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex 3062 ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j) 3063 : +ex*(ph(il,inb(il))-ph(il,inb(il)+1)) 3064 : /(ph(il,inb(il)-1)-ph(il,inb(il))) 3065 else 3052 3066 ex=0.1*ment(il,inb(il),inb(il)) 3053 3067 : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j)) … … 3057 3071 : +ex*(ph(il,inb(il))-ph(il,inb(il)+1)) 3058 3072 : /(ph(il,inb(il)-1)-ph(il,inb(il))) 3073 ENDIF !cvflag grav 3059 3074 ENDIF !iflag 3060 3075 enddo … … 3122 3137 c in order to ensure moisture positivity 3123 3138 DO il = 1,ncum 3139 alpha_qpos(il)=1. 3124 3140 IF (iflag(il) .le. 1) THEN 3125 alpha_qpos(il) = max(1. , -delt*fr(il,1)/ 3141 if (fr(il,1) .le. 0.) then 3142 alpha_qpos(il) = max(alpha_qpos(il) , 3143 : (-delt*fr(il,1))/ 3126 3144 : (s_wake(il)*rr_wake(il,1)+(1.-s_wake(il))*rr(il,1))) 3145 end if 3127 3146 ENDIF 3128 3147 ENDDO … … 3130 3149 DO il = 1,ncum 3131 3150 IF (iflag(il) .le. 1) THEN 3132 alpha_qpos(il) = max(alpha_qpos(il) , -delt*fr(il,i)/ 3151 IF (fr(il,i) .le. 0.) THEN 3152 alpha_qpos1(il)=max(1. , (-delt*fr(il,i))/ 3133 3153 : (s_wake(il)*rr_wake(il,i)+(1.-s_wake(il))*rr(il,i))) 3154 IF (alpha_qpos1(il) .ge. alpha_qpos(il)) 3155 : alpha_qpos(il)=alpha_qpos1(il) 3156 ENDIF 3134 3157 ENDIF 3135 3158 ENDDO -
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/cv3p1_closure.F
r1299 r1322 63 63 integer nsupmax(nloc) 64 64 real supcrit,temp(nloc,nd) 65 real P1(nloc),Pmin(nloc) 65 real P1(nloc),Pmin(nloc),plfc(nloc) 66 66 real asupmax0(nloc) 67 67 logical ok(nloc) … … 385 385 CALL cv3_cine (nloc,ncum,nd,icb,inb 386 386 : ,pbase,plcl,p,ph,tv,tvp 387 : ,cina,cinb )387 : ,cina,cinb,plfc) 388 388 c 389 389 DO il = 1,ncum … … 495 495 do k= 1,nl 496 496 do il = 1,ncum 497 !IM IF (k .ge. icb(il) .and. k .le. inb(il)) THEN 498 IF (k .ge. icb(il)+1 .and. k .le. inb(il)) THEN 497 !old IF (k .ge. icb(il) .and. k .le. inb(il)) THEN 498 !IM IF (k .ge. icb(il)+1 .and. k .le. inb(il)) THEN 499 IF (k .ge. icb(il) .and. k .le. inb(il) !cor jyg 500 $ .and. icb(il)+1 .le. inb(il)) THEN !cor jyg 499 501 cbmflim(il) = cbmflim(il)+MLIM(il,k) 500 502 ENDIF … … 548 550 do il = 1,ncum 549 551 IF ( k .ge. icb(il)+1 .AND. k .le. inb(il)) THEN 550 sig(il,k) = beta*sig(il,k)+(1.-beta)*coef(il)*siglim(il,k) 551 cc sig(il,k) = beta*sig(il,k)+siglim(il,k) 552 w0(il,k) = beta*w0(il,k) +(1.-beta)*wlim(il,k) 553 AMU=SIG(il,k)*W0(il,k) 552 amu=beta*sig(il,k)*w0(il,k)+ 553 : (1.-beta)*coef(il)*siglim(il,k)*wlim(il,k) 554 w0(il,k) = wlim(il,k) 555 w0(il,k) =max(w0(il,k),1.e-10) 556 sig(il,k)=amu/w0(il,k) 557 sig(il,k)=min(sig(il,k),1.) 554 558 cc amu = 0.5*(SIG(il,k)+sigold(il,k))*W0(il,k) 555 559 M(il,k)=AMU*0.007*P(il,k)*(PH(il,k)-PH(il,k+1))/TV(il,k) -
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/ini_wake.F
r970 r1322 1 ! 2 ! $Id$ 3 ! 1 4 SUBROUTINE INI_WAKE(wape,fip,it_wape_prescr, 2 5 : wape_prescr, fip_prescr, alp_bl_prescr, ale_bl_prescr) … … 23 26 c wape_prescr : valeur prescrite de la WAPE. 24 27 c fip_prescr : valeur prescrite de la FIP. 28 c ale_bl_prescr : valeur prescrite de la Ale de PBL. 29 c alp_bl_prescr : valeur prescrite de la Alp de PBL. 25 30 c 26 31 c Variables internes … … 29 34 c w = WAPE lue 30 35 c f = FIP lue 36 c alebl = Ale de PBL lue 37 c alpbl = Alp de PBL lue 31 38 c 39 include 'iniprint.h' 32 40 cdeclarations 33 41 real ale_bl_prescr 34 42 real alp_bl_prescr 35 43 real it 36 cCR: on rajoute ale et alp de la PBL precrits37 c open (99,file='wake.data',form='formatted')38 c read (99,*) it39 c read (99,*) w40 c read (99,*) f41 c read (99,*) u42 c read (99,*) p43 c close (99)44 44 45 45 ! FH A mettre si besoin dans physiq.def … … 48 48 w=4. 49 49 f=0.1 50 u=0.151 p=4.50 alebl=4. 51 alpbl=0.1 52 52 c 53 print *,' it,w ',it,w 53 cCR: on rajoute ale et alp de la PBL precrits 54 open (99,file='wake.data',form='formatted',status='old',err=902) 55 read (99,*) it 56 read (99,*) w 57 read (99,*) f 58 read (99,*,end=901) alebl 59 read (99,*,end=901) alpbl 60 901 close (99) 61 902 continue 62 c 63 write(lunout,*)' it,wape ',it,wape 54 64 it_wape_prescr = it 55 65 if (w .lt. 0) then … … 61 71 endif 62 72 c 63 print *,' u,p ',u,p64 al p_bl_prescr=u65 al e_bl_prescr=p73 write(lunout,*)' alebl, alpbl ',alebl,alpbl 74 ale_bl_prescr=alebl 75 alp_bl_prescr=alpbl 66 76 print *,'Initialisation de la poche : WAPE, FIP imposees =' 67 77 $ ,wape_prescr, fip_prescr -
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phyetat0.F
r1311 r1322 1 1 ! 2 ! $ Header$2 ! $Id$ 3 3 ! 4 4 c … … 228 228 $ 'coherente ', i, zmasq(i), pctsrf(i, is_ter) 229 229 $ ,pctsrf(i, is_lic) 230 WRITE(*,*) 'Je force la coherence zmasq=fractint' 230 231 zmasq(i) = fractint(i) 231 232 ENDIF … … 238 239 $ 'coherente ', i, zmasq(i) , pctsrf(i, is_oce) 239 240 $ ,pctsrf(i, is_sic) 241 WRITE(*,*) 'Je force la coherence zmasq=fractint' 240 242 zmasq(i) = fractint(i) 241 243 ENDIF … … 987 989 PRINT*,'(ecart-type) wake_cstar:', xmin, xmax 988 990 c 991 c wake_pe 992 c 993 CALL get_field("WAKE_PE",wake_pe,found) 994 IF (.NOT. found) THEN 995 PRINT*, "phyetat0: Le champ <WAKE_PE> est absent" 996 PRINT*, "Depart legerement fausse. Mais je continue" 997 wake_pe=0. 998 ENDIF 999 xmin = 1.0E+20 1000 xmax = -1.0E+20 1001 xmin = MINval(wake_pe) 1002 xmax = MAXval(wake_pe) 1003 PRINT*,'(ecart-type) wake_pe:', xmin, xmax 1004 c 989 1005 c wake_fip 990 1006 c -
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phyredem.F
r1302 r1322 299 299 !!!!!!!!!!!!!!!!!!!! FIN TKE PBL !!!!!!!!!!!!!!!!!!!!!!!!! 300 300 cIM ajout zmax0, f0, ema_work1, ema_work2 301 cIM wake_deltat, wake_deltaq, wake_s, wake_cstar, wake_ fip301 cIM wake_deltat, wake_deltaq, wake_s, wake_cstar, wake_pe, wake_fip 302 302 303 303 CALL put_field("ZMAX0","ZMAX0",zmax0) … … 318 318 CALL put_field("WAKE_CSTAR","WAKE_CSTAR",wake_cstar) 319 319 320 CALL put_field("WAKE_PE","WAKE_PE",wake_pe) 321 320 322 CALL put_field("WAKE_FIP","WAKE_FIP",wake_fip) 321 323 -
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phys_state_var_mod.F90
r1299 r1322 1 ! 2 ! $Id$ 3 ! 1 4 MODULE phys_state_var_mod 2 5 ! Variables sauvegardees pour le startphy.nc … … 169 172 ! wake_Cstar : vitesse d'etalement de la poche 170 173 ! wake_s : fraction surfacique occupee par la poche froide 174 ! wake_pe : wake potential energy - WAPE 171 175 ! wake_fip : Gust Front Impinging power - ALP 172 176 ! dt_wake, dq_wake: LS tendencies due to wake … … 179 183 REAL,ALLOCATABLE,SAVE :: wake_s(:) 180 184 !$OMP THREADPRIVATE(wake_s) 185 REAL,ALLOCATABLE,SAVE :: wake_pe(:) 186 !$OMP THREADPRIVATE(wake_pe) 181 187 REAL,ALLOCATABLE,SAVE :: wake_fip(:) 182 188 !$OMP THREADPRIVATE(wake_fip) … … 370 376 ALLOCATE(wght_th(klon,klev)) 371 377 ALLOCATE(wake_deltat(klon,klev), wake_deltaq(klon,klev)) 372 ALLOCATE(wake_Cstar(klon), wake_s(klon), wake_fip(klon)) 378 ALLOCATE(wake_Cstar(klon), wake_s(klon)) 379 ALLOCATE(wake_pe(klon), wake_fip(klon)) 373 380 ALLOCATE(dt_wake(klon,klev), dq_wake(klon,klev)) 374 381 ALLOCATE(pfrac_impa(klon,klev), pfrac_nucl(klon,klev)) … … 464 471 deallocate(lalim_conv, wght_th) 465 472 deallocate(wake_deltat, wake_deltaq) 466 deallocate(wake_Cstar, wake_s, wake_ fip)473 deallocate(wake_Cstar, wake_s, wake_pe, wake_fip) 467 474 deallocate(dt_wake, dq_wake) 468 475 deallocate(pfrac_impa, pfrac_nucl) -
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/physiq.F
r1311 r1322 625 625 REAL q_undi(klon,klev) ! humidite moyenne dans la zone non perturbee 626 626 c 627 REAL wake_pe(klon) ! Wake potential energy - WAPE 627 cjyg 628 ccc REAL wake_pe(klon) ! Wake potential energy - WAPE 628 629 629 630 REAL wake_gfl(klon) ! Gust Front Length … … 641 642 REAL dt_a(klon,klev) 642 643 REAL dq_a(klon,klev) 644 REAL, SAVE :: alp_offset 645 c$OMP THREADPRIVATE(alp_offset) 646 643 647 c 644 648 cRR:fin declarations poches froides … … 1188 1192 . iflag_thermals_ed,iflag_thermals_optflux, 1189 1193 c nv flags pour la convection et les poches froides 1190 . iflag_coupl,iflag_clos,iflag_wake, read_climoz) 1194 . iflag_coupl,iflag_clos,iflag_wake, read_climoz, 1195 & alp_offset) 1191 1196 call phys_state_var_init(read_climoz) 1192 1197 print*, '=================================================' … … 2021 2026 c 2022 2027 ccalcul de ale_wake et alp_wake 2023 do i = 1,klon 2024 if (iflag_wake.eq.1) then 2025 ale_wake(i) = 0.5*wake_cstar(i)**2 2026 alp_wake(i) = wake_fip(i) 2027 else 2028 ale_wake(i) = 0. 2029 alp_wake(i) = 0. 2030 endif 2031 enddo 2028 if (iflag_wake.eq.1) then 2029 if (itap .le. it_wape_prescr) then 2030 do i = 1,klon 2031 ale_wake(i) = wape_prescr 2032 alp_wake(i) = fip_prescr 2033 enddo 2034 else 2035 do i = 1,klon 2036 cjyg ALE=WAPE au lieu de ALE = 1/2 Cstar**2 2037 ccc ale_wake(i) = 0.5*wake_cstar(i)**2 2038 ale_wake(i) = wake_pe(i) 2039 alp_wake(i) = wake_fip(i) 2040 enddo 2041 endif 2042 else 2043 do i = 1,klon 2044 ale_wake(i) = 0. 2045 alp_wake(i) = 0. 2046 enddo 2047 endif 2032 2048 ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees 2033 2049 cdans le thermique sinon 2034 2050 if (iflag_coupl.eq.0) then 2035 if (debut) print*,'ALE et ALP imposes' 2051 if (debut.and.prt_level.gt.9) 2052 $ WRITE(lunout,*)'ALE et ALP imposes' 2036 2053 do i = 1,klon 2037 2054 con ne couple que ale … … 2046 2063 do i = 1,klon 2047 2064 ALE(i) = max(ale_wake(i),Ale_bl(i)) 2048 ALP(i) = alp_wake(i) + Alp_bl(i) 2065 c avant ALP(i) = alp_wake(i) + Alp_bl(i) 2066 ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb 2049 2067 c write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i) 2050 2068 c write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i) … … 2289 2307 2290 2308 endif 2309 c 2310 c=================================================================== 2311 cJYG 2312 IF (ip_ebil_phy.ge.2) THEN 2313 ztit='after wake' 2314 CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime 2315 e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay 2316 s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) 2317 call diagphy(airephy,ztit,ip_ebil_phy 2318 e , zero_v, zero_v, zero_v, zero_v, zero_v 2319 e , zero_v, zero_v, zero_v, ztsol 2320 e , d_h_vcol, d_qt, d_ec 2321 s , fs_bound, fq_bound ) 2322 END IF 2323 2291 2324 c print*,'apres callwake iflag_cldcon=', iflag_cldcon 2292 2325 c … … 2402 2435 e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay 2403 2436 s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) 2437 call diagphy(airephy,ztit,ip_ebil_phy 2438 e , zero_v, zero_v, zero_v, zero_v, zero_v 2439 e , zero_v, zero_v, zero_v, ztsol 2440 e , d_h_vcol, d_qt, d_ec 2441 s , fs_bound, fq_bound ) 2404 2442 END IF 2405 2443 … … 2450 2488 enddo 2451 2489 tau_overturning_th(:)=zmax_th(:)/max(0.5*wmax_th(:),0.1) 2452 print*,'TAU TH OK ',tau_overturning_th(1),detr_therm(1,3) 2490 if(prt_level.ge.9) 2491 & write(lunout,*)'TAU TH OK ', 2492 & tau_overturning_th(1),detr_therm(1,3) 2453 2493 2454 2494 c On impose que l'air autour de la fraction couverte par le thermique … … 2561 2601 2562 2602 ! les ratqs sont une combinaison de ratqss et ratqsc 2563 ! write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs 2603 if(prt_level.ge.9) 2604 $ write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs 2564 2605 2565 2606 if (tau_ratqs>1.e-10) then … … 2785 2826 e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay 2786 2827 s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) 2828 call diagphy(airephy,ztit,ip_ebil_phy 2829 e , zero_v, zero_v, zero_v, zero_v, zero_v 2830 e , zero_v, zero_v, zero_v, ztsol 2831 e , d_h_vcol, d_qt, d_ec 2832 s , fs_bound, fq_bound ) 2787 2833 END IF 2788 2834 c … … 2992 3038 itaprad = itaprad + 1 2993 3039 2994 if (iflag_radia.eq.0 ) then3040 if (iflag_radia.eq.0 .and. prt_level.ge.9) then 2995 3041 print *,'--------------------------------------------------' 2996 3042 print *,'>>>> ATTENTION rayonnement desactive pour ce cas' … … 3195 3241 e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay 3196 3242 s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) 3243 call diagphy(airephy,ztit,ip_ebil_phy 3244 e , zero_v, zero_v, zero_v, zero_v, zero_v 3245 e , zero_v, zero_v, zero_v, ztsol 3246 e , d_h_vcol, d_qt, d_ec 3247 s , fs_bound, fq_bound ) 3197 3248 END IF 3198 3249 c … … 3262 3313 IF (offline) THEN 3263 3314 3264 print*,'Attention on met a 0 les thermiques pour phystoke' 3265 call phystokenc ( 3315 IF (prt_level.ge.9) 3316 $ print*,'Attention on met a 0 les thermiques pour phystoke' 3317 call phystokenc ( 3266 3318 I nlon,klev,pdtphys,rlon,rlat, 3267 3319 I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, -
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/wake.F
r1277 r1322 1 Subroutine WAKE (p,ph,ppi,dtime,sigd_con 1 ! 2 ! $Id$ 3 ! 4 Subroutine WAKE (p,ph,pi,dtime,sigd_con 2 5 : ,te0,qe0,omgb 3 6 : ,dtdwn,dqdwn,amdwn,amup,dta,dqa … … 24 27 c 25 28 use dimphy 29 IMPLICIT none 30 c============================================================================ 31 C 32 C 33 C But : Decrire le comportement des poches froides apparaissant dans les 34 C grands systemes convectifs, et fournir l'energie disponible pour 35 C le declenchement de nouvelles colonnes convectives. 36 C 37 C Variables d'etat : deltatw : ecart de temperature wake-undisturbed area 38 C deltaqw : ecart d'humidite wake-undisturbed area 39 C sigmaw : fraction d'aire occupee par la poche. 40 C 41 C Variable de sortie : 42 c 43 c wape : WAke Potential Energy 44 c fip : Front Incident Power (W/m2) - ALP 45 c gfl : Gust Front Length per unit area (m-1) 46 C dtls : large scale temperature tendency due to wake 47 C dqls : large scale humidity tendency due to wake 48 C hw : hauteur de la poche 49 C dp_omgb : vertical gradient of large scale omega 50 C wdens : densite de poches 51 C omgbdth: flux of Delta_Theta transported by LS omega 52 C dtKE : differential heating (wake - unpertubed) 53 C dqKE : differential moistening (wake - unpertubed) 54 C omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s) 55 C dp_deltomg : vertical gradient of omg (s-1) 56 C spread : spreading term in dt_wake and dq_wake 57 C deltatw : updated temperature difference (T_w-T_u). 58 C deltaqw : updated humidity difference (q_w-q_u). 59 C sigmaw : updated wake fractional area. 60 C d_deltat_gw : delta T tendency due to GW 61 c 62 C Variables d'entree : 63 c 64 c aire : aire de la maille 65 c te0 : temperature dans l'environnement (K) 66 C qe0 : humidite dans l'environnement (kg/kg) 67 C omgb : vitesse verticale moyenne sur la maille (Pa/s) 68 C dtdwn: source de chaleur due aux descentes (K/s) 69 C dqdwn: source d'humidite due aux descentes (kg/kg/s) 70 C dta : source de chaleur due courants satures et detrain (K/s) 71 C dqa : source d'humidite due aux courants satures et detra (kg/kg/s) 72 C amdwn: flux de masse total des descentes, par unite de 73 C surface de la maille (kg/m2/s) 74 C amup : flux de masse total des ascendances, par unite de 75 C surface de la maille (kg/m2/s) 76 C p : pressions aux milieux des couches (Pa) 77 C ph : pressions aux interfaces (Pa) 78 C pi : (p/p_0)**kapa (adim) 79 C dtime: increment temporel (s) 80 c 81 C Variables internes : 82 c 83 c rhow : masse volumique de la poche froide 84 C rho : environment density at P levels 85 C rhoh : environment density at Ph levels 86 C te : environment temperature | may change within 87 C qe : environment humidity | sub-time-stepping 88 C the : environment potential temperature 89 C thu : potential temperature in undisturbed area 90 C tu : temperature in undisturbed area 91 C qu : humidity in undisturbed area 92 C dp_omgb: vertical gradient og LS omega 93 C omgbw : wake average vertical omega 94 C dp_omgbw: vertical gradient of omgbw 95 C omgbdq : flux of Delta_q transported by LS omega 96 C dth : potential temperature diff. wake-undist. 97 C th1 : first pot. temp. for vertical advection (=thu) 98 C th2 : second pot. temp. for vertical advection (=thw) 99 C q1 : first humidity for vertical advection 100 C q2 : second humidity for vertical advection 101 C d_deltatw : terme de redistribution pour deltatw 102 C d_deltaqw : terme de redistribution pour deltaqw 103 C deltatw0 : deltatw initial 104 C deltaqw0 : deltaqw initial 105 C hw0 : hw initial 106 C sigmaw0: sigmaw initial 107 C amflux : horizontal mass flux through wake boundary 108 C wdens_ref: initial number of wakes per unit area (3D) or per 109 C unit length (2D), at the beginning of each time step 110 C Tgw : 1 sur la période de onde de gravité 111 c Cgw : vitesse de propagation de onde de gravité 112 c LL : distance entre 2 poches 113 114 c------------------------------------------------------------------------- 115 c Déclaration de variables 116 c------------------------------------------------------------------------- 117 118 #include "dimensions.h" 119 #include "YOMCST.h" 120 #include "cvthermo.h" 121 #include "iniprint.h" 122 123 c Arguments en entree 124 c-------------------- 125 126 REAL, dimension(klon,klev) :: p, pi 127 REAL, dimension(klon,klev+1) :: ph, omgb 128 REAL dtime 129 REAL, dimension(klon,klev) :: te0,qe0 130 REAL, dimension(klon,klev) :: dtdwn, dqdwn 131 REAL, dimension(klon,klev) :: wdtPBL,wdqPBL 132 REAL, dimension(klon,klev) :: udtPBL,udqPBL 133 REAL, dimension(klon,klev) :: amdwn, amup 134 REAL, dimension(klon,klev) :: dta, dqa 135 REAL, dimension(klon) :: sigd_con 136 137 c Sorties 138 c-------- 139 140 REAL, dimension(klon,klev) :: deltatw, deltaqw, dth 141 REAL, dimension(klon,klev) :: tu, qu 142 REAL, dimension(klon,klev) :: dtls, dqls 143 REAL, dimension(klon,klev) :: dtKE, dqKE 144 REAL, dimension(klon,klev) :: dtPBL, dqPBL 145 REAL, dimension(klon,klev) :: spread 146 REAL, dimension(klon,klev) :: d_deltatgw 147 REAL, dimension(klon,klev) :: d_deltatw2, d_deltaqw2 148 REAL, dimension(klon,klev+1) :: omgbdth, omg 149 REAL, dimension(klon,klev) :: dp_omgb, dp_deltomg 150 REAL, dimension(klon,klev) :: d_deltat_gw 151 REAL, dimension(klon) :: hw, sigmaw, wape, fip, gfl, Cstar 152 REAL, dimension(klon) :: wdens 153 INTEGER, dimension(klon) :: ktopw 154 155 c Variables internes 156 c------------------- 157 158 c Variables à fixer 159 REAL ALON 160 REAL coefgw 161 REAL :: wdens_ref 162 REAL stark 163 REAL alpk 164 REAL delta_t_min 165 INTEGER nsub 166 REAL dtimesub 167 REAL sigmad, hwmin,wapecut 168 REAL :: sigmaw_max 169 REAL :: dens_rate 170 REAL wdens0 171 cIM 080208 172 LOGICAL, dimension(klon) :: gwake 173 174 c Variables de sauvegarde 175 REAL, dimension(klon,klev) :: deltatw0 176 REAL, dimension(klon,klev) :: deltaqw0 177 REAL, dimension(klon,klev) :: te, qe 178 REAL, dimension(klon) :: sigmaw0, sigmaw1 179 180 c Variables pour les GW 181 REAL, DIMENSION(klon) :: LL 182 REAL, dimension(klon,klev) :: N2 183 REAL, dimension(klon,klev) :: Cgw 184 REAL, dimension(klon,klev) :: Tgw 185 186 c Variables liées au calcul de hw 187 REAL, DIMENSION(klon) :: ptop_provis, ptop, ptop_new 188 REAL, DIMENSION(klon) :: sum_dth 189 REAL, DIMENSION(klon) :: dthmin 190 REAL, DIMENSION(klon) :: z, dz, hw0 191 INTEGER, DIMENSION(klon) :: ktop, kupper 192 193 c Sub-timestep tendencies and related variables 194 REAL d_deltatw(klon,klev),d_deltaqw(klon,klev) 195 REAL d_te(klon,klev),d_qe(klon,klev) 196 REAL d_sigmaw(klon),alpha(klon) 197 REAL q0_min(klon),q1_min(klon) 198 LOGICAL wk_adv(klon), OK_qx_qw(klon) 199 REAL epsilon 200 DATA epsilon/1.e-15/ 201 202 c Autres variables internes 203 INTEGER isubstep, k, i 204 205 REAL, DIMENSION(klon) :: sum_thu, sum_tu, sum_qu,sum_thvu 206 REAL, DIMENSION(klon) :: sum_dq, sum_rho 207 REAL, DIMENSION(klon) :: sum_dtdwn, sum_dqdwn 208 REAL, DIMENSION(klon) :: av_thu, av_tu, av_qu, av_thvu 209 REAL, DIMENSION(klon) :: av_dth, av_dq, av_rho 210 REAL, DIMENSION(klon) :: av_dtdwn, av_dqdwn 211 212 REAL, DIMENSION(klon,klev) :: rho, rhow 213 REAL, DIMENSION(klon,klev+1) :: rhoh 214 REAL, DIMENSION(klon,klev) :: rhow_moyen 215 REAL, DIMENSION(klon,klev) :: zh 216 REAL, DIMENSION(klon,klev+1) :: zhh 217 REAL, DIMENSION(klon,klev) :: epaisseur1, epaisseur2 218 219 REAL, DIMENSION(klon,klev) :: the, thu 220 221 ! REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw 222 223 REAL, DIMENSION(klon,klev+1) :: omgbw 224 REAL, DIMENSION(klon) :: pupper 225 REAL, DIMENSION(klon) :: omgtop 226 REAL, DIMENSION(klon,klev) :: dp_omgbw 227 REAL, DIMENSION(klon) :: ztop, dztop 228 REAL, DIMENSION(klon,klev) :: alpha_up 229 230 REAL, dimension(klon) :: RRe1, RRe2 231 REAL :: RRd1, RRd2 232 REAL, DIMENSION(klon,klev) :: Th1, Th2, q1, q2 233 REAL, DIMENSION(klon,klev) :: D_Th1, D_Th2, D_dth 234 REAL, DIMENSION(klon,klev) :: D_q1, D_q2, D_dq 235 REAL, DIMENSION(klon,klev) :: omgbdq 236 237 REAL, dimension(klon) :: ff, gg 238 REAL, dimension(klon) :: wape2, Cstar2, heff 239 240 REAL, DIMENSION(klon,klev) :: Crep 241 REAL Crep_upper, Crep_sol 242 243 REAL, DIMENSION(klon,klev) :: ppi 244 245 ccc nrlmd 246 real, dimension(klon) :: death_rate,nat_rate 247 real, dimension(klon,klev) :: entr 248 real, dimension(klon,klev) :: detr 249 250 C------------------------------------------------------------------------- 251 c Initialisations 252 c------------------------------------------------------------------------- 253 254 c print*, 'wake initialisations' 255 256 c Essais d'initialisation avec sigmaw = 0.02 et hw = 10. 257 c------------------------------------------------------------------------- 258 259 DATA wapecut,sigmad, hwmin /5.,.02,10./ 260 ccc nrlmd 261 DATA sigmaw_max /0.4/ 262 DATA dens_rate /0.1/ 263 ccc 264 C Longueur de maille (en m) 265 c------------------------------------------------------------------------- 266 267 c ALON = 3.e5 268 ALON = 1.e6 269 270 271 C Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei) 272 c 273 c coefgw : Coefficient pour les ondes de gravité 274 c stark : Coefficient k dans Cstar=k*sqrt(2*WAPE) 275 c wdens : Densité de poche froide par maille 276 c------------------------------------------------------------------------- 277 278 ccc nrlmd coefgw=10 279 c coefgw=1 280 c wdens0 = 1.0/(alon**2) 281 ccc nrlmd wdens = 1.0/(alon**2) 282 ccc nrlmd stark = 0.50 283 cCRtest 284 ccc nrlmd alpk=0.1 285 c alpk = 1.0 286 c alpk = 0.5 287 c alpk = 0.05 288 c 289 stark = 0.33 290 Alpk = 0.25 291 wdens_ref = 8.e-12 292 coefgw = 4. 293 Crep_upper=0.9 294 Crep_sol=1.0 295 296 ccc nrlmd Lecture du fichier wake_param.data 297 c$$$ OPEN(99,file='wake_param.data',status='old', 298 c$$$ $ form='formatted',err=9999) 299 c$$$ READ(99,*,end=9998) stark 300 c$$$ READ(99,*,end=9998) Alpk 301 c$$$ READ(99,*,end=9998) wdens_ref 302 c$$$ READ(99,*,end=9998) coefgw 303 c$$$9998 Continue 304 c$$$ CLOSE(99) 305 c$$$9999 Continue 306 c 307 c Initialisation de toutes des densites a wdens_ref. 308 c Les densites peuvent evoluer si les poches debordent 309 c (voir au tout debut de la boucle sur les substeps) 310 wdens = wdens_ref 311 c 312 c print*,'stark',stark 313 c print*,'alpk',alpk 314 c print*,'wdens',wdens 315 c print*,'coefgw',coefgw 316 ccc 317 C Minimum value for |T_wake - T_undist|. Used for wake top definition 318 c------------------------------------------------------------------------- 319 320 delta_t_min = 0.2 321 322 C 1. - Save initial values and initialize tendencies 323 C -------------------------------------------------- 324 325 DO k=1,klev 326 DO i=1, klon 327 ppi(i,k)=pi(i,k) 328 deltatw0(i,k) = deltatw(i,k) 329 deltaqw0(i,k)= deltaqw(i,k) 330 te(i,k) = te0(i,k) 331 qe(i,k) = qe0(i,k) 332 dtls(i,k) = 0. 333 dqls(i,k) = 0. 334 d_deltat_gw(i,k)=0. 335 d_te(i,k) = 0. 336 d_qe(i,k) = 0. 337 d_deltatw(i,k) = 0. 338 d_deltaqw(i,k) = 0. 339 !IM 060508 beg 340 d_deltatw2(i,k)=0. 341 d_deltaqw2(i,k)=0. 342 !IM 060508 end 343 ENDDO 344 ENDDO 345 c sigmaw1=sigmaw 346 c IF (sigd_con.GT.sigmaw1) THEN 347 c print*, 'sigmaw,sigd_con', sigmaw, sigd_con 348 c ENDIF 349 DO i=1, klon 350 cc sigmaw(i) = amax1(sigmaw(i),sigd_con(i)) 351 sigmaw(i) = amax1(sigmaw(i),sigmad) 352 sigmaw(i) = amin1(sigmaw(i),0.99) 353 sigmaw0(i) = sigmaw(i) 354 wape(i) = 0. 355 wape2(i) = 0. 356 d_sigmaw(i) = 0. 357 ktopw(i) = 0 358 ENDDO 359 C 360 C 361 C 2. - Prognostic part 362 C -------------------- 363 C 364 C 365 C 2.1 - Undisturbed area and Wake integrals 366 C --------------------------------------------------------- 367 368 DO i=1, klon 369 z(i) = 0. 370 ktop(i)=0 371 kupper(i) = 0 372 sum_thu(i) = 0. 373 sum_tu(i) = 0. 374 sum_qu(i) = 0. 375 sum_thvu(i) = 0. 376 sum_dth(i) = 0. 377 sum_dq(i) = 0. 378 sum_rho(i) = 0. 379 sum_dtdwn(i) = 0. 380 sum_dqdwn(i) = 0. 381 382 av_thu(i) = 0. 383 av_tu(i) =0. 384 av_qu(i) =0. 385 av_thvu(i) = 0. 386 av_dth(i) = 0. 387 av_dq(i) = 0. 388 av_rho(i) =0. 389 av_dtdwn(i) =0. 390 av_dqdwn(i) = 0. 391 ENDDO 392 c 393 c Distance between wakes 394 DO i = 1,klon 395 LL(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i)) 396 ENDDO 397 C Potential temperatures and humidity 398 c---------------------------------------------------------- 399 DO k =1,klev 400 DO i=1, klon 401 ! write(*,*)'wake 1',i,k,rd,te(i,k) 402 rho(i,k) = p(i,k)/(rd*te(i,k)) 403 ! write(*,*)'wake 2',rho(i,k) 404 IF(k .eq. 1) THEN 405 ! write(*,*)'wake 3',i,k,rd,te(i,k) 406 rhoh(i,k) = ph(i,k)/(rd*te(i,k)) 407 ! write(*,*)'wake 4',i,k,rd,te(i,k) 408 zhh(i,k)=0 409 ELSE 410 ! write(*,*)'wake 5',rd,(te(i,k)+te(i,k-1)) 411 rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1))) 412 ! write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1) 413 zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1) 414 ENDIF 415 ! write(*,*)'wake 7',ppi(i,k) 416 the(i,k) = te(i,k)/ppi(i,k) 417 thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k) 418 tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i) 419 qu(i,k) = qe(i,k) - deltaqw(i,k)*sigmaw(i) 420 ! write(*,*)'wake 8',(rd*(te(i,k)+deltatw(i,k))) 421 rhow(i,k) = p(i,k)/(rd*(te(i,k)+deltatw(i,k))) 422 dth(i,k) = deltatw(i,k)/ppi(i,k) 423 ENDDO 424 ENDDO 425 426 DO k = 1, klev-1 427 DO i=1, klon 428 IF(k.eq.1) THEN 429 N2(i,k)=0 430 ELSE 431 N2(i,k)=amax1(0.,-RG**2/the(i,k)*rho(i,k)*(the(i,k+1)- 432 $ the(i,k-1))/(p(i,k+1)-p(i,k-1))) 433 ENDIF 434 ZH(i,k)=(zhh(i,k)+zhh(i,k+1))/2 435 436 Cgw(i,k)=sqrt(N2(i,k))*ZH(i,k) 437 Tgw(i,k)=coefgw*Cgw(i,k)/LL(i) 438 ENDDO 439 ENDDO 440 441 DO i=1, klon 442 N2(i,klev)=0 443 ZH(i,klev)=0 444 Cgw(i,klev)=0 445 Tgw(i,klev)=0 446 ENDDO 447 448 c Calcul de la masse volumique moyenne de la colonne (bdlmd) 449 c----------------------------------------------------------------- 450 451 DO k=1,klev 452 DO i=1, klon 453 epaisseur1(i,k)=0. 454 epaisseur2(i,k)=0. 455 ENDDO 456 ENDDO 457 458 DO i=1, klon 459 epaisseur1(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1. 460 epaisseur2(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1. 461 rhow_moyen(i,1) = rhow(i,1) 462 ENDDO 463 464 DO k = 2, klev 465 DO i=1, klon 466 epaisseur1(i,k)= -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) +1. 467 epaisseur2(i,k)=epaisseur2(i,k-1)+epaisseur1(i,k) 468 rhow_moyen(i,k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+ 469 $ rhow(i,k)*epaisseur1(i,k))/epaisseur2(i,k) 470 ENDDO 471 ENDDO 472 473 C 474 C Choose an integration bound well above wake top 475 c----------------------------------------------------------------- 476 c 477 C Pupper = 50000. ! melting level 478 c Pupper = 60000. 479 c Pupper = 80000. ! essais pour case_e 480 DO i = 1,klon 481 Pupper(i) = 0.6*ph(i,1) 482 Pupper(i) = max(Pupper(i), 45000.) 483 ccc Pupper(i) = 60000. 484 ENDDO 485 486 C 487 C Determine Wake top pressure (Ptop) from buoyancy integral 488 C -------------------------------------------------------- 489 c 490 c-1/ Pressure of the level where dth becomes less than delta_t_min. 491 492 DO i=1,klon 493 ptop_provis(i)=ph(i,1) 494 ENDDO 495 DO k= 2,klev 496 DO i=1,klon 497 c 498 cIM v3JYG; ptop_provis(i).LT. ph(i,1) 499 c 500 IF (dth(i,k) .GT. -delta_t_min .and. 501 $ dth(i,k-1).LT. -delta_t_min .and. 502 $ ptop_provis(i).EQ. ph(i,1)) THEN 503 ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) 504 $ - (dth(i,k-1)+delta_t_min)*p(i,k)) / 505 $ (dth(i,k) - dth(i,k-1)) 506 ENDIF 507 ENDDO 508 ENDDO 509 510 c-2/ dth integral 511 512 DO i=1,klon 513 sum_dth(i) = 0. 514 dthmin(i) = -delta_t_min 515 z(i) = 0. 516 ENDDO 517 518 DO k = 1,klev 519 DO i=1,klon 520 dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg) 521 IF (dz(i) .gt. 0) THEN 522 z(i) = z(i)+dz(i) 523 sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i) 524 dthmin(i) = amin1(dthmin(i),dth(i,k)) 525 ENDIF 526 ENDDO 527 ENDDO 528 529 c-3/ height of triangle with area= sum_dth and base = dthmin 530 531 DO i=1,klon 532 hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5) 533 hw0(i) = amax1(hwmin,hw0(i)) 534 ENDDO 535 536 c-4/ now, get Ptop 537 538 DO i=1,klon 539 z(i) = 0. 540 ptop(i) = ph(i,1) 541 ENDDO 542 543 DO k = 1,klev 544 DO i=1,klon 545 dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg),hw0(i)-z(i)) 546 IF (dz(i) .gt. 0) THEN 547 z(i) = z(i)+dz(i) 548 ptop(i) = ph(i,k)-rho(i,k)*rg*dz(i) 549 ENDIF 550 ENDDO 551 ENDDO 552 553 554 C-5/ Determination de ktop et kupper 555 556 DO k=klev,1,-1 557 DO i=1,klon 558 IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k 559 IF (ph(i,k+1) .lt. pupper(i)) kupper(i)=k 560 ENDDO 561 ENDDO 562 563 c On evite kupper = 1 564 DO i=1,klon 565 kupper(i) = max(kupper(i),2) 566 ENDDO 567 568 569 c-6/ Correct ktop and ptop 570 571 DO i = 1,klon 572 ptop_new(i)=ptop(i) 573 ENDDO 574 DO k= klev,2,-1 575 DO i=1,klon 576 IF (k .LE. ktop(i) .and. 577 $ ptop_new(i) .EQ. ptop(i) .and. 578 $ dth(i,k) .GT. -delta_t_min .and. 579 $ dth(i,k-1).LT. -delta_t_min) THEN 580 ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) 581 $ - (dth(i,k-1)+delta_t_min)*p(i,k)) / 582 $ (dth(i,k) - dth(i,k-1)) 583 ENDIF 584 ENDDO 585 ENDDO 586 587 DO i=1,klon 588 ptop(i) = ptop_new(i) 589 ENDDO 590 591 DO k=klev,1,-1 592 DO i=1,klon 593 IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k 594 ENDDO 595 ENDDO 596 c 597 c-5/ Set deltatw & deltaqw to 0 above kupper 598 c 599 DO k = 1,klev 600 DO i=1,klon 601 IF (k.GE. kupper(i)) THEN 602 deltatw(i,k) = 0. 603 deltaqw(i,k) = 0. 604 ENDIF 605 ENDDO 606 ENDDO 607 c 608 C 609 C Vertical gradient of LS omega 610 C 611 DO k = 1,klev 612 DO i=1,klon 613 IF (k.LE. kupper(i)) THEN 614 dp_omgb(i,k) = (omgb(i,k+1) - omgb(i,k))/(ph(i,k+1)-ph(i,k)) 615 ENDIF 616 ENDDO 617 ENDDO 618 C 619 C Integrals (and wake top level number) 620 C -------------------------------------- 621 C 622 C Initialize sum_thvu to 1st level virt. pot. temp. 623 624 DO i=1,klon 625 z(i) = 1. 626 dz(i) = 1. 627 sum_thvu(i) = thu(i,1)*(1.+eps*qu(i,1))*dz(i) 628 sum_dth(i) = 0. 629 ENDDO 630 631 DO k = 1,klev 632 DO i=1,klon 633 dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) 634 IF (dz(i) .GT. 0) THEN 635 z(i) = z(i)+dz(i) 636 sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i) 637 sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i) 638 sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i) 639 sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i) 640 sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i) 641 sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i) 642 sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i) 643 sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i) 644 sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i) 645 ENDIF 646 ENDDO 647 ENDDO 648 c 649 DO i=1,klon 650 hw0(i) = z(i) 651 ENDDO 652 c 653 C 654 C 2.1 - WAPE and mean forcing computation 655 C --------------------------------------- 656 C 657 C --------------------------------------- 658 C 659 C Means 660 661 DO i=1,klon 662 av_thu(i) = sum_thu(i)/hw0(i) 663 av_tu(i) = sum_tu(i)/hw0(i) 664 av_qu(i) = sum_qu(i)/hw0(i) 665 av_thvu(i) = sum_thvu(i)/hw0(i) 666 c av_thve = sum_thve/hw0 667 av_dth(i) = sum_dth(i)/hw0(i) 668 av_dq(i) = sum_dq(i)/hw0(i) 669 av_rho(i) = sum_rho(i)/hw0(i) 670 av_dtdwn(i) = sum_dtdwn(i)/hw0(i) 671 av_dqdwn(i) = sum_dqdwn(i)/hw0(i) 672 673 wape(i) = - rg*hw0(i)*(av_dth(i) 674 $ + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)* 675 $ av_dq(i) ))/av_thvu(i) 676 ENDDO 677 C 678 C 2.2 Prognostic variable update 679 C ------------------------------ 680 C 681 C Filter out bad wakes 682 683 DO k = 1,klev 684 DO i=1,klon 685 IF ( wape(i) .LT. 0.) THEN 686 deltatw(i,k) = 0. 687 deltaqw(i,k) = 0. 688 dth(i,k) = 0. 689 ENDIF 690 ENDDO 691 ENDDO 692 c 693 DO i=1,klon 694 IF ( wape(i) .LT. 0.) THEN 695 wape(i) = 0. 696 Cstar(i) = 0. 697 hw(i) = hwmin 698 sigmaw(i) = amax1(sigmad,sigd_con(i)) 699 fip(i) = 0. 700 gwake(i) = .FALSE. 701 ELSE 702 Cstar(i) = stark*sqrt(2.*wape(i)) 703 gwake(i) = .TRUE. 704 ENDIF 705 ENDDO 706 707 c 708 c Check qx and qw positivity 709 c -------------------------- 710 DO i = 1,klon 711 q0_min(i)=min( (qe(i,1)-sigmaw(i)*deltaqw(i,1)), 712 $ (qe(i,1)+(1.-sigmaw(i))*deltaqw(i,1)) ) 713 ENDDO 714 DO k = 2,klev 715 DO i = 1,klon 716 q1_min(i)=min( (qe(i,k)-sigmaw(i)*deltaqw(i,k)), 717 $ (qe(i,k)+(1.-sigmaw(i))*deltaqw(i,k)) ) 718 IF (q1_min(i).le.q0_min(i)) THEN 719 q0_min(i)=q1_min(i) 720 ENDIF 721 ENDDO 722 ENDDO 723 c 724 DO i = 1,klon 725 OK_qx_qw(i) = q0_min(i) .GE. 0. 726 alpha(i) = 1. 727 ENDDO 728 c 729 CC ----------------------------------------------------------------- 730 C Sub-time-stepping 731 C ----------------- 732 C 733 nsub=10 734 dtimesub=dtime/nsub 735 c 736 c------------------------------------------------------------ 737 DO isubstep = 1,nsub 738 c------------------------------------------------------------ 739 c 740 c wk_adv is the logical flag enabling wake evolution in the time advance loop 741 DO i = 1,klon 742 wk_adv(i) = OK_qx_qw(i) .AND. alpha(i) .GE. 1. 743 ENDDO 744 c 745 ccc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement négatif de ktop à kupper -------- 746 ccc On calcule pour cela une densité wdens0 pour laquelle on aurait un entrainement nul --- 747 DO i=1,klon 748 cc print *,' isubstep,wk_adv(i),cstar(i),wape(i) ', 749 cc $ isubstep,wk_adv(i),cstar(i),wape(i) 750 IF (wk_adv(i) .AND. cstar(i).GT.0.01) THEN 751 omg(i,kupper(i)+1) = - Rg*amdwn(i,kupper(i)+1)/sigmaw(i) 752 $ + Rg*amup(i,kupper(i)+1)/(1.-sigmaw(i)) 753 wdens0 = ( sigmaw(i) / (4.*3.14) ) * 754 $ ( (1.-sigmaw(i)) * omg(i,kupper(i)+1) / 755 $ ( (ph(i,1)-pupper(i)) * cstar(i) ) ) **(2) 756 IF ( wdens(i) .LE. wdens0*1.1 ) THEN 757 wdens(i) = wdens0 758 ENDIF 759 cc print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i) 760 cc $ ,ph(i,1)-pupper(i)', 761 cc $ omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i) 762 cc $ ,ph(i,1)-pupper(i) 763 ENDIF 764 ENDDO 765 766 ccc nrlmd 767 768 DO i=1,klon 769 IF (wk_adv(i)) THEN 770 gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i)) 771 sigmaw(i)=amin1(sigmaw(i),sigmaw_max) 772 ENDIF 773 ENDDO 774 DO i=1,klon 775 IF (wk_adv(i)) THEN 776 ccc nrlmd Introduction du taux de mortalité des poches et test sur sigmaw_max=0.4 777 ccc d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub 778 IF (sigmaw(i).ge.sigmaw_max) THEN 779 death_rate(i)=gfl(i)*Cstar(i)/sigmaw(i) 780 ELSE 781 death_rate(i)=0. 782 END IF 783 d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub 784 $ - death_rate(i)*sigmaw(i)*dtimesub 785 c $ - nat_rate(i)*sigmaw(i)*dtimesub 786 cc print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i), 787 cc $ death_rate(i),ktop(i),kupper(i)', 788 cc $ d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i), 789 cc $ death_rate(i),ktop(i),kupper(i) 790 791 c sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub 792 c sigmaw(i) =min(sigmaw(i),0.99) !!!!!!!! 793 c wdens = wdens0/(10.*sigmaw) 794 c sigmaw =max(sigmaw,sigd_con) 795 c sigmaw =max(sigmaw,sigmad) 796 ENDIF 797 ENDDO 798 C 799 C 800 c calcul de la difference de vitesse verticale poche - zone non perturbee 801 cIM 060208 differences par rapport au code initial; init. a 0 dp_deltomg 802 cIM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit 803 cIM 060208 au niveau k=1..? 804 DO k= 1,klev 805 DO i = 1,klon 806 if (wk_adv(i)) THEN !!! nrlmd 807 dp_deltomg(i,k)=0. 808 end if 809 ENDDO 810 ENDDO 811 DO k= 1,klev+1 812 DO i = 1,klon 813 if (wk_adv(i)) THEN !!! nrlmd 814 omg(i,k)=0. 815 end if 816 ENDDO 817 ENDDO 818 c 819 DO i=1,klon 820 IF (wk_adv(i)) THEN 821 z(i)= 0. 822 omg(i,1) = 0. 823 dp_deltomg(i,1) = -(gfl(i)*Cstar(i))/(sigmaw(i) * (1-sigmaw(i))) 824 ENDIF 825 ENDDO 826 c 827 DO k= 2,klev 828 DO i = 1,klon 829 IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN 830 dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg) 831 z(i) = z(i)+dz(i) 832 dp_deltomg(i,k)= dp_deltomg(i,1) 833 omg(i,k)= dp_deltomg(i,1)*z(i) 834 ENDIF 835 ENDDO 836 ENDDO 837 c 838 DO i = 1,klon 839 IF (wk_adv(i)) THEN 840 dztop(i)=-(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg) 841 ztop(i) = z(i)+dztop(i) 842 omgtop(i)=dp_deltomg(i,1)*ztop(i) 843 ENDIF 844 ENDDO 845 c 846 c ----------------- 847 c From m/s to Pa/s 848 c ----------------- 849 c 850 DO i=1,klon 851 IF (wk_adv(i)) THEN 852 omgtop(i) = -rho(i,ktop(i))*rg*omgtop(i) 853 dp_deltomg(i,1) = omgtop(i)/(ptop(i)-ph(i,1)) 854 ENDIF 855 ENDDO 856 c 857 DO k= 1,klev 858 DO i = 1,klon 859 IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN 860 omg(i,k) = - rho(i,k)*rg*omg(i,k) 861 dp_deltomg(i,k) = dp_deltomg(i,1) 862 ENDIF 863 ENDDO 864 ENDDO 865 c 866 c raccordement lineaire de omg de ptop a pupper 867 868 DO i=1,klon 869 IF ( wk_adv(i) .AND. kupper(i) .GT. ktop(i)) THEN 870 omg(i,kupper(i)+1) = - Rg*amdwn(i,kupper(i)+1)/sigmaw(i) 871 $ + Rg*amup(i,kupper(i)+1)/(1.-sigmaw(i)) 872 dp_deltomg(i,kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ 873 $ (ptop(i)-pupper(i)) 874 ENDIF 875 ENDDO 876 c 877 cc DO i=1,klon 878 cc print*,'Pente entre 0 et kupper (référence)' 879 cc $ ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1)) 880 cc print*,'Pente entre ktop et kupper' 881 cc $ ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i)) 882 cc ENDDO 883 cc 884 DO k= 1,klev 885 DO i = 1,klon 886 IF( wk_adv(i) .AND. k .GT. ktop(i) .AND. k .LE. kupper(i)) THEN 887 dp_deltomg(i,k) = dp_deltomg(i,kupper(i)) 888 omg(i,k) = omgtop(i)+(ph(i,k)-ptop(i))*dp_deltomg(i,kupper(i)) 889 ENDIF 890 ENDDO 891 ENDDO 892 ccc nrlmd 893 cc DO i=1,klon 894 cc print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1) 895 cc END DO 896 ccc 897 c 898 c 899 c-- Compute wake average vertical velocity omgbw 900 c 901 c 902 DO k = 1,klev+1 903 DO i=1,klon 904 IF ( wk_adv(i)) THEN 905 omgbw(i,k) = omgb(i,k)+(1.-sigmaw(i))*omg(i,k) 906 ENDIF 907 ENDDO 908 ENDDO 909 c-- and its vertical gradient dp_omgbw 910 c 911 DO k = 1,klev 912 DO i=1,klon 913 IF ( wk_adv(i)) THEN 914 dp_omgbw(i,k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k)) 915 ENDIF 916 ENDDO 917 ENDDO 918 C 919 c-- Upstream coefficients for omgb velocity 920 c-- (alpha_up(k) is the coefficient of the value at level k) 921 c-- (1-alpha_up(k) is the coefficient of the value at level k-1) 922 DO k = 1,klev 923 DO i=1,klon 924 IF ( wk_adv(i)) THEN 925 alpha_up(i,k) = 0. 926 IF (omgb(i,k) .GT. 0.) alpha_up(i,k) = 1. 927 ENDIF 928 ENDDO 929 ENDDO 930 931 c Matrix expressing [The,deltatw] from [Th1,Th2] 932 933 DO i=1,klon 934 IF ( wk_adv(i)) THEN 935 RRe1(i) = 1.-sigmaw(i) 936 RRe2(i) = sigmaw(i) 937 ENDIF 938 ENDDO 939 RRd1 = -1. 940 RRd2 = 1. 941 c 942 c-- Get [Th1,Th2], dth and [q1,q2] 943 c 944 DO k= 1,klev 945 DO i = 1,klon 946 IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN 947 dth(i,k) = deltatw(i,k)/ppi(i,k) 948 Th1(i,k) = the(i,k) - sigmaw(i) *dth(i,k) ! undisturbed area 949 Th2(i,k) = the(i,k) + (1.-sigmaw(i))*dth(i,k) ! wake 950 q1(i,k) = qe(i,k) - sigmaw(i) *deltaqw(i,k) ! undisturbed area 951 q2(i,k) = qe(i,k) + (1.-sigmaw(i))*deltaqw(i,k) ! wake 952 ENDIF 953 ENDDO 954 ENDDO 955 956 DO i=1,klon 957 if (wk_adv(i)) then !!! nrlmd 958 D_Th1(i,1) = 0. 959 D_Th2(i,1) = 0. 960 D_dth(i,1) = 0. 961 D_q1(i,1) = 0. 962 D_q2(i,1) = 0. 963 D_dq(i,1) = 0. 964 end if 965 ENDDO 966 967 DO k= 2,klev 968 DO i = 1,klon 969 IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN 970 D_Th1(i,k) = Th1(i,k-1)-Th1(i,k) 971 D_Th2(i,k) = Th2(i,k-1)-Th2(i,k) 972 D_dth(i,k) = dth(i,k-1)-dth(i,k) 973 D_q1(i,k) = q1(i,k-1)-q1(i,k) 974 D_q2(i,k) = q2(i,k-1)-q2(i,k) 975 D_dq(i,k) = deltaqw(i,k-1)-deltaqw(i,k) 976 ENDIF 977 ENDDO 978 ENDDO 979 980 DO i=1,klon 981 IF( wk_adv(i)) THEN 982 omgbdth(i,1) = 0. 983 omgbdq(i,1) = 0. 984 ENDIF 985 ENDDO 986 987 DO k= 2,klev 988 DO i = 1,klon 989 IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN ! loop on interfaces 990 omgbdth(i,k) = omgb(i,k)*( dth(i,k-1) - dth(i,k)) 991 omgbdq(i,k) = omgb(i,k)*(deltaqw(i,k-1) - deltaqw(i,k)) 992 ENDIF 993 ENDDO 994 ENDDO 995 c 996 c----------------------------------------------------------------- 997 DO k= 1,klev 998 DO i = 1,klon 999 IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN 1000 c----------------------------------------------------------------- 1001 c 1002 c Compute redistribution (advective) term 1003 c 1004 d_deltatw(i,k) = 1005 $ dtimesub/(Ph(i,k)-Ph(i,k+1))*( 1006 $ RRd1*omg(i,k )*sigmaw(i) *D_Th1(i,k) 1007 $ -RRd2*omg(i,k+1)*(1.-sigmaw(i))*D_Th2(i,k+1) 1008 $ -(1.-alpha_up(i,k))*omgbdth(i,k) - alpha_up(i,k+1)* 1009 $ omgbdth(i,k+1))*ppi(i,k) 1010 c print*,'d_deltatw=',d_deltatw(i,k) 1011 c 1012 d_deltaqw(i,k) = 1013 $ dtimesub/(Ph(i,k)-Ph(i,k+1))*( 1014 $ RRd1*omg(i,k )*sigmaw(i) *D_q1(i,k) 1015 $ -RRd2*omg(i,k+1)*(1.-sigmaw(i))*D_q2(i,k+1) 1016 $ -(1.-alpha_up(i,k))*omgbdq(i,k) - alpha_up(i,k+1)* 1017 $ omgbdq(i,k+1)) 1018 c print*,'d_deltaqw=',d_deltaqw(i,k) 1019 c 1020 c and increment large scale tendencies 1021 c 1022 1023 c 1024 C 1025 CC ----------------------------------------------------------------- 1026 d_te(i,k) = dtimesub*( 1027 $ ( RRe1(i)*omg(i,k )*sigmaw(i) *D_Th1(i,k) 1028 $ -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_Th2(i,k+1) ) 1029 $ /(Ph(i,k)-Ph(i,k+1)) 1030 ccc nrlmd $ -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*dp_deltomg(i,k) 1031 $ -sigmaw(i)*(1.-sigmaw(i))*dth(i,k) 1032 $ *(omg(i,k)-omg(i,k+1))/(Ph(i,k)-Ph(i,k+1)) 1033 ccc 1034 $ )*ppi(i,k) 1035 c 1036 d_qe(i,k) = dtimesub*( 1037 $ ( RRe1(i)*omg(i,k )*sigmaw(i) *D_q1(i,k) 1038 $ -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_q2(i,k+1) ) 1039 $ /(Ph(i,k)-Ph(i,k+1)) 1040 ccc nrlmd $ -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*dp_deltomg(i,k) 1041 $ -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k) 1042 $ *(omg(i,k)-omg(i,k+1))/(Ph(i,k)-Ph(i,k+1)) 1043 ccc 1044 $ ) 1045 ccc nrlmd 1046 ELSE IF(wk_adv(i) .AND. k .EQ. kupper(i)) THEN 1047 d_te(i,k) = dtimesub*( 1048 $ ( RRe1(i)*omg(i,k )*sigmaw(i) *D_Th1(i,k) 1049 $ /(Ph(i,k)-Ph(i,k+1))) 1050 $ )*ppi(i,k) 1051 1052 d_qe(i,k) = dtimesub*( 1053 $ ( RRe1(i)*omg(i,k )*sigmaw(i) *D_q1(i,k) 1054 $ /(Ph(i,k)-Ph(i,k+1))) 1055 $ ) 1056 1057 ENDIF 1058 ccc 1059 ENDDO 1060 ENDDO 1061 c------------------------------------------------------------------ 1062 C 1063 C Increment state variables 1064 1065 DO k= 1,klev 1066 DO i = 1,klon 1067 ccc nrlmd IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN 1068 IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN 1069 ccc 1070 1071 1072 c 1073 c Coefficient de répartition 1074 1075 Crep(i,k)=Crep_sol*(ph(i,kupper(i))-ph(i,k))/(ph(i,kupper(i)) 1076 $ -ph(i,1)) 1077 Crep(i,k)=Crep(i,k)+Crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)- 1078 $ ph(i,kupper(i))) 1079 1080 1081 c Reintroduce compensating subsidence term. 1082 1083 c dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw 1084 c dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)) 1085 c . /(1-sigmaw) 1086 c dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw 1087 c dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)) 1088 c . /(1-sigmaw) 1089 c 1090 c dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw 1091 c dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k)) 1092 c . /(1-sigmaw) 1093 c dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw 1094 c dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k)) 1095 c . /(1-sigmaw) 1096 1097 dtKE(i,k)=(dtdwn(i,k)/sigmaw(i) - dta(i,k)/(1.-sigmaw(i))) 1098 dqKE(i,k)=(dqdwn(i,k)/sigmaw(i) - dqa(i,k)/(1.-sigmaw(i))) 1099 c print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k) 1100 c 1101 dtPBL(i,k)=(wdtPBL(i,k)/sigmaw(i) - udtPBL(i,k)/(1.-sigmaw(i))) 1102 dqPBL(i,k)=(wdqPBL(i,k)/sigmaw(i) - udqPBL(i,k)/(1.-sigmaw(i))) 1103 c print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k) 1104 c 1105 ccc nrlmd Prise en compte du taux de mortalité 1106 ccc Définitions de entr, detr 1107 detr(i,k)=0. 1108 1109 entr(i,k)=detr(i,k)+gfl(i)*cstar(i)+ 1110 $ sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i,k) 1111 1112 spread(i,k) = (entr(i,k)-detr(i,k))/sigmaw(i) 1113 ccc spread(i,k) = (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/ 1114 ccc $ sigmaw(i) 1115 1116 1117 c ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU Jingmei 1118 1119 ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k), 1120 ! & Tgw(i,k),deltatw(i,k) 1121 d_deltat_gw(i,k)=d_deltat_gw(i,k)-Tgw(i,k)*deltatw(i,k)* 1122 $ dtimesub 1123 ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k) 1124 ff(i)=d_deltatw(i,k)/dtimesub 1125 1126 c Sans GW 1127 c 1128 c deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k)) 1129 c 1130 c GW formule 1 1131 c 1132 c deltatw(k) = deltatw(k)+dtimesub* 1133 c $ (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k)) 1134 c 1135 c GW formule 2 1136 1137 IF (dtimesub*Tgw(i,k).lt.1.e-10) THEN 1138 d_deltatw(i,k) = dtimesub* 1139 $ (ff(i)+dtKE(i,k)+dtPBL(i,k) 1140 ccc $ -spread(i,k)*deltatw(i,k) 1141 $ - entr(i,k)*deltatw(i,k)/sigmaw(i) 1142 $ - (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k) 1143 $ / (1.-sigmaw(i)) 1144 ccc 1145 $ -Tgw(i,k)*deltatw(i,k)) 1146 ELSE 1147 d_deltatw(i,k) = 1/Tgw(i,k)*(1-exp(-dtimesub* 1148 $ Tgw(i,k)))* 1149 $ (ff(i)+dtKE(i,k)+dtPBL(i,k) 1150 ccc $ -spread(i,k)*deltatw(i,k) 1151 $ - entr(i,k)*deltatw(i,k)/sigmaw(i) 1152 $ - (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k) 1153 $ / (1.-sigmaw(i)) 1154 ccc 1155 $ -Tgw(i,k)*deltatw(i,k)) 1156 ENDIF 1157 1158 dth(i,k) = deltatw(i,k)/ppi(i,k) 1159 1160 gg(i)=d_deltaqw(i,k)/dtimesub 1161 1162 d_deltaqw(i,k) = dtimesub*(gg(i)+ dqKE(i,k)+dqPBL(i,k) 1163 ccc $ -spread(i,k)*deltaqw(i,k)) 1164 $ - entr(i,k)*deltaqw(i,k)/sigmaw(i) 1165 $ - (death_rate(i)*sigmaw(i)+detr(i,k))*deltaqw(i,k) 1166 $ /(1.-sigmaw(i))) 1167 ccc 1168 1169 ccc nrlmd 1170 ccc d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k) 1171 ccc d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k) 1172 ccc 1173 ENDIF 1174 ENDDO 1175 ENDDO 1176 1177 C 1178 C Scale tendencies so that water vapour remains positive in w and x. 1179 C 1180 call wake_vec_modulation(klon,klev,wk_adv,epsilon,qe,d_qe,deltaqw, 1181 $ d_deltaqw,sigmaw,d_sigmaw,alpha) 1182 c 1183 ccc nrlmd 1184 cc print*,'alpha' 1185 cc do i=1,klon 1186 cc print*,alpha(i) 1187 cc end do 1188 ccc 1189 DO k = 1,klev 1190 DO i = 1,klon 1191 IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN 1192 d_te(i,k)=alpha(i)*d_te(i,k) 1193 d_qe(i,k)=alpha(i)*d_qe(i,k) 1194 d_deltatw(i,k)=alpha(i)*d_deltatw(i,k) 1195 d_deltaqw(i,k)=alpha(i)*d_deltaqw(i,k) 1196 d_deltat_gw(i,k)=alpha(i)*d_deltat_gw(i,k) 1197 ENDIF 1198 ENDDO 1199 ENDDO 1200 DO i = 1,klon 1201 IF( wk_adv(i)) THEN 1202 d_sigmaw(i)=alpha(i)*d_sigmaw(i) 1203 ENDIF 1204 ENDDO 1205 1206 C Update large scale variables and wake variables 1207 cIM 060208 manque DO i + remplace DO k=1,kupper(i) 1208 cIM 060208 DO k = 1,kupper(i) 1209 DO k= 1,klev 1210 DO i = 1,klon 1211 IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN 1212 dtls(i,k)=dtls(i,k)+d_te(i,k) 1213 dqls(i,k)=dqls(i,k)+d_qe(i,k) 1214 ccc nrlmd 1215 d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k) 1216 d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k) 1217 ccc 1218 ENDIF 1219 ENDDO 1220 ENDDO 1221 DO k= 1,klev 1222 DO i = 1,klon 1223 IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN 1224 te(i,k) = te0(i,k) + dtls(i,k) 1225 qe(i,k) = qe0(i,k) + dqls(i,k) 1226 the(i,k) = te(i,k)/ppi(i,k) 1227 deltatw(i,k) = deltatw(i,k)+d_deltatw(i,k) 1228 deltaqw(i,k) = deltaqw(i,k)+d_deltaqw(i,k) 1229 dth(i,k) = deltatw(i,k)/ppi(i,k) 1230 cc print*,'k,qx,qw',k,qe(i,k)-sigmaw(i)*deltaqw(i,k) 1231 cc $ ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k) 1232 ENDIF 1233 ENDDO 1234 ENDDO 1235 DO i = 1,klon 1236 IF( wk_adv(i)) THEN 1237 sigmaw(i) = sigmaw(i)+d_sigmaw(i) 1238 ENDIF 1239 ENDDO 1240 c 1241 C 1242 c Determine Ptop from buoyancy integral 1243 c --------------------------------------- 1244 c 1245 c- 1/ Pressure of the level where dth changes sign. 1246 c 1247 DO i=1,klon 1248 IF ( wk_adv(i)) THEN 1249 Ptop_provis(i)=ph(i,1) 1250 ENDIF 1251 ENDDO 1252 c 1253 DO k= 2,klev 1254 DO i=1,klon 1255 IF ( wk_adv(i) .AND. 1256 $ Ptop_provis(i) .EQ. ph(i,1) .AND. 1257 $ dth(i,k) .GT. -delta_t_min .and. 1258 $ dth(i,k-1).LT. -delta_t_min) THEN 1259 Ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) 1260 $ - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k) 1261 $ - dth(i,k-1)) 1262 ENDIF 1263 ENDDO 1264 ENDDO 1265 c 1266 c- 2/ dth integral 1267 c 1268 DO i=1,klon 1269 if (wk_adv(i)) then !!! nrlmd 1270 sum_dth(i) = 0. 1271 dthmin(i) = -delta_t_min 1272 z(i) = 0. 1273 end if 1274 ENDDO 1275 1276 DO k = 1,klev 1277 DO i=1,klon 1278 IF ( wk_adv(i)) THEN 1279 dz(i) = -(amax1(ph(i,k+1),Ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg) 1280 IF (dz(i) .gt. 0) THEN 1281 z(i) = z(i)+dz(i) 1282 sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i) 1283 dthmin(i) = amin1(dthmin(i),dth(i,k)) 1284 ENDIF 1285 ENDIF 1286 ENDDO 1287 ENDDO 1288 c 1289 c- 3/ height of triangle with area= sum_dth and base = dthmin 1290 1291 DO i=1,klon 1292 IF ( wk_adv(i)) THEN 1293 hw(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5) 1294 hw(i) = amax1(hwmin,hw(i)) 1295 ENDIF 1296 ENDDO 1297 c 1298 c- 4/ now, get Ptop 1299 c 1300 DO i=1,klon 1301 if (wk_adv(i)) then !!! nrlmd 1302 ktop(i) = 0 1303 z(i)=0. 1304 end if 1305 ENDDO 1306 c 1307 DO k = 1,klev 1308 DO i=1,klon 1309 IF ( wk_adv(i)) THEN 1310 dz(i) = amin1(-(ph(i,k+1)-Ph(i,k))/(rho(i,k)*rg),hw(i)-z(i)) 1311 IF (dz(i) .gt. 0) THEN 1312 z(i) = z(i)+dz(i) 1313 Ptop(i) = Ph(i,k)-rho(i,k)*rg*dz(i) 1314 ktop(i) = k 1315 ENDIF 1316 ENDIF 1317 ENDDO 1318 ENDDO 1319 c 1320 c 4.5/Correct ktop and ptop 1321 c 1322 DO i=1,klon 1323 IF ( wk_adv(i)) THEN 1324 Ptop_new(i)=ptop(i) 1325 ENDIF 1326 ENDDO 1327 c 1328 DO k= klev,2,-1 1329 DO i=1,klon 1330 cIM v3JYG; IF (k .GE. ktop(i) 1331 IF ( wk_adv(i) .AND. 1332 $ k .LE. ktop(i) .AND. 1333 $ ptop_new(i) .EQ. ptop(i) .AND. 1334 $ dth(i,k) .GT. -delta_t_min .and. 1335 $ dth(i,k-1).LT. -delta_t_min) THEN 1336 Ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) 1337 $ - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k) 1338 $ - dth(i,k-1)) 1339 ENDIF 1340 ENDDO 1341 ENDDO 1342 c 1343 c 1344 DO i=1,klon 1345 IF ( wk_adv(i)) THEN 1346 ptop(i) = ptop_new(i) 1347 ENDIF 1348 ENDDO 1349 1350 DO k=klev,1,-1 1351 DO i=1,klon 1352 if (wk_adv(i)) then !!! nrlmd 1353 IF (ph(i,k+1) .LT. ptop(i)) ktop(i)=k 1354 end if 1355 ENDDO 1356 ENDDO 1357 c 1358 c 5/ Set deltatw & deltaqw to 0 above kupper 1359 c 1360 DO k = 1,klev 1361 DO i=1,klon 1362 IF ( wk_adv(i) .AND. k .GE. kupper(i)) THEN 1363 deltatw(i,k) = 0. 1364 deltaqw(i,k) = 0. 1365 ENDIF 1366 ENDDO 1367 ENDDO 1368 c 1369 C 1370 c-------------Cstar computation--------------------------------- 1371 DO i=1, klon 1372 if (wk_adv(i)) then !!! nrlmd 1373 sum_thu(i) = 0. 1374 sum_tu(i) = 0. 1375 sum_qu(i) = 0. 1376 sum_thvu(i) = 0. 1377 sum_dth(i) = 0. 1378 sum_dq(i) = 0. 1379 sum_rho(i) = 0. 1380 sum_dtdwn(i) = 0. 1381 sum_dqdwn(i) = 0. 1382 1383 av_thu(i) = 0. 1384 av_tu(i) =0. 1385 av_qu(i) =0. 1386 av_thvu(i) = 0. 1387 av_dth(i) = 0. 1388 av_dq(i) = 0. 1389 av_rho(i) =0. 1390 av_dtdwn(i) =0. 1391 av_dqdwn(i) = 0. 1392 end if 1393 ENDDO 1394 C 1395 C Integrals (and wake top level number) 1396 C -------------------------------------- 1397 C 1398 C Initialize sum_thvu to 1st level virt. pot. temp. 1399 1400 DO i=1,klon 1401 if (wk_adv(i)) then !!! nrlmd 1402 z(i) = 1. 1403 dz(i) = 1. 1404 sum_thvu(i) = thu(i,1)*(1.+eps*qu(i,1))*dz(i) 1405 sum_dth(i) = 0. 1406 end if 1407 ENDDO 1408 1409 DO k = 1,klev 1410 DO i=1,klon 1411 if (wk_adv(i)) then !!! nrlmd 1412 dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) 1413 IF (dz(i) .GT. 0) THEN 1414 z(i) = z(i)+dz(i) 1415 sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i) 1416 sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i) 1417 sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i) 1418 sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i) 1419 sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i) 1420 sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i) 1421 sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i) 1422 sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i) 1423 sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i) 1424 ENDIF 1425 end if 1426 ENDDO 1427 ENDDO 1428 c 1429 DO i=1,klon 1430 if (wk_adv(i)) then !!! nrlmd 1431 hw0(i) = z(i) 1432 end if 1433 ENDDO 1434 c 1435 C 1436 C - WAPE and mean forcing computation 1437 C --------------------------------------- 1438 C 1439 C --------------------------------------- 1440 C 1441 C Means 1442 1443 DO i=1,klon 1444 if (wk_adv(i)) then !!! nrlmd 1445 av_thu(i) = sum_thu(i)/hw0(i) 1446 av_tu(i) = sum_tu(i)/hw0(i) 1447 av_qu(i) = sum_qu(i)/hw0(i) 1448 av_thvu(i) = sum_thvu(i)/hw0(i) 1449 av_dth(i) = sum_dth(i)/hw0(i) 1450 av_dq(i) = sum_dq(i)/hw0(i) 1451 av_rho(i) = sum_rho(i)/hw0(i) 1452 av_dtdwn(i) = sum_dtdwn(i)/hw0(i) 1453 av_dqdwn(i) = sum_dqdwn(i)/hw0(i) 1454 c 1455 wape(i) = - rg*hw0(i)*(av_dth(i) 1456 $ + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)* 1457 $ av_dq(i) ))/av_thvu(i) 1458 end if 1459 ENDDO 1460 C 1461 C Filter out bad wakes 1462 1463 DO k = 1,klev 1464 DO i=1,klon 1465 if (wk_adv(i)) then !!! nrlmd 1466 IF ( wape(i) .LT. 0.) THEN 1467 deltatw(i,k) = 0. 1468 deltaqw(i,k) = 0. 1469 dth(i,k) = 0. 1470 ENDIF 1471 end if 1472 ENDDO 1473 ENDDO 1474 c 1475 DO i=1,klon 1476 if (wk_adv(i)) then !!! nrlmd 1477 IF ( wape(i) .LT. 0.) THEN 1478 wape(i) = 0. 1479 Cstar(i) = 0. 1480 hw(i) = hwmin 1481 sigmaw(i) = max(sigmad,sigd_con(i)) 1482 fip(i) = 0. 1483 gwake(i) = .FALSE. 1484 ELSE 1485 Cstar(i) = stark*sqrt(2.*wape(i)) 1486 gwake(i) = .TRUE. 1487 ENDIF 1488 end if 1489 ENDDO 1490 1491 ENDDO ! end sub-timestep loop 1492 C 1493 C ----------------------------------------------------------------- 1494 c Get back to tendencies per second 1495 c 1496 DO k = 1,klev 1497 DO i=1,klon 1498 1499 ccc nrlmd IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN 1500 IF ( OK_qx_qw(i) .AND. k .LE. kupper(i)) THEN 1501 ccc 1502 dtls(i,k) = dtls(i,k)/dtime 1503 dqls(i,k) = dqls(i,k)/dtime 1504 d_deltatw2(i,k)=d_deltatw2(i,k)/dtime 1505 d_deltaqw2(i,k)=d_deltaqw2(i,k)/dtime 1506 d_deltat_gw(i,k) = d_deltat_gw(i,k)/dtime 1507 cc print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k) 1508 cc $ ,death_rate(i)*sigmaw(i) 1509 ENDIF 1510 ENDDO 1511 ENDDO 1512 1513 c 1514 c---------------------------------------------------------- 1515 c Determine wake final state; recompute wape, cstar, ktop; 1516 c filter out bad wakes. 1517 c---------------------------------------------------------- 1518 c 1519 C 2.1 - Undisturbed area and Wake integrals 1520 C --------------------------------------------------------- 1521 1522 DO i=1,klon 1523 ccc nrlmd if (wk_adv(i)) then !!! nrlmd 1524 if (OK_qx_qw(i)) then 1525 ccc 1526 z(i) = 0. 1527 sum_thu(i) = 0. 1528 sum_tu(i) = 0. 1529 sum_qu(i) = 0. 1530 sum_thvu(i) = 0. 1531 sum_dth(i) = 0. 1532 sum_dq(i) = 0. 1533 sum_rho(i) = 0. 1534 sum_dtdwn(i) = 0. 1535 sum_dqdwn(i) = 0. 1536 1537 av_thu(i) = 0. 1538 av_tu(i) =0. 1539 av_qu(i) =0. 1540 av_thvu(i) = 0. 1541 av_dth(i) = 0. 1542 av_dq(i) = 0. 1543 av_rho(i) =0. 1544 av_dtdwn(i) =0. 1545 av_dqdwn(i) = 0. 1546 end if 1547 ENDDO 1548 C Potential temperatures and humidity 1549 c---------------------------------------------------------- 1550 1551 DO k =1,klev 1552 DO i=1,klon 1553 ccc nrlmd IF ( wk_adv(i)) THEN 1554 if (OK_qx_qw(i)) then 1555 ccc 1556 rho(i,k) = p(i,k)/(rd*te(i,k)) 1557 IF(k .eq. 1) THEN 1558 rhoh(i,k) = ph(i,k)/(rd*te(i,k)) 1559 zhh(i,k)=0 1560 ELSE 1561 rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1))) 1562 zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1) 1563 ENDIF 1564 the(i,k) = te(i,k)/ppi(i,k) 1565 thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k) 1566 tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i) 1567 qu(i,k) = qe(i,k) - deltaqw(i,k)*sigmaw(i) 1568 rhow(i,k) = p(i,k)/(rd*(te(i,k)+deltatw(i,k))) 1569 dth(i,k) = deltatw(i,k)/ppi(i,k) 1570 ENDIF 1571 ENDDO 1572 ENDDO 1573 1574 C Integrals (and wake top level number) 1575 C ----------------------------------------------------------- 1576 1577 C Initialize sum_thvu to 1st level virt. pot. temp. 1578 1579 DO i=1,klon 1580 ccc nrlmd IF ( wk_adv(i)) THEN 1581 if (OK_qx_qw(i)) then 1582 ccc 1583 z(i) = 1. 1584 dz(i) = 1. 1585 sum_thvu(i) = thu(i,1)*(1.+eps*qu(i,1))*dz(i) 1586 sum_dth(i) = 0. 1587 ENDIF 1588 ENDDO 1589 1590 DO k = 1,klev 1591 DO i=1,klon 1592 ccc nrlmd IF ( wk_adv(i)) THEN 1593 if (OK_qx_qw(i)) then 1594 ccc 1595 dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) 1596 IF (dz(i) .GT. 0) THEN 1597 z(i) = z(i)+dz(i) 1598 sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i) 1599 sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i) 1600 sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i) 1601 sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i) 1602 sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i) 1603 sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i) 1604 sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i) 1605 sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i) 1606 sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i) 1607 ENDIF 1608 ENDIF 1609 ENDDO 1610 ENDDO 1611 c 1612 DO i=1,klon 1613 ccc nrlmd IF ( wk_adv(i)) THEN 1614 if (OK_qx_qw(i)) then 1615 ccc 1616 hw0(i) = z(i) 1617 ENDIF 1618 ENDDO 1619 c 1620 C - WAPE and mean forcing computation 1621 C------------------------------------------------------------- 1622 1623 C Means 1624 1625 DO i=1, klon 1626 ccc nrlmd IF ( wk_adv(i)) THEN 1627 if (OK_qx_qw(i)) then 1628 ccc 1629 av_thu(i) = sum_thu(i)/hw0(i) 1630 av_tu(i) = sum_tu(i)/hw0(i) 1631 av_qu(i) = sum_qu(i)/hw0(i) 1632 av_thvu(i) = sum_thvu(i)/hw0(i) 1633 av_dth(i) = sum_dth(i)/hw0(i) 1634 av_dq(i) = sum_dq(i)/hw0(i) 1635 av_rho(i) = sum_rho(i)/hw0(i) 1636 av_dtdwn(i) = sum_dtdwn(i)/hw0(i) 1637 av_dqdwn(i) = sum_dqdwn(i)/hw0(i) 1638 1639 wape2(i) = - rg*hw0(i)*(av_dth(i) 1640 $ + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i) 1641 $ + av_dth(i)*av_dq(i) ))/av_thvu(i) 1642 ENDIF 1643 ENDDO 1644 1645 C Prognostic variable update 1646 C ------------------------------------------------------------ 1647 1648 C Filter out bad wakes 1649 c 1650 DO k = 1,klev 1651 DO i=1,klon 1652 ccc nrlmd IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN 1653 if (OK_qx_qw(i) .AND. wape2(i) .LT. 0.) then 1654 ccc 1655 deltatw(i,k) = 0. 1656 deltaqw(i,k) = 0. 1657 dth(i,k) = 0. 1658 ENDIF 1659 ENDDO 1660 ENDDO 1661 c 1662 1663 DO i=1, klon 1664 ccc nrlmd IF ( wk_adv(i)) THEN 1665 if (OK_qx_qw(i)) then 1666 ccc 1667 IF ( wape2(i) .LT. 0.) THEN 1668 wape2(i) = 0. 1669 Cstar2(i) = 0. 1670 hw(i) = hwmin 1671 sigmaw(i) = amax1(sigmad,sigd_con(i)) 1672 fip(i) = 0. 1673 gwake(i) = .FALSE. 1674 ELSE 1675 if(prt_level.ge.10) print*,'wape2>0' 1676 Cstar2(i) = stark*sqrt(2.*wape2(i)) 1677 gwake(i) = .TRUE. 1678 ENDIF 1679 ENDIF 1680 ENDDO 1681 c 1682 DO i=1, klon 1683 ccc nrlmd IF ( wk_adv(i)) THEN 1684 if (OK_qx_qw(i)) then 1685 ccc 1686 ktopw(i) = ktop(i) 1687 ENDIF 1688 ENDDO 1689 c 1690 DO i=1, klon 1691 ccc nrlmd IF ( wk_adv(i)) THEN 1692 if (OK_qx_qw(i)) then 1693 ccc 1694 IF (ktopw(i) .gt. 0 .and. gwake(i)) then 1695 1696 Cjyg1 Utilisation d'un h_efficace constant ( ~ feeding layer) 1697 ccc heff = 600. 1698 C Utilisation de la hauteur hw 1699 cc heff = 0.7*hw 1700 heff(i) = hw(i) 1701 1702 FIP(i) = 0.5*rho(i,ktopw(i))*Cstar2(i)**3*heff(i)*2* 1703 $ sqrt(sigmaw(i)*wdens(i)*3.14) 1704 FIP(i) = alpk * FIP(i) 1705 Cjyg2 1706 ELSE 1707 FIP(i) = 0. 1708 ENDIF 1709 ENDIF 1710 ENDDO 1711 c 1712 C Limitation de sigmaw 1713 1714 ccc nrlmd 1715 c DO i=1,klon 1716 c IF (OK_qx_qw(i)) THEN 1717 c IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max 1718 c ENDIF 1719 c ENDDO 1720 ccc 1721 DO k = 1,klev 1722 DO i=1, klon 1723 1724 ccc nrlmd On maintient désormais constant sigmaw en régime permanent 1725 ccc IF ((sigmaw(i).GT.sigmaw_max).or. 1726 IF ( ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or. 1727 $ (ktopw(i).le.2) .OR. 1728 $ .not. OK_qx_qw(i) ) THEN 1729 ccc 1730 dtls(i,k) = 0. 1731 dqls(i,k) = 0. 1732 deltatw(i,k) = 0. 1733 deltaqw(i,k) = 0. 1734 ENDIF 1735 ENDDO 1736 ENDDO 1737 c 1738 ccc nrlmd On maintient désormais constant sigmaw en régime permanent 1739 DO i=1, klon 1740 IF ( ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or. 1741 $ (ktopw(i).le.2) .OR. 1742 $ .not. OK_qx_qw(i) ) THEN 1743 wape(i) = 0. 1744 cstar(i)=0. 1745 hw(i) = hwmin 1746 sigmaw(i) = sigmad 1747 fip(i) = 0. 1748 ELSE 1749 wape(i) = wape2(i) 1750 cstar(i)=cstar2(i) 1751 ENDIF 1752 cc print*,'wape wape2 ktopw OK_qx_qw =', 1753 cc $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i) 1754 ENDDO 1755 c 1756 c 1757 RETURN 1758 END 1759 1760 SUBROUTINE wake_vec_modulation(nlon,nl,wk_adv,epsilon,qe,d_qe, 1761 $ deltaqw,d_deltaqw,sigmaw,d_sigmaw,alpha) 1762 c------------------------------------------------------ 1763 cDtermination du coefficient alpha tel que les tendances 1764 c corriges alpha*d_G, pour toutes les grandeurs G, correspondent 1765 c a une humidite positive dans la zone (x) et dans la zone (w). 1766 c------------------------------------------------------ 1767 c 1768 1769 c Input 1770 REAL qe(nlon,nl),d_qe(nlon,nl) 1771 REAL deltaqw(nlon,nl),d_deltaqw(nlon,nl) 1772 REAL sigmaw(nlon),d_sigmaw(nlon) 1773 LOGICAL wk_adv(nlon) 1774 INTEGER nl,nlon 1775 c Output 1776 REAL alpha(nlon) 1777 c Internal variables 1778 REAL zeta(nlon,nl) 1779 REAL alpha1(nlon) 1780 REAL x,a,b,c,discrim 1781 REAL epsilon 1782 ! DATA epsilon/1.e-15/ 1783 c 1784 DO k=1,nl 1785 DO i = 1,nlon 1786 IF (wk_adv(i)) THEN 1787 IF ((deltaqw(i,k)+d_deltaqw(i,k)).ge.0.) then 1788 zeta(i,k)=0. 1789 ELSE 1790 zeta(i,k)=1. 1791 END IF 1792 ENDIF 1793 ENDDO 1794 DO i = 1,nlon 1795 IF (wk_adv(i)) THEN 1796 x = qe(i,k)+(zeta(i,k)-sigmaw(i))*deltaqw(i,k) 1797 $ + d_qe(i,k)+(zeta(i,k)-sigmaw(i))*d_deltaqw(i,k) 1798 $ - d_sigmaw(i)*(deltaqw(i,k)+d_deltaqw(i,k)) 1799 a = -d_sigmaw(i)*d_deltaqw(i,k) 1800 b = d_qe(i,k)+(zeta(i,k)-sigmaw(i))*d_deltaqw(i,k) 1801 $ - deltaqw(i,k)*d_sigmaw(i) 1802 c = qe(i,k)+(zeta(i,k)-sigmaw(i))*deltaqw(i,k)+epsilon 1803 discrim = b*b-4.*a*c 1804 c print*, 'x, a, b, c, discrim', x, a, b, c, discrim 1805 IF (a+b .GE. 0.) THEN !! Condition suffisante pour la positivité de ovap 1806 alpha1(i)=1. 1807 ELSE 1808 IF (x .GE. 0.) THEN 1809 alpha1(i)=1. 1810 ELSE 1811 IF (a .GT. 0.) THEN 1812 alpha1(i)=0.9*min( (2.*c)/(-b+sqrt(discrim)), 1813 $ (-b+sqrt(discrim))/(2.*a) ) 1814 ELSE IF (a .eq. 0.) then 1815 alpha1(i)=0.9*(-c/b) 1816 ELSE 1817 c print*,'a,b,c discrim',a,b,c discrim 1818 alpha1(i)=0.9*max( (2.*c)/(-b+sqrt(discrim)), 1819 $ (-b+sqrt(discrim))/(2.*a) ) 1820 ENDIF 1821 ENDIF 1822 ENDIF 1823 alpha(i) = min(alpha(i),alpha1(i)) 1824 ENDIF 1825 ENDDO 1826 ENDDO 1827 ! 1828 return 1829 end 1830 1831 Subroutine WAKE_scal (p,ph,ppi,dtime,sigd_con 1832 : ,te0,qe0,omgb 1833 : ,dtdwn,dqdwn,amdwn,amup,dta,dqa 1834 : ,wdtPBL,wdqPBL,udtPBL,udqPBL 1835 o ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl 1836 o ,dtls,dqls 1837 o ,ktopw,omgbdth,dp_omgb,wdens 1838 o ,tu,qu 1839 o ,dtKE,dqKE 1840 o ,dtPBL,dqPBL 1841 o ,omg,dp_deltomg,spread 1842 o ,Cstar,d_deltat_gw 1843 o ,d_deltatw2,d_deltaqw2) 1844 1845 *************************************************************** 1846 * * 1847 * WAKE * 1848 * retour a un Pupper fixe * 1849 * * 1850 * written by : GRANDPEIX Jean-Yves 09/03/2000 * 1851 * modified by : ROEHRIG Romain 01/29/2007 * 1852 *************************************************************** 1853 c 1854 USE dimphy 26 1855 IMPLICIT none 27 1856 c============================================================================ … … 113 1942 114 1943 #include "dimensions.h" 115 #include "YOMCST.h"116 #include "cvthermo.h"117 #include "iniprint.h"118 119 c Arguments en entree120 c--------------------121 122 REAL, dimension(klon,klev) :: p, ppi123 REAL, dimension(klon,klev+1) :: ph, omgb124 REAL dtime125 REAL, dimension(klon,klev) :: te0,qe0126 REAL, dimension(klon,klev) :: dtdwn, dqdwn127 REAL, dimension(klon,klev) :: wdtPBL,wdqPBL128 REAL, dimension(klon,klev) :: udtPBL,udqPBL129 REAL, dimension(klon,klev) :: amdwn, amup130 REAL, dimension(klon,klev) :: dta, dqa131 REAL, dimension(klon) :: sigd_con132 133 c Sorties134 c--------135 136 REAL, dimension(klon,klev) :: deltatw, deltaqw, dth137 REAL, dimension(klon,klev) :: tu, qu138 REAL, dimension(klon,klev) :: dtls, dqls139 REAL, dimension(klon,klev) :: dtKE, dqKE140 REAL, dimension(klon,klev) :: dtPBL, dqPBL141 REAL, dimension(klon,klev) :: spread142 REAL, dimension(klon,klev) :: d_deltatgw143 REAL, dimension(klon,klev) :: d_deltatw2, d_deltaqw2144 REAL, dimension(klon,klev+1) :: omgbdth, omg145 REAL, dimension(klon,klev) :: dp_omgb, dp_deltomg146 REAL, dimension(klon,klev) :: d_deltat_gw147 REAL, dimension(klon) :: hw, sigmaw, wape, fip, gfl, Cstar148 INTEGER, dimension(klon) :: ktopw149 150 c Variables internes151 c-------------------152 153 c Variables à fixer154 REAL ALON155 REAL coefgw156 REAL :: wdens0, wdens157 REAL stark158 REAL alpk159 REAL delta_t_min160 INTEGER nsub161 REAL dtimesub162 REAL sigmad, hwmin163 REAL :: sigmaw_max164 cIM 080208165 LOGICAL, dimension(klon) :: gwake166 167 c Variables de sauvegarde168 REAL, dimension(klon,klev) :: deltatw0169 REAL, dimension(klon,klev) :: deltaqw0170 REAL, dimension(klon,klev) :: te, qe171 REAL, dimension(klon) :: sigmaw0, sigmaw1172 173 c Variables pour les GW174 REAL, DIMENSION(klon) :: LL175 REAL, dimension(klon,klev) :: N2176 REAL, dimension(klon,klev) :: Cgw177 REAL, dimension(klon,klev) :: Tgw178 179 c Variables liées au calcul de hw180 REAL, DIMENSION(klon) :: ptop_provis, ptop, ptop_new181 REAL, DIMENSION(klon) :: sum_dth182 REAL, DIMENSION(klon) :: dthmin183 REAL, DIMENSION(klon) :: z, dz, hw0184 INTEGER, DIMENSION(klon) :: ktop, kupper185 186 c Sub-timestep tendencies and related variables187 REAL d_deltatw(klon,klev),d_deltaqw(klon,klev)188 REAL d_te(klon,klev),d_qe(klon,klev)189 REAL d_sigmaw(klon),alpha(klon)190 REAL q0_min(klon),q1_min(klon)191 LOGICAL wk_adv(klon), OK_qx_qw(klon)192 193 c Autres variables internes194 INTEGER isubstep, k, i195 196 REAL, DIMENSION(klon) :: sum_thu, sum_tu, sum_qu,sum_thvu197 REAL, DIMENSION(klon) :: sum_dq, sum_rho198 REAL, DIMENSION(klon) :: sum_dtdwn, sum_dqdwn199 REAL, DIMENSION(klon) :: av_thu, av_tu, av_qu, av_thvu200 REAL, DIMENSION(klon) :: av_dth, av_dq, av_rho201 REAL, DIMENSION(klon) :: av_dtdwn, av_dqdwn202 203 REAL, DIMENSION(klon,klev) :: rho, rhow204 REAL, DIMENSION(klon,klev+1) :: rhoh205 REAL, DIMENSION(klon,klev) :: rhow_moyen206 REAL, DIMENSION(klon,klev) :: zh207 REAL, DIMENSION(klon,klev+1) :: zhh208 REAL, DIMENSION(klon,klev) :: epaisseur1, epaisseur2209 210 REAL, DIMENSION(klon,klev) :: the, thu211 212 ! REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw213 214 REAL, DIMENSION(klon,klev+1) :: omgbw215 REAL, DIMENSION(klon) :: pupper216 REAL, DIMENSION(klon) :: omgtop217 REAL, DIMENSION(klon,klev) :: dp_omgbw218 REAL, DIMENSION(klon) :: ztop, dztop219 REAL, DIMENSION(klon,klev) :: alpha_up220 221 REAL, dimension(klon) :: RRe1, RRe2222 REAL :: RRd1, RRd2223 REAL, DIMENSION(klon,klev) :: Th1, Th2, q1, q2, T1224 REAL, DIMENSION(klon,klev) :: D_Th1, D_Th2, D_dth225 REAL, DIMENSION(klon,klev) :: D_q1, D_q2, D_dq226 REAL, DIMENSION(klon,klev) :: omgbdq227 228 REAL, dimension(klon) :: ff, gg229 REAL, dimension(klon) :: wape2, Cstar2, heff230 231 REAL, DIMENSION(klon,klev) :: Crep232 REAL Crep_upper, Crep_sol233 234 C-------------------------------------------------------------------------235 c Initialisations236 c-------------------------------------------------------------------------237 238 c print*, 'wake initialisations'239 240 c Essais d'initialisation avec sigmaw = 0.02 et hw = 10.241 c-------------------------------------------------------------------------242 243 DATA sigmad, hwmin /.02,10./244 245 C Longueur de maille (en m)246 c-------------------------------------------------------------------------247 248 c ALON = 3.e5249 ALON = 1.e6250 251 252 C Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)253 c254 c coefgw : Coefficient pour les ondes de gravité255 c stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)256 c wdens : Densité de poche froide par maille257 c-------------------------------------------------------------------------258 259 coefgw=10260 c coefgw=1261 c wdens0 = 1.0/(alon**2)262 wdens = 1.0/(alon**2)263 stark = 0.50264 cCRtest265 alpk=0.1266 c alpk = 1.0267 c alpk = 0.5268 c alpk = 0.05269 Crep_upper=0.9270 Crep_sol=1.0271 272 C Minimum value for |T_wake - T_undist|. Used for wake top definition273 c-------------------------------------------------------------------------274 275 delta_t_min = 0.2276 277 C 1. - Save initial values and initialize tendencies278 C --------------------------------------------------279 280 DO k=1,klev281 DO i=1, klon282 deltatw0(i,k) = deltatw(i,k)283 deltaqw0(i,k)= deltaqw(i,k)284 te(i,k) = te0(i,k)285 qe(i,k) = qe0(i,k)286 dtls(i,k) = 0.287 dqls(i,k) = 0.288 d_deltat_gw(i,k)=0.289 d_te(i,k) = 0.290 d_qe(i,k) = 0.291 d_deltatw(i,k) = 0.292 d_deltaqw(i,k) = 0.293 !IM 060508 beg294 d_deltatw2(i,k)=0.295 d_deltaqw2(i,k)=0.296 !IM 060508 end297 ENDDO298 ENDDO299 c sigmaw1=sigmaw300 c IF (sigd_con.GT.sigmaw1) THEN301 c print*, 'sigmaw,sigd_con', sigmaw, sigd_con302 c ENDIF303 DO i=1, klon304 cc sigmaw(i) = amax1(sigmaw(i),sigd_con(i))305 sigmaw(i) = amax1(sigmaw(i),sigmad)306 sigmaw(i) = amin1(sigmaw(i),0.99)307 sigmaw0(i) = sigmaw(i)308 wape(i) = 0.309 wape2(i) = 0.310 d_sigmaw(i) = 0.311 ktopw(i) = 0312 ENDDO313 C314 C315 C 2. - Prognostic part316 C --------------------317 C318 C319 C 2.1 - Undisturbed area and Wake integrals320 C ---------------------------------------------------------321 322 DO i=1, klon323 z(i) = 0.324 ktop(i)=0325 kupper(i) = 0326 sum_thu(i) = 0.327 sum_tu(i) = 0.328 sum_qu(i) = 0.329 sum_thvu(i) = 0.330 sum_dth(i) = 0.331 sum_dq(i) = 0.332 sum_rho(i) = 0.333 sum_dtdwn(i) = 0.334 sum_dqdwn(i) = 0.335 336 av_thu(i) = 0.337 av_tu(i) =0.338 av_qu(i) =0.339 av_thvu(i) = 0.340 av_dth(i) = 0.341 av_dq(i) = 0.342 av_rho(i) =0.343 av_dtdwn(i) =0.344 av_dqdwn(i) = 0.345 ENDDO346 c347 c Distance between wakes348 DO i = 1,klon349 LL(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens)350 ENDDO351 C Potential temperatures and humidity352 c----------------------------------------------------------353 DO k =1,klev354 DO i=1, klon355 rho(i,k) = p(i,k)/(rd*te(i,k))356 IF(k .eq. 1) THEN357 rhoh(i,k) = ph(i,k)/(rd*te(i,k))358 zhh(i,k)=0359 ELSE360 rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1)))361 zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1)362 ENDIF363 the(i,k) = te(i,k)/ppi(i,k)364 thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k)365 tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i)366 qu(i,k) = qe(i,k) - deltaqw(i,k)*sigmaw(i)367 rhow(i,k) = p(i,k)/(rd*(tu(i,k)+deltatw(i,k)))368 dth(i,k) = deltatw(i,k)/ppi(i,k)369 ENDDO370 ENDDO371 372 DO k = 1, klev-1373 DO i=1, klon374 IF(k.eq.1) THEN375 N2(i,k)=0376 ELSE377 N2(i,k)=amax1(0.,-RG**2/the(i,k)*rho(i,k)*(the(i,k+1)-378 $ the(i,k-1))/(p(i,k+1)-p(i,k-1)))379 ENDIF380 ZH(i,k)=(zhh(i,k)+zhh(i,k+1))/2381 382 Cgw(i,k)=sqrt(N2(i,k))*ZH(i,k)383 Tgw(i,k)=coefgw*Cgw(i,k)/LL(i)384 ENDDO385 ENDDO386 387 DO i=1, klon388 N2(i,klev)=0389 ZH(i,klev)=0390 Cgw(i,klev)=0391 Tgw(i,klev)=0392 ENDDO393 394 c Calcul de la masse volumique moyenne de la colonne (bdlmd)395 c-----------------------------------------------------------------396 397 DO k=1,klev398 DO i=1, klon399 epaisseur1(i,k)=0.400 epaisseur2(i,k)=0.401 ENDDO402 ENDDO403 404 DO i=1, klon405 epaisseur1(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1.406 epaisseur2(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1.407 rhow_moyen(i,1) = rhow(i,1)408 ENDDO409 410 DO k = 2, klev411 DO i=1, klon412 epaisseur1(i,k)= -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) +1.413 epaisseur2(i,k)=epaisseur2(i,k-1)+epaisseur1(i,k)414 rhow_moyen(i,k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+415 $ rhow(i,k)*epaisseur1(i,k))/epaisseur2(i,k)416 ENDDO417 ENDDO418 419 C420 C Choose an integration bound well above wake top421 c-----------------------------------------------------------------422 c423 C Pupper = 50000. ! melting level424 c Pupper = 60000.425 c Pupper = 80000. ! essais pour case_e426 DO i = 1,klon427 ccc Pupper(i) = 0.6*ph(i,1)428 Pupper(i) = 60000.429 ENDDO430 431 C432 C Determine Wake top pressure (Ptop) from buoyancy integral433 C --------------------------------------------------------434 c435 c-1/ Pressure of the level where dth becomes less than delta_t_min.436 437 DO i=1,klon438 ptop_provis(i)=ph(i,1)439 ENDDO440 DO k= 2,klev441 DO i=1,klon442 c443 cIM v3JYG; ptop_provis(i).LT. ph(i,1)444 c445 IF (dth(i,k) .GT. -delta_t_min .and.446 $ dth(i,k-1).LT. -delta_t_min .and.447 $ ptop_provis(i).EQ. ph(i,1)) THEN448 ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)449 $ - (dth(i,k-1)+delta_t_min)*p(i,k)) /450 $ (dth(i,k) - dth(i,k-1))451 ENDIF452 ENDDO453 ENDDO454 455 c-2/ dth integral456 457 DO i=1,klon458 sum_dth(i) = 0.459 dthmin(i) = -delta_t_min460 z(i) = 0.461 ENDDO462 463 DO k = 1,klev464 DO i=1,klon465 dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg)466 IF (dz(i) .gt. 0) THEN467 z(i) = z(i)+dz(i)468 sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)469 dthmin(i) = amin1(dthmin(i),dth(i,k))470 ENDIF471 ENDDO472 ENDDO473 474 c-3/ height of triangle with area= sum_dth and base = dthmin475 476 DO i=1,klon477 hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)478 hw0(i) = amax1(hwmin,hw0(i))479 ENDDO480 481 c-4/ now, get Ptop482 483 DO i=1,klon484 z(i) = 0.485 ptop(i) = ph(i,1)486 ENDDO487 488 DO k = 1,klev489 DO i=1,klon490 dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg),hw0(i)-z(i))491 IF (dz(i) .gt. 0) THEN492 z(i) = z(i)+dz(i)493 ptop(i) = ph(i,k)-rho(i,k)*rg*dz(i)494 ENDIF495 ENDDO496 ENDDO497 498 499 C-5/ Determination de ktop et kupper500 501 DO k=klev,1,-1502 DO i=1,klon503 IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k504 IF (ph(i,k+1) .lt. pupper(i)) kupper(i)=k505 ENDDO506 ENDDO507 508 c-6/ Correct ktop and ptop509 510 DO i = 1,klon511 ptop_new(i)=ptop(i)512 ENDDO513 DO k= klev,2,-1514 DO i=1,klon515 IF (k .LE. ktop(i) .and.516 $ ptop_new(i) .EQ. ptop(i) .and.517 $ dth(i,k) .GT. -delta_t_min .and.518 $ dth(i,k-1).LT. -delta_t_min) THEN519 ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)520 $ - (dth(i,k-1)+delta_t_min)*p(i,k)) /521 $ (dth(i,k) - dth(i,k-1))522 ENDIF523 ENDDO524 ENDDO525 526 DO i=1,klon527 ptop(i) = ptop_new(i)528 ENDDO529 530 DO k=klev,1,-1531 DO i=1,klon532 IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k533 ENDDO534 ENDDO535 c536 c-5/ Set deltatw & deltaqw to 0 above kupper537 c538 DO k = 1,klev539 DO i=1,klon540 IF (k.GE. kupper(i)) THEN541 deltatw(i,k) = 0.542 deltaqw(i,k) = 0.543 ENDIF544 ENDDO545 ENDDO546 c547 C548 C Vertical gradient of LS omega549 C550 DO k = 1,klev551 DO i=1,klon552 IF (k.LE. kupper(i)) THEN553 dp_omgb(i,k) = (omgb(i,k+1) - omgb(i,k))/(ph(i,k+1)-ph(i,k))554 ENDIF555 ENDDO556 ENDDO557 C558 C Integrals (and wake top level number)559 C --------------------------------------560 C561 C Initialize sum_thvu to 1st level virt. pot. temp.562 563 DO i=1,klon564 z(i) = 1.565 dz(i) = 1.566 sum_thvu(i) = thu(i,1)*(1.+eps*qu(i,1))*dz(i)567 sum_dth(i) = 0.568 ENDDO569 570 DO k = 1,klev571 DO i=1,klon572 dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)573 IF (dz(i) .GT. 0) THEN574 z(i) = z(i)+dz(i)575 sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i)576 sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i)577 sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i)578 sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i)579 sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)580 sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)581 sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)582 sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)583 sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)584 ENDIF585 ENDDO586 ENDDO587 c588 DO i=1,klon589 hw0(i) = z(i)590 ENDDO591 c592 C593 C 2.1 - WAPE and mean forcing computation594 C ---------------------------------------595 C596 C ---------------------------------------597 C598 C Means599 600 DO i=1,klon601 av_thu(i) = sum_thu(i)/hw0(i)602 av_tu(i) = sum_tu(i)/hw0(i)603 av_qu(i) = sum_qu(i)/hw0(i)604 av_thvu(i) = sum_thvu(i)/hw0(i)605 c av_thve = sum_thve/hw0606 av_dth(i) = sum_dth(i)/hw0(i)607 av_dq(i) = sum_dq(i)/hw0(i)608 av_rho(i) = sum_rho(i)/hw0(i)609 av_dtdwn(i) = sum_dtdwn(i)/hw0(i)610 av_dqdwn(i) = sum_dqdwn(i)/hw0(i)611 612 wape(i) = - rg*hw0(i)*(av_dth(i)613 $ + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)*614 $ av_dq(i) ))/av_thvu(i)615 ENDDO616 C617 C 2.2 Prognostic variable update618 C ------------------------------619 C620 C Filter out bad wakes621 622 DO k = 1,klev623 DO i=1,klon624 IF ( wape(i) .LT. 0.) THEN625 deltatw(i,k) = 0.626 deltaqw(i,k) = 0.627 dth(i,k) = 0.628 ENDIF629 ENDDO630 ENDDO631 c632 DO i=1,klon633 IF ( wape(i) .LT. 0.) THEN634 wape(i) = 0.635 Cstar(i) = 0.636 hw(i) = hwmin637 sigmaw(i) = amax1(sigmad,sigd_con(i))638 fip(i) = 0.639 gwake(i) = .FALSE.640 ELSE641 Cstar(i) = stark*sqrt(2.*wape(i))642 gwake(i) = .TRUE.643 ENDIF644 ENDDO645 646 c647 c Check qx and qw positivity648 c --------------------------649 DO i = 1,klon650 q0_min(i)=min( (qe(i,1)-sigmaw(i)*deltaqw(i,1)),651 $ (qe(i,1)+(1.-sigmaw(i))*deltaqw(i,1)) )652 ENDDO653 DO k = 2,klev654 DO i = 1,klon655 q1_min(i)=min( (qe(i,k)-sigmaw(i)*deltaqw(i,k)),656 $ (qe(i,k)+(1.-sigmaw(i))*deltaqw(i,k)) )657 IF (q1_min(i).le.q0_min(i)) THEN658 q0_min(i)=q1_min(i)659 ENDIF660 ENDDO661 ENDDO662 c663 DO i = 1,klon664 OK_qx_qw(i) = q0_min(i) .GE. 0.665 alpha(i) = 1.666 ENDDO667 c668 CC -----------------------------------------------------------------669 C Sub-time-stepping670 C -----------------671 C672 nsub=10673 dtimesub=dtime/nsub674 c675 c------------------------------------------------------------676 DO isubstep = 1,nsub677 c------------------------------------------------------------678 c679 c wk_adv is the logical flag enabling wake evolution in the time advance loop680 DO i = 1,klon681 wk_adv(i) = OK_qx_qw(i) .AND. alpha(i) .GE. 1.682 ENDDO683 c684 DO i=1,klon685 IF (wk_adv(i)) THEN686 gfl(i) = 2.*sqrt(3.14*wdens*sigmaw(i))687 ENDIF688 ENDDO689 DO i=1,klon690 IF (wk_adv(i)) THEN691 d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub692 c sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub693 c sigmaw(i) =min(sigmaw(i),0.99) !!!!!!!!694 c wdens = wdens0/(10.*sigmaw)695 c sigmaw =max(sigmaw,sigd_con)696 c sigmaw =max(sigmaw,sigmad)697 ENDIF698 ENDDO699 C700 C701 c calcul de la difference de vitesse verticale poche - zone non perturbee702 cIM 060208 differences par rapport au code initial; init. a 0 dp_deltomg703 cIM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit704 cIM 060208 au niveau k=1..?705 DO k= 1,klev706 DO i = 1,klon707 dp_deltomg(i,k)=0.708 ENDDO709 ENDDO710 DO k= 1,klev+1711 DO i = 1,klon712 omg(i,k)=0.713 ENDDO714 ENDDO715 c716 DO i=1,klon717 IF (wk_adv(i)) THEN718 z(i)= 0.719 omg(i,1) = 0.720 dp_deltomg(i,1) = -(gfl(i)*Cstar(i))/(sigmaw(i) * (1-sigmaw(i)))721 ENDIF722 ENDDO723 c724 DO k= 2,klev725 DO i = 1,klon726 IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN727 dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)728 z(i) = z(i)+dz(i)729 dp_deltomg(i,k)= dp_deltomg(i,1)730 omg(i,k)= dp_deltomg(i,1)*z(i)731 ENDIF732 ENDDO733 ENDDO734 c735 DO i = 1,klon736 IF (wk_adv(i)) THEN737 dztop(i)=-(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)738 ztop(i) = z(i)+dztop(i)739 omgtop(i)=dp_deltomg(i,1)*ztop(i)740 ENDIF741 ENDDO742 c743 c -----------------744 c From m/s to Pa/s745 c -----------------746 c747 DO i=1,klon748 IF (wk_adv(i)) THEN749 omgtop(i) = -rho(i,ktop(i))*rg*omgtop(i)750 dp_deltomg(i,1) = omgtop(i)/(ptop(i)-ph(i,1))751 ENDIF752 ENDDO753 c754 DO k= 1,klev755 DO i = 1,klon756 IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN757 omg(i,k) = - rho(i,k)*rg*omg(i,k)758 dp_deltomg(i,k) = dp_deltomg(i,1)759 ENDIF760 ENDDO761 ENDDO762 c763 c raccordement lineaire de omg de ptop a pupper764 765 DO i=1,klon766 IF ( wk_adv(i) .AND. kupper(i) .GT. ktop(i)) THEN767 omg(i,kupper(i)+1) = - Rg*amdwn(i,kupper(i)+1)/sigmaw(i)768 $ + Rg*amup(i,kupper(i)+1)/(1.-sigmaw(i))769 dp_deltomg(i,kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/770 $ (ptop(i)-pupper(i))771 ENDIF772 ENDDO773 c774 DO k= 1,klev775 DO i = 1,klon776 IF( wk_adv(i) .AND. k .GT. ktop(i) .AND. k .LE. kupper(i)) THEN777 dp_deltomg(i,k) = dp_deltomg(i,kupper(i))778 omg(i,k) = omgtop(i)+(ph(i,k)-ptop(i))*dp_deltomg(i,kupper(i))779 ENDIF780 ENDDO781 ENDDO782 c783 c784 c-- Compute wake average vertical velocity omgbw785 c786 c787 DO k = 1,klev+1788 DO i=1,klon789 IF ( wk_adv(i)) THEN790 omgbw(i,k) = omgb(i,k)+(1.-sigmaw(i))*omg(i,k)791 ENDIF792 ENDDO793 ENDDO794 c-- and its vertical gradient dp_omgbw795 c796 DO k = 1,klev797 DO i=1,klon798 IF ( wk_adv(i)) THEN799 dp_omgbw(i,k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))800 ENDIF801 ENDDO802 ENDDO803 C804 c-- Upstream coefficients for omgb velocity805 c-- (alpha_up(k) is the coefficient of the value at level k)806 c-- (1-alpha_up(k) is the coefficient of the value at level k-1)807 DO k = 1,klev808 DO i=1,klon809 IF ( wk_adv(i)) THEN810 alpha_up(i,k) = 0.811 IF (omgb(i,k) .GT. 0.) alpha_up(i,k) = 1.812 ENDIF813 ENDDO814 ENDDO815 816 c Matrix expressing [The,deltatw] from [Th1,Th2]817 818 DO i=1,klon819 IF ( wk_adv(i)) THEN820 RRe1(i) = 1.-sigmaw(i)821 RRe2(i) = sigmaw(i)822 ENDIF823 ENDDO824 RRd1 = -1.825 RRd2 = 1.826 c827 c-- Get [Th1,Th2], dth and [q1,q2]828 c829 DO k= 1,klev830 DO i = 1,klon831 IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN832 dth(i,k) = deltatw(i,k)/ppi(i,k)833 Th1(i,k) = the(i,k) - sigmaw(i) *dth(i,k) ! undisturbed area834 Th2(i,k) = the(i,k) + (1.-sigmaw(i))*dth(i,k) ! wake835 q1(i,k) = qe(i,k) - sigmaw(i) *deltaqw(i,k) ! undisturbed area836 q2(i,k) = qe(i,k) + (1.-sigmaw(i))*deltaqw(i,k) ! wake837 T1(i,k) = te(i,k) - sigmaw(i)*deltatw(i,k)! undisturb itlmd838 ENDIF839 ENDDO840 ENDDO841 842 DO i=1,klon843 D_Th1(i,1) = 0. !!!itlmd : ne pas mettre if wk_adv cf nrlmd?844 D_Th2(i,1) = 0.845 D_dth(i,1) = 0.846 D_q1(i,1) = 0.847 D_q2(i,1) = 0.848 D_dq(i,1) = 0.849 ENDDO850 851 DO k= 2,klev852 DO i = 1,klon853 IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN854 D_Th1(i,k) = Th1(i,k-1)-Th1(i,k)855 D_Th2(i,k) = Th2(i,k-1)-Th2(i,k)856 D_dth(i,k) = dth(i,k-1)-dth(i,k)857 D_q1(i,k) = q1(i,k-1)-q1(i,k)858 D_q2(i,k) = q2(i,k-1)-q2(i,k)859 D_dq(i,k) = deltaqw(i,k-1)-deltaqw(i,k)860 ENDIF861 ENDDO862 ENDDO863 864 DO i=1,klon865 IF( wk_adv(i)) THEN866 omgbdth(i,1) = 0.867 omgbdq(i,1) = 0.868 ENDIF869 ENDDO870 871 DO k= 2,klev872 DO i = 1,klon873 IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN ! loop on interfaces874 omgbdth(i,k) = omgb(i,k)*( dth(i,k-1) - dth(i,k))875 omgbdq(i,k) = omgb(i,k)*(deltaqw(i,k-1) - deltaqw(i,k))876 ENDIF877 ENDDO878 ENDDO879 c880 c-----------------------------------------------------------------881 DO k= 1,klev882 DO i = 1,klon883 IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN884 c-----------------------------------------------------------------885 c886 c Compute redistribution (advective) term887 c888 d_deltatw(i,k) =889 $ dtimesub/(Ph(i,k)-Ph(i,k+1))*(890 $ RRd1*omg(i,k )*sigmaw(i) *D_Th1(i,k)891 $ -RRd2*omg(i,k+1)*(1.-sigmaw(i))*D_Th2(i,k+1)892 $ -(1.-alpha_up(i,k))*omgbdth(i,k) - alpha_up(i,k+1)*893 $ omgbdth(i,k+1))*ppi(i,k)894 c print*,'d_deltatw=',d_deltatw(i,k)895 c896 d_deltaqw(i,k) =897 $ dtimesub/(Ph(i,k)-Ph(i,k+1))*(898 $ RRd1*omg(i,k )*sigmaw(i) *D_q1(i,k)899 $ -RRd2*omg(i,k+1)*(1.-sigmaw(i))*D_q2(i,k+1)900 $ -(1.-alpha_up(i,k))*omgbdq(i,k) - alpha_up(i,k+1)*901 $ omgbdq(i,k+1))902 c print*,'d_deltaqw=',d_deltaqw(i,k)903 c904 c and increment large scale tendencies905 c906 907 c908 C909 CC -----------------------------------------------------------------910 d_te(i,k) = dtimesub*(911 $ ( RRe1(i)*omg(i,k )*sigmaw(i) *D_Th1(i,k)912 $ -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_Th2(i,k+1) )913 $ /(Ph(i,k)-Ph(i,k+1))914 $ -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*915 $ (omg(i,k)-omg(i,k+1))/(Ph(i,k)-Ph(i,k+1)) !instead of dp_deltomg(i,k)916 $ )*ppi(i,k)917 c918 d_qe(i,k) = dtimesub*(919 $ ( RRe1(i)*omg(i,k )*sigmaw(i) *D_q1(i,k)920 $ -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_q2(i,k+1) )921 $ /(Ph(i,k)-Ph(i,k+1))922 $ -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*923 $ (omg(i,k)-omg(i,k+1))/(Ph(i,k)-Ph(i,k+1))!instead of dp_deltomg(i,k)924 $ )925 ELSE IF(wk_adv(i) .AND. k .EQ. kupper(i)) THEN ! corr pour conserver l'eau926 927 d_te(i,k) = dtimesub*(928 $ ( RRe1(i)*omg(i,k )*sigmaw(i) *D_Th1(i,k))929 $ /(Ph(i,k)-Ph(i,k+1))930 $ )*ppi(i,k)931 932 d_qe(i,k) = dtimesub*(933 $ ( RRe1(i)*omg(i,k )*sigmaw(i) *D_q1(i,k))934 $ /(Ph(i,k)-Ph(i,k+1))935 $ )936 ENDIF937 938 c-------------------------------------------------------------------939 ENDDO940 ENDDO941 c------------------------------------------------------------------942 C943 C Increment state variables944 945 DO k= 1,klev946 DO i = 1,klon947 IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN948 c949 c Coefficient de répartition950 951 Crep(i,k)=Crep_sol*(ph(i,kupper(i))-ph(i,k))/(ph(i,kupper(i))952 $ -ph(i,1))953 Crep(i,k)=Crep(i,k)+Crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)-954 $ ph(i,kupper(i)))955 956 957 c Reintroduce compensating subsidence term.958 959 c dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw960 c dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))961 c . /(1-sigmaw)962 c dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw963 c dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))964 c . /(1-sigmaw)965 c966 c dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw967 c dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k))968 c . /(1-sigmaw)969 c dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw970 c dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k))971 c . /(1-sigmaw)972 973 dtKE(i,k)=(dtdwn(i,k)/sigmaw(i) - dta(i,k)/(1.-sigmaw(i)))974 dqKE(i,k)=(dqdwn(i,k)/sigmaw(i) - dqa(i,k)/(1.-sigmaw(i)))975 c print*,'dtKE=',dtKE(k)976 c print*,'dqKE=',dqKE(k)977 c978 dtPBL(i,k)=(wdtPBL(i,k)/sigmaw(i) - udtPBL(i,k)/(1.-sigmaw(i)))979 dqPBL(i,k)=(wdqPBL(i,k)/sigmaw(i) - udqPBL(i,k)/(1.-sigmaw(i)))980 c981 spread(i,k) = (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/982 $ sigmaw(i)983 984 985 c ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU Jingmei986 987 d_deltat_gw(i,k)=d_deltat_gw(i,k)-Tgw(i,k)*deltatw(i,k)*988 $ dtimesub989 ff(i)=d_deltatw(i,k)/dtimesub990 991 c Sans GW992 c993 c deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k))994 c995 c GW formule 1996 c997 c deltatw(k) = deltatw(k)+dtimesub*998 c $ (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))999 c1000 c GW formule 21001 1002 IF (dtimesub*Tgw(i,k).lt.1.e-10) THEN1003 d_deltatw(i,k) = dtimesub*1004 $ (ff(i)+dtKE(i,k)+dtPBL(i,k)1005 $ - spread(i,k)*deltatw(i,k)-Tgw(i,k)*deltatw(i,k))1006 ELSE1007 d_deltatw(i,k) = 1/Tgw(i,k)*(1-exp(-dtimesub*1008 $ Tgw(i,k)))*1009 $ (ff(i)+dtKE(i,k)+dtPBL(i,k)1010 $ - spread(i,k)*deltatw(i,k)-Tgw(i,k)*deltatw(i,k))1011 ENDIF1012 1013 dth(i,k) = deltatw(i,k)/ppi(i,k)1014 1015 gg(i)=d_deltaqw(i,k)/dtimesub1016 1017 d_deltaqw(i,k) = dtimesub*(gg(i)+ dqKE(i,k)+dqPBL(i,k)1018 $ - spread(i,k)*deltaqw(i,k))1019 1020 d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)1021 d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)1022 ENDIF1023 ENDDO1024 ENDDO1025 1026 C1027 C Scale tendencies so that water vapour remains positive in w and x.1028 C1029 call wake_vec_modulation(klon,klev,wk_adv,qe,d_qe,deltaqw,1030 $ d_deltaqw,sigmaw,d_sigmaw,alpha)1031 c1032 DO k = 1,klev1033 DO i = 1,klon1034 IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN1035 d_te(i,k)=alpha(i)*d_te(i,k)1036 d_qe(i,k)=alpha(i)*d_qe(i,k)1037 d_deltatw(i,k)=alpha(i)*d_deltatw(i,k)1038 d_deltaqw(i,k)=alpha(i)*d_deltaqw(i,k)1039 d_deltat_gw(i,k)=alpha(i)*d_deltat_gw(i,k)1040 ENDIF1041 ENDDO1042 ENDDO1043 DO i = 1,klon1044 IF( wk_adv(i)) THEN1045 d_sigmaw(i)=alpha(i)*d_sigmaw(i)1046 ENDIF1047 ENDDO1048 1049 C Update large scale variables and wake variables1050 cIM 060208 manque DO i + remplace DO k=1,kupper(i)1051 cIM 060208 DO k = 1,kupper(i)1052 DO k= 1,klev1053 DO i = 1,klon1054 IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN1055 dtls(i,k)=dtls(i,k)+d_te(i,k)1056 dqls(i,k)=dqls(i,k)+d_qe(i,k)1057 ENDIF1058 ENDDO1059 ENDDO1060 DO k= 1,klev1061 DO i = 1,klon1062 IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN1063 te(i,k) = te0(i,k) + dtls(i,k)1064 qe(i,k) = qe0(i,k) + dqls(i,k)1065 the(i,k) = te(i,k)/ppi(i,k)1066 deltatw(i,k) = deltatw(i,k)+d_deltatw(i,k)1067 deltaqw(i,k) = deltaqw(i,k)+d_deltaqw(i,k)1068 dth(i,k) = deltatw(i,k)/ppi(i,k)1069 ENDIF1070 ENDDO1071 ENDDO1072 DO i = 1,klon1073 IF( wk_adv(i)) THEN1074 sigmaw(i) = sigmaw(i)+d_sigmaw(i)1075 ENDIF1076 ENDDO1077 c1078 C1079 c Determine Ptop from buoyancy integral1080 c ---------------------------------------1081 c1082 c- 1/ Pressure of the level where dth changes sign.1083 c1084 DO i=1,klon1085 IF ( wk_adv(i)) THEN1086 Ptop_provis(i)=ph(i,1)1087 ENDIF1088 ENDDO1089 c1090 DO k= 2,klev1091 DO i=1,klon1092 IF ( wk_adv(i) .AND.1093 $ Ptop_provis(i) .EQ. ph(i,1) .AND.1094 $ dth(i,k) .GT. -delta_t_min .and.1095 $ dth(i,k-1).LT. -delta_t_min) THEN1096 Ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)1097 $ - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k)1098 $ - dth(i,k-1))1099 ENDIF1100 ENDDO1101 ENDDO1102 c1103 c- 2/ dth integral1104 c1105 DO i=1,klon1106 sum_dth(i) = 0.1107 dthmin(i) = -delta_t_min1108 z(i) = 0.1109 ENDDO1110 1111 DO k = 1,klev1112 DO i=1,klon1113 IF ( wk_adv(i)) THEN1114 dz(i) = -(amax1(ph(i,k+1),Ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg)1115 IF (dz(i) .gt. 0) THEN1116 z(i) = z(i)+dz(i)1117 sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)1118 dthmin(i) = amin1(dthmin(i),dth(i,k))1119 ENDIF1120 ENDIF1121 ENDDO1122 ENDDO1123 c1124 c- 3/ height of triangle with area= sum_dth and base = dthmin1125 1126 DO i=1,klon1127 IF ( wk_adv(i)) THEN1128 hw(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)1129 hw(i) = amax1(hwmin,hw(i))1130 ENDIF1131 ENDDO1132 c1133 c- 4/ now, get Ptop1134 c1135 DO i=1,klon1136 ktop(i) = 01137 z(i)=0.1138 ENDDO1139 c1140 DO k = 1,klev1141 DO i=1,klon1142 IF ( wk_adv(i)) THEN1143 dz(i) = amin1(-(ph(i,k+1)-Ph(i,k))/(rho(i,k)*rg),hw(i)-z(i))1144 IF (dz(i) .gt. 0) THEN1145 z(i) = z(i)+dz(i)1146 Ptop(i) = Ph(i,k)-rho(i,k)*rg*dz(i)1147 ktop(i) = k1148 ENDIF1149 ENDIF1150 ENDDO1151 ENDDO1152 c1153 c 4.5/Correct ktop and ptop1154 c1155 DO i=1,klon1156 IF ( wk_adv(i)) THEN1157 Ptop_new(i)=ptop(i)1158 ENDIF1159 ENDDO1160 c1161 DO k= klev,2,-11162 DO i=1,klon1163 cIM v3JYG; IF (k .GE. ktop(i)1164 IF ( wk_adv(i) .AND.1165 $ k .LE. ktop(i) .AND.1166 $ ptop_new(i) .EQ. ptop(i) .AND.1167 $ dth(i,k) .GT. -delta_t_min .and.1168 $ dth(i,k-1).LT. -delta_t_min) THEN1169 Ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)1170 $ - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k)1171 $ - dth(i,k-1))1172 ENDIF1173 ENDDO1174 ENDDO1175 c1176 c1177 DO i=1,klon1178 IF ( wk_adv(i)) THEN1179 ptop(i) = ptop_new(i)1180 ENDIF1181 ENDDO1182 1183 DO k=klev,1,-11184 DO i=1,klon1185 IF (ph(i,k+1) .LT. ptop(i)) ktop(i)=k1186 ENDDO1187 ENDDO1188 c1189 c 5/ Set deltatw & deltaqw to 0 above kupper1190 c1191 DO k = 1,klev1192 DO i=1,klon1193 IF ( wk_adv(i) .AND. k .GE. kupper(i)) THEN1194 deltatw(i,k) = 0.1195 deltaqw(i,k) = 0.1196 ENDIF1197 ENDDO1198 ENDDO1199 c1200 C1201 c-------------Cstar computation---------------------------------1202 DO i=1, klon1203 sum_thu(i) = 0.1204 sum_tu(i) = 0.1205 sum_qu(i) = 0.1206 sum_thvu(i) = 0.1207 sum_dth(i) = 0.1208 sum_dq(i) = 0.1209 sum_rho(i) = 0.1210 sum_dtdwn(i) = 0.1211 sum_dqdwn(i) = 0.1212 1213 av_thu(i) = 0.1214 av_tu(i) =0.1215 av_qu(i) =0.1216 av_thvu(i) = 0.1217 av_dth(i) = 0.1218 av_dq(i) = 0.1219 av_rho(i) =0.1220 av_dtdwn(i) =0.1221 av_dqdwn(i) = 0.1222 ENDDO1223 C1224 C Integrals (and wake top level number)1225 C --------------------------------------1226 C1227 C Initialize sum_thvu to 1st level virt. pot. temp.1228 1229 DO i=1,klon1230 z(i) = 1.1231 dz(i) = 1.1232 sum_thvu(i) = thu(i,1)*(1.+eps*qu(i,1))*dz(i)1233 sum_dth(i) = 0.1234 ENDDO1235 1236 DO k = 1,klev1237 DO i=1,klon1238 dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)1239 IF (dz(i) .GT. 0) THEN1240 z(i) = z(i)+dz(i)1241 sum_thu(i) = sum_thu(i) + th1(i,k)*dz(i)1242 sum_tu(i) = sum_tu(i) + t1(i,k)*dz(i)1243 sum_qu(i) = sum_qu(i) + q1(i,k)*dz(i)1244 sum_thvu(i) = sum_thvu(i) + th1(i,k)*(1.+eps*q1(i,k))*dz(i)!itlmd1245 1246 sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)1247 sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)1248 sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)1249 sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)1250 sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)1251 ENDIF1252 ENDDO1253 ENDDO1254 c1255 DO i=1,klon1256 hw0(i) = z(i)1257 ENDDO1258 c1259 C1260 C - WAPE and mean forcing computation1261 C ---------------------------------------1262 C1263 C ---------------------------------------1264 C1265 C Means1266 1267 DO i=1,klon1268 av_thu(i) = sum_thu(i)/hw0(i)1269 av_tu(i) = sum_tu(i)/hw0(i)1270 av_qu(i) = sum_qu(i)/hw0(i)1271 av_thvu(i) = sum_thvu(i)/hw0(i)1272 av_dth(i) = sum_dth(i)/hw0(i)1273 av_dq(i) = sum_dq(i)/hw0(i)1274 av_rho(i) = sum_rho(i)/hw0(i)1275 av_dtdwn(i) = sum_dtdwn(i)/hw0(i)1276 av_dqdwn(i) = sum_dqdwn(i)/hw0(i)1277 c1278 wape(i) = - rg*hw0(i)*(av_dth(i)1279 $ + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)*1280 $ av_dq(i) ))/av_thvu(i)1281 ENDDO1282 C1283 C Filter out bad wakes1284 1285 DO k = 1,klev1286 DO i=1,klon1287 IF ( wape(i) .LT. 0.) THEN1288 deltatw(i,k) = 0.1289 deltaqw(i,k) = 0.1290 dth(i,k) = 0.1291 ENDIF1292 ENDDO1293 ENDDO1294 c1295 DO i=1,klon1296 IF ( wape(i) .LT. 0.) THEN1297 wape(i) = 0.1298 Cstar(i) = 0.1299 hw(i) = hwmin1300 sigmaw(i) = max(sigmad,sigd_con(i))1301 fip(i) = 0.1302 gwake(i) = .FALSE.1303 ELSE1304 Cstar(i) = stark*sqrt(2.*wape(i))1305 gwake(i) = .TRUE.1306 ENDIF1307 ENDDO1308 1309 ENDDO ! end sub-timestep loop1310 C1311 C -----------------------------------------------------------------1312 c Get back to tendencies per second1313 c1314 DO k = 1,klev1315 DO i=1,klon1316 IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN !! corr conservation eau1317 dtls(i,k) = dtls(i,k)/dtime1318 dqls(i,k) = dqls(i,k)/dtime1319 d_deltatw2(i,k)=d_deltatw2(i,k)/dtime1320 d_deltaqw2(i,k)=d_deltaqw2(i,k)/dtime1321 d_deltat_gw(i,k) = d_deltat_gw(i,k)/dtime1322 ENDIF1323 ENDDO1324 ENDDO1325 c1326 c1327 c----------------------------------------------------------1328 c Determine wake final state; recompute wape, cstar, ktop;1329 c filter out bad wakes.1330 c----------------------------------------------------------1331 c1332 C 2.1 - Undisturbed area and Wake integrals1333 C ---------------------------------------------------------1334 1335 DO i=1,klon1336 z(i) = 0.1337 sum_thu(i) = 0.1338 sum_tu(i) = 0.1339 sum_qu(i) = 0.1340 sum_thvu(i) = 0.1341 sum_dth(i) = 0.1342 sum_dq(i) = 0.1343 sum_rho(i) = 0.1344 sum_dtdwn(i) = 0.1345 sum_dqdwn(i) = 0.1346 1347 av_thu(i) = 0.1348 av_tu(i) =0.1349 av_qu(i) =0.1350 av_thvu(i) = 0.1351 av_dth(i) = 0.1352 av_dq(i) = 0.1353 av_rho(i) =0.1354 av_dtdwn(i) =0.1355 av_dqdwn(i) = 0.1356 ENDDO1357 C Potential temperatures and humidity1358 c----------------------------------------------------------1359 1360 DO k =1,klev1361 DO i=1,klon1362 IF ( wk_adv(i)) THEN1363 rho(i,k) = p(i,k)/(rd*te(i,k))1364 IF(k .eq. 1) THEN1365 rhoh(i,k) = ph(i,k)/(rd*te(i,k))1366 zhh(i,k)=01367 ELSE1368 rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1)))1369 zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1)1370 ENDIF1371 the(i,k) = te(i,k)/ppi(i,k)1372 thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k)1373 tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i)1374 qu(i,k) = qe(i,k) - deltaqw(i,k)*sigmaw(i)1375 rhow(i,k) = p(i,k)/(rd*(tu(i,k)+deltatw(i,k)))1376 dth(i,k) = deltatw(i,k)/ppi(i,k)1377 ENDIF1378 ENDDO1379 ENDDO1380 1381 C Integrals (and wake top level number)1382 C -----------------------------------------------------------1383 1384 C Initialize sum_thvu to 1st level virt. pot. temp.1385 1386 DO i=1,klon1387 IF ( wk_adv(i)) THEN1388 z(i) = 1.1389 dz(i) = 1.1390 sum_thvu(i) = thu(i,1)*(1.+eps*qu(i,1))*dz(i)1391 sum_dth(i) = 0.1392 ENDIF1393 ENDDO1394 1395 DO k = 1,klev1396 DO i=1,klon1397 IF ( wk_adv(i)) THEN1398 dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)1399 IF (dz(i) .GT. 0) THEN1400 z(i) = z(i)+dz(i)1401 sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i)1402 sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i)1403 sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i)1404 sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i)1405 sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)1406 sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)1407 sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)1408 sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)1409 sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)1410 ENDIF1411 ENDIF1412 ENDDO1413 ENDDO1414 c1415 DO i=1,klon1416 IF ( wk_adv(i)) THEN1417 hw0(i) = z(i)1418 ENDIF1419 ENDDO1420 c1421 C - WAPE and mean forcing computation1422 C-------------------------------------------------------------1423 1424 C Means1425 1426 DO i=1, klon1427 IF ( wk_adv(i)) THEN1428 av_thu(i) = sum_thu(i)/hw0(i)1429 av_tu(i) = sum_tu(i)/hw0(i)1430 av_qu(i) = sum_qu(i)/hw0(i)1431 av_thvu(i) = sum_thvu(i)/hw0(i)1432 av_dth(i) = sum_dth(i)/hw0(i)1433 av_dq(i) = sum_dq(i)/hw0(i)1434 av_rho(i) = sum_rho(i)/hw0(i)1435 av_dtdwn(i) = sum_dtdwn(i)/hw0(i)1436 av_dqdwn(i) = sum_dqdwn(i)/hw0(i)1437 1438 wape2(i) = - rg*hw0(i)*(av_dth(i)1439 $ + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+1440 $ av_dth(i)*av_dq(i) ))/av_thvu(i)1441 ENDIF1442 ENDDO1443 1444 C Prognostic variable update1445 C ------------------------------------------------------------1446 1447 C Filter out bad wakes1448 c1449 DO k = 1,klev1450 DO i=1,klon1451 IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN1452 deltatw(i,k) = 0.1453 deltaqw(i,k) = 0.1454 dth(i,k) = 0.1455 ENDIF1456 ENDDO1457 ENDDO1458 c1459 1460 DO i=1, klon1461 IF ( wk_adv(i)) THEN1462 IF ( wape2(i) .LT. 0.) THEN1463 wape2(i) = 0.1464 Cstar2(i) = 0.1465 hw(i) = hwmin1466 sigmaw(i) = amax1(sigmad,sigd_con(i))1467 fip(i) = 0.1468 gwake(i) = .FALSE.1469 ELSE1470 if(prt_level.ge.10) print*,'wape2>0'1471 Cstar2(i) = stark*sqrt(2.*wape2(i))1472 gwake(i) = .TRUE.1473 ENDIF1474 ENDIF1475 ENDDO1476 c1477 DO i=1, klon1478 IF ( wk_adv(i)) THEN1479 ktopw(i) = ktop(i)1480 ENDIF1481 ENDDO1482 c1483 DO i=1, klon1484 IF ( wk_adv(i)) THEN1485 IF (ktopw(i) .gt. 0 .and. gwake(i)) then1486 1487 Cjyg1 Utilisation d'un h_efficace constant ( ~ feeding layer)1488 ccc heff = 600.1489 C Utilisation de la hauteur hw1490 cc heff = 0.7*hw1491 heff(i) = hw(i)1492 1493 FIP(i) = 0.5*rho(i,ktopw(i))*Cstar2(i)**3*heff(i)*2*1494 $ sqrt(sigmaw(i)*wdens*3.14)1495 FIP(i) = alpk * FIP(i)1496 Cjyg21497 ELSE1498 FIP(i) = 0.1499 ENDIF1500 ENDIF1501 ENDDO1502 c1503 C Limitation de sigmaw1504 c1505 C sécurité : si le wake occuppe plus de 90 % de la surface de la maille,1506 C alors il disparait en se mélangeant à la partie undisturbed1507 c1508 sigmaw_max = 0.91509 DO k = 1,klev1510 DO i=1, klon1511 c correction NICOLAS $ ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN1512 ! print*,'wape wape2 ktopw OK_qx_qw =',1513 ! $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i)1514 IF ((sigmaw(i).GT.sigmaw_max).or.1515 $ ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or.1516 $ (ktopw(i).le.2) .OR.1517 $ .not. OK_qx_qw(i)) THEN1518 cIM cf NR/JYG 251108 $ ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN1519 ccc IF (sigmaw(i).GT.0.9) THEN1520 dtls(i,k) = 0.1521 dqls(i,k) = 0.1522 deltatw(i,k) = 0.1523 deltaqw(i,k) = 0.1524 ENDIF1525 ENDDO1526 ENDDO1527 c1528 DO i=1, klon1529 IF ( (sigmaw(i).GT.sigmaw_max).or.1530 $ ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or.1531 $ (ktopw(i).le.2) .OR.1532 $ .not. OK_qx_qw(i)) THEN1533 ! correction NICOLAS $ ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN1534 ccc IF (sigmaw(i).GT.0.9) THEN1535 wape(i) = 0.1536 cstar(i)= 0. !!corr itlmd1537 hw(i) = hwmin1538 sigmaw(i) = sigmad1539 fip(i) = 0.1540 ELSE1541 wape(i) = wape2(i)1542 cstar(i)= cstar2(i) !!corr itlmd1543 ENDIF1544 ENDDO1545 c1546 c1547 RETURN1548 END1549 1550 SUBROUTINE wake_vec_modulation(nlon,nl,wk_adv,qe,d_qe,1551 $ deltaqw,d_deltaqw,sigmaw,d_sigmaw,alpha)1552 c------------------------------------------------------1553 cDtermination du coefficient alpha tel que les tendances1554 c corriges alpha*d_G, pour toutes les grandeurs G, correspondent1555 c a une humidite positive dans la zone (x) et dans la zone (w).1556 c------------------------------------------------------1557 c1558 1559 c Input1560 REAL qe(nlon,nl),d_qe(nlon,nl)1561 REAL deltaqw(nlon,nl),d_deltaqw(nlon,nl)1562 REAL sigmaw(nlon),d_sigmaw(nlon)1563 LOGICAL wk_adv(nlon)1564 INTEGER nl,nlon1565 c Output1566 REAL alpha(nlon)1567 c Internal variables1568 REAL alpha1(nlon)1569 REAL x,a,b,c,discrim,zeta(nlon)1570 REAL epsilon1571 DATA epsilon/1.e-15/1572 c1573 DO k=1,nl1574 DO i = 1,nlon1575 IF (wk_adv(i)) THEN1576 IF ((deltaqw(i,k)+d_deltaqw(i,k)).ge.0.) then1577 zeta(i)=0.1578 ELSE1579 zeta(i)=1.1580 END IF1581 ENDIF1582 ENDDO1583 DO i = 1,nlon1584 IF (wk_adv(i)) THEN1585 x = qe(i,k)+(zeta(i)-sigmaw(i))*deltaqw(i,k)1586 $ +d_qe(i,k)+(zeta(i)-sigmaw(i))*d_deltaqw(i,k)1587 $ -d_sigmaw(i)*(deltaqw(i,k)+d_deltaqw(i,k))1588 a=-d_sigmaw(i)*d_deltaqw(i,k)1589 b=d_qe(i,k)+(zeta(i)-sigmaw(i))*d_deltaqw(i,k)1590 $ -deltaqw(i,k)*d_sigmaw(i)1591 c=qe(i,k)+(zeta(i)-sigmaw(i))*deltaqw(i,k)-epsilon1592 ! c=qe(i,k)+(zeta(i)-sigmaw(i))*deltaqw(i,k)1593 1594 discrim=b*b-4.*a*c1595 ! print*,'ZETA *********************'1596 ! print*,'zeta sigmaw ',zeta(:)1597 ! print*,'SIGMA *********************'1598 ! print*,'sigmaw ',sigmaw(:)1599 1600 ! print*,' x ************************'1601 ! print*,'x ',x1602 ! print*,' a+b ************************'1603 ! print*,'a+b ',a+b1604 1605 ! print*,'a b c delta zeta ',a,b,c,discrim1606 IF (a+b .GE. 0.) THEN1607 alpha1(i)=1.1608 ELSE1609 IF (x .GE. 0.) THEN1610 alpha1(i)=1.1611 ELSE1612 ! IF (a .GE. 0.) THEN1613 IF (a .GT. 0.) THEN1614 ! print*,'a b c delta zeta ',a,b,c,discrim,zeta(i)1615 ! print*,'-b+sqrt(discrim) ',-b+sqrt(discrim)1616 alpha1(i)=0.9*min( (2.*c)/(-b+sqrt(discrim)),1617 $ (-b+sqrt(discrim))/(2.*a) )1618 ELSE IF (a.eq.0.) THEN1619 alpha1(i)=0.9*(-c/b)1620 ELSE1621 ! print*,'a b c delta zeta ',a,b,c,discrim,zeta(i)1622 ! print*,'-b+sqrt(discrim) ',-b+sqrt(discrim)1623 alpha1(i)=0.9*max( (2.*c)/(-b+sqrt(discrim)),1624 $ (-b+sqrt(discrim))/(2.*a) )1625 ENDIF1626 ENDIF1627 ENDIF1628 ENDIF1629 ENDDO1630 ENDDO1631 c1632 DO i = 1,nlon1633 IF (wk_adv(i)) THEN1634 alpha(i) = min(alpha(i),alpha1(i))1635 ENDIF1636 ENDDO1637 c1638 return1639 end1640 1641 Subroutine WAKE_scal (p,ph,ppi,dtime,sigd_con1642 : ,te0,qe0,omgb1643 : ,dtdwn,dqdwn,amdwn,amup,dta,dqa1644 : ,wdtPBL,wdqPBL,udtPBL,udqPBL1645 o ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl1646 o ,dtls,dqls1647 o ,ktopw,omgbdth,dp_omgb,wdens1648 o ,tu,qu1649 o ,dtKE,dqKE1650 o ,dtPBL,dqPBL1651 o ,omg,dp_deltomg,spread1652 o ,Cstar,d_deltat_gw1653 o ,d_deltatw2,d_deltaqw2)1654 1655 ***************************************************************1656 * *1657 * WAKE *1658 * retour a un Pupper fixe *1659 * *1660 * written by : GRANDPEIX Jean-Yves 09/03/2000 *1661 * modified by : ROEHRIG Romain 01/29/2007 *1662 ***************************************************************1663 c1664 USE dimphy1665 IMPLICIT none1666 c============================================================================1667 C1668 C1669 C But : Decrire le comportement des poches froides apparaissant dans les1670 C grands systemes convectifs, et fournir l'energie disponible pour1671 C le declenchement de nouvelles colonnes convectives.1672 C1673 C Variables d'etat : deltatw : ecart de temperature wake-undisturbed area1674 C deltaqw : ecart d'humidite wake-undisturbed area1675 C sigmaw : fraction d'aire occupee par la poche.1676 C1677 C Variable de sortie :1678 c1679 c wape : WAke Potential Energy1680 c fip : Front Incident Power (W/m2) - ALP1681 c gfl : Gust Front Length per unit area (m-1)1682 C dtls : large scale temperature tendency due to wake1683 C dqls : large scale humidity tendency due to wake1684 C hw : hauteur de la poche1685 C dp_omgb : vertical gradient of large scale omega1686 C omgbdth: flux of Delta_Theta transported by LS omega1687 C dtKE : differential heating (wake - unpertubed)1688 C dqKE : differential moistening (wake - unpertubed)1689 C omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)1690 C dp_deltomg : vertical gradient of omg (s-1)1691 C spread : spreading term in dt_wake and dq_wake1692 C deltatw : updated temperature difference (T_w-T_u).1693 C deltaqw : updated humidity difference (q_w-q_u).1694 C sigmaw : updated wake fractional area.1695 C d_deltat_gw : delta T tendency due to GW1696 c1697 C Variables d'entree :1698 c1699 c aire : aire de la maille1700 c te0 : temperature dans l'environnement (K)1701 C qe0 : humidite dans l'environnement (kg/kg)1702 C omgb : vitesse verticale moyenne sur la maille (Pa/s)1703 C dtdwn: source de chaleur due aux descentes (K/s)1704 C dqdwn: source d'humidite due aux descentes (kg/kg/s)1705 C dta : source de chaleur due courants satures et detrain (K/s)1706 C dqa : source d'humidite due aux courants satures et detra (kg/kg/s)1707 C amdwn: flux de masse total des descentes, par unite de1708 C surface de la maille (kg/m2/s)1709 C amup : flux de masse total des ascendances, par unite de1710 C surface de la maille (kg/m2/s)1711 C p : pressions aux milieux des couches (Pa)1712 C ph : pressions aux interfaces (Pa)1713 C ppi : (p/p_0)**kapa (adim)1714 C dtime: increment temporel (s)1715 c1716 C Variables internes :1717 c1718 c rhow : masse volumique de la poche froide1719 C rho : environment density at P levels1720 C rhoh : environment density at Ph levels1721 C te : environment temperature | may change within1722 C qe : environment humidity | sub-time-stepping1723 C the : environment potential temperature1724 C thu : potential temperature in undisturbed area1725 C tu : temperature in undisturbed area1726 C qu : humidity in undisturbed area1727 C dp_omgb: vertical gradient og LS omega1728 C omgbw : wake average vertical omega1729 C dp_omgbw: vertical gradient of omgbw1730 C omgbdq : flux of Delta_q transported by LS omega1731 C dth : potential temperature diff. wake-undist.1732 C th1 : first pot. temp. for vertical advection (=thu)1733 C th2 : second pot. temp. for vertical advection (=thw)1734 C q1 : first humidity for vertical advection1735 C q2 : second humidity for vertical advection1736 C d_deltatw : terme de redistribution pour deltatw1737 C d_deltaqw : terme de redistribution pour deltaqw1738 C deltatw0 : deltatw initial1739 C deltaqw0 : deltaqw initial1740 C hw0 : hw initial1741 C sigmaw0: sigmaw initial1742 C amflux : horizontal mass flux through wake boundary1743 C wdens : number of wakes per unit area (3D) or per1744 C unit length (2D)1745 C Tgw : 1 sur la période de onde de gravité1746 c Cgw : vitesse de propagation de onde de gravité1747 c LL : distance entre 2 poches1748 1749 c-------------------------------------------------------------------------1750 c Déclaration de variables1751 c-------------------------------------------------------------------------1752 1753 #include "dimensions.h"1754 1944 cccc#include "dimphy.h" 1755 1945 #include "YOMCST.h"
Note: See TracChangeset
for help on using the changeset viewer.
