Changeset 1943
- Timestamp:
- Jan 22, 2014, 10:51:36 AM (11 years ago)
- Location:
- LMDZ5/trunk/libf
- Files:
-
- 8 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/trunk/libf/phy1d/lmdz1d.F
r1914 r1943 357 357 c Le numero du jour est dans "day". L heure est traitee separement. 358 358 c La date complete est dans "daytime" (l'unite est le jour). 359 fnday=nday 359 if (nday>0) then 360 fnday=nday 361 else 362 fnday=-nday/float(day_step) 363 endif 364 360 365 c Special case for arm_cu which lasts less than one day : 53100s !! (MPL 20111026) 361 366 IF(forcing_type .EQ. 61) fnday=53100./86400. … … 370 375 call ymds2ju(annee_ref,mois,day_ref,heure,day) 371 376 day_ini = day 372 day_end = day_ini + nday377 day_end = day_ini + fnday 373 378 374 379 IF (forcing_type .eq.2) THEN … … 463 468 !! mpl et jyg le 22/08/2012 : 464 469 !! pour que les cas a flux de surface imposes marchent 465 IF(.NOT.ok_flux_surf ) THEN470 IF(.NOT.ok_flux_surf.or.max(abs(wtsurf),abs(wqsurf))>0.) THEN 466 471 fsens=-wtsurf*rcpd*rho(1) 467 472 flat=-wqsurf*rlvtt*rho(1) 468 473 print *,'Flux: ok_flux wtsurf wqsurf',ok_flux_surf,wtsurf,wqsurf 469 474 ENDIF 475 print*,'Flux sol ',fsens,flat 470 476 !! ok_flux_surf=.false. 471 477 !! fsens=-wtsurf*rcpd*rho(1) -
LMDZ5/trunk/libf/phylmd/calltherm.F90
r1907 r1943 174 174 do isplit=1,nsplit_thermals 175 175 176 if (iflag_thermals .eq.1) then177 CALL thermcell_2002(klon,klev,zdt &176 if (iflag_thermals>=1000) then 177 CALL thermcell_2002(klon,klev,zdt,iflag_thermals & 178 178 & ,pplay,paprs,pphi & 179 179 & ,u_seri,v_seri,t_seri,q_seri & 180 180 & ,d_u_the,d_v_the,d_t_the,d_q_the & 181 & ,zfm_therm,zentr_therm &181 & ,zfm_therm,zentr_therm,fraca,zw2 & 182 182 & ,r_aspect_thermals,30.,w2di_thermals & 183 183 & ,tau_thermals) … … 273 273 flag_bidouille_stratocu=iflag_thermals<=12.or.iflag_thermals==14.or.iflag_thermals==16.or.iflag_thermals==18 274 274 275 if (iflag_thermals<=12) then 276 lmax=1 275 ! Calcul a posteriori du niveau max des thermiques pour les schémas qui 276 ! ne la sortent pas. 277 if (iflag_thermals<=12.or.iflag_thermals>=1000) then 278 lmax(:)=1 277 279 do k=1,klev-1 278 280 zdetr_therm(:,k)=zentr_therm(:,k)+zfm_therm(:,k)-zfm_therm(:,k+1) 281 enddo 282 do k=1,klev-1 283 do i=1,klon 284 if (zfm_therm(i,k+1)>0.) lmax(i)=k 285 enddo 279 286 enddo 280 287 endif -
LMDZ5/trunk/libf/phylmd/hbtm.F
r1907 r1943 67 67 INTEGER isommet 68 68 cum PARAMETER (isommet=klev) ! limite max sommet pbl 69 REAL vk 70 PARAMETER (vk=0.35) ! Von Karman => passer a .41 ! cf U.Olgstrom 71 REAL ricr 72 PARAMETER (ricr=0.4) 73 REAL fak 74 PARAMETER (fak=8.5) ! b calcul du Prandtl et de dTetas 75 REAL fakn 76 PARAMETER (fakn=7.2) ! a 77 REAL onet 78 PARAMETER (onet=1.0/3.0) 79 REAL t_coup 80 PARAMETER(t_coup=273.15) 81 REAL zkmin 82 PARAMETER (zkmin=0.01) 83 REAL betam 84 PARAMETER (betam=15.0) ! pour Phim / h dans la S.L stable 85 REAL betah 86 PARAMETER (betah=15.0) 87 REAL betas 88 PARAMETER (betas=5.0) ! Phit dans la S.L. stable (mais 2 formes / z/OBL<>1 89 REAL sffrac 90 PARAMETER (sffrac=0.1) ! S.L. = z/h < .1 91 REAL binm 92 PARAMETER (binm=betam*sffrac) 93 REAL binh 94 PARAMETER (binh=betah*sffrac) 95 REAL ccon 96 PARAMETER (ccon=fak*sffrac*vk) 97 c 69 REAL, PARAMETER :: vk=0.35 ! Von Karman => passer a .41 ! cf U.Olgstrom 70 REAL, PARAMETER :: ricr=0.4 71 REAL, PARAMETER :: fak=8.5 ! b calcul du Prandtl et de dTetas 72 REAL, PARAMETER :: fakn=7.2 ! a 73 REAL, PARAMETER :: onet=1.0/3.0 74 REAL, PARAMETER:: t_coup=273.15 75 REAL, PARAMETER :: zkmin=0.01 76 REAL, PARAMETER :: betam=15.0 ! pour Phim / h dans la S.L stable 77 REAL, PARAMETER :: betah=15.0 78 REAL, PARAMETER :: betas=5.0 ! Phit dans la S.L. stable (mais 2 formes / z/OBL<>1 79 REAL, PARAMETER :: sffrac=0.1 ! S.L. = z/h < .1 80 REAL, PARAMETER :: usmin=1.e-12 81 REAL, PARAMETER :: binm=betam*sffrac 82 REAL, PARAMETER :: binh=betah*sffrac 83 REAL, PARAMETER :: ccon=fak*sffrac*vk 84 REAL, PARAMETER :: b1=70.,b2=20. 85 REAL, PARAMETER :: zref=2. ! Niveau de ref a 2m peut eventuellement 86 c etre choisi a 10m 98 87 REAL q_star,t_star 99 REAL b1,b2,b212,b2sr ! Lambert correlations T' q' avec T* q* 100 PARAMETER (b1=70.,b2=20.) 88 REAL b212,b2sr ! Lambert correlations T' q' avec T* q* 101 89 c 102 90 REAL z(klon,klev) 103 91 cAM REAL pcfm(klon,klev), pcfh(klon,klev) 104 cAM105 REAL zref106 PARAMETER (zref=2.) ! Niveau de ref a 2m peut eventuellement107 c etre choisi a 10m108 cMA109 c110 92 INTEGER i, k, j 111 93 REAL zxt … … 129 111 cAM REAL cgq(klon,2:klev) ! counter-gradient term for constituents 130 112 cAM REAL cgs(klon,2:klev) ! counter-gradient star (cg/flux) 131 REAL obklen(klon) ! Monin-Obukhov lengh113 REAL unsobklen(klon) ! Monin-Obukhov lengh 132 114 cAM REAL ztvd, ztvu, 133 115 REAL zdu2 … … 345 327 plcl(i) = 6000. 346 328 c Lambda = -u*^3 / (alpha.g.kvon.<w'Theta'v> 347 obklen(i) = -t(i,1)*ustar(i)**3/(RG*vk*heatv(i))329 unsobklen(i) = -RG*vk*heatv(i)/(t(i,1)*max(ustar(i),usmin)**3) 348 330 trmb1(i) = 0. 349 331 trmb2(i) = 0. … … 420 402 DO i = 1, knon 421 403 IF (check(i)) THEN 422 phiminv(i) = (1.-binm*pblh(i) /obklen(i))**onet404 phiminv(i) = (1.-binm*pblh(i)*unsobklen(i))**onet 423 405 c *************************************************** 424 406 c Wm ? et W* ? c'est la formule pour z/h < .1 … … 621 603 zxt=(Th_th(i)-zref*0.5*RG/RCPD/(1.+RVTMP2*qT_th(i))) 622 604 . *(1.+RETV*qT_th(i)) 623 phiminv(i) = (1. - binm*pblh(i) /obklen(i))**onet624 phihinv(i) = sqrt(1. - binh*pblh(i) /obklen(i))605 phiminv(i) = (1. - binm*pblh(i)*unsobklen(i))**onet 606 phihinv(i) = sqrt(1. - binh*pblh(i)*unsobklen(i)) 625 607 wm(i) = ustar(i)*phiminv(i) 626 608 fak2(i) = wm(i)*pblh(i)*vk … … 657 639 658 640 zh(i) = zmzp/pblh(i) 659 zl(i) = zmzp /obklen(i)641 zl(i) = zmzp*unsobklen(i) 660 642 zzh(i) = 0. 661 643 IF (zh(i).LE.1.0) zzh(i) = (1. - zh(i))**2 -
LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90
r1924 r1943 41 41 REAL, SAVE, ALLOCATABLE :: d_u_ajs(:,:), d_v_ajs(:,:) 42 42 !$OMP THREADPRIVATE(d_u_ajs, d_v_ajs) 43 REAL, SAVE, ALLOCATABLE :: d_u_ajsb(:,:), d_v_ajsb(:,:)44 !$OMP THREADPRIVATE(d_u_ajsb, d_v_ajsb)45 43 REAL, SAVE, ALLOCATABLE :: d_t_eva(:,:),d_q_eva(:,:) 46 44 !$OMP THREADPRIVATE(d_t_eva,d_q_eva) … … 303 301 allocate(d_t_ajs(klon,klev),d_q_ajs(klon,klev)) 304 302 allocate(d_u_ajs(klon,klev),d_v_ajs(klon,klev)) 305 allocate(d_u_ajsb(klon,klev),d_v_ajsb(klon,klev))306 303 allocate(d_t_eva(klon,klev),d_q_eva(klon,klev)) 307 304 allocate(d_t_lscst(klon,klev),d_q_lscst(klon,klev)) … … 456 453 deallocate(d_t_ajs,d_q_ajs) 457 454 deallocate(d_u_ajs,d_v_ajs) 458 deallocate(d_u_ajsb,d_v_ajsb)459 455 deallocate(d_t_eva,d_q_eva) 460 456 deallocate(d_t_lscst,d_q_lscst) -
LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90
r1938 r1943 206 206 rneb, rnebjn, zx_rh, frugs, agesno, d_t_dyn, d_q_dyn, & 207 207 d_u_dyn, d_v_dyn, d_t_con, d_t_ajsb, d_t_ajs, & 208 d_u_ajs b, d_u_ajs, d_v_ajsb, d_v_ajs, &208 d_u_ajs, d_v_ajs, & 209 209 d_u_con, d_v_con, d_q_con, d_q_ajs, d_t_lsc, & 210 210 d_t_eva, d_q_lsc, beta_prec, d_t_lscth, & … … 1046 1046 CALL histwrite_phy(o_dtthe, zx_tmp_fi3d) 1047 1047 IF (vars_defined) THEN 1048 zx_tmp_fi3d(1:klon,1:klev)=d_u_ajs(1:klon,1:klev)/pdtphys - & 1049 d_u_ajsb(1:klon,1:klev)/pdtphys 1048 zx_tmp_fi3d(1:klon,1:klev)=d_u_ajs(1:klon,1:klev)/pdtphys 1050 1049 ENDIF 1051 1050 CALL histwrite_phy(o_duthe, zx_tmp_fi3d) 1052 1051 IF (vars_defined) THEN 1053 zx_tmp_fi3d(1:klon,1:klev)=d_v_ajs(1:klon,1:klev)/pdtphys - & 1054 d_v_ajsb(1:klon,1:klev)/pdtphys 1052 zx_tmp_fi3d(1:klon,1:klev)=d_v_ajs(1:klon,1:klev)/pdtphys 1055 1053 ENDIF 1056 1054 CALL histwrite_phy(o_dvthe, zx_tmp_fi3d) -
LMDZ5/trunk/libf/phylmd/thermcell_dq.F90
r1907 r1943 53 53 cfl=max(cfl,fm(ig,k)/zzm) 54 54 if (entr(ig,k).gt.zzm) then 55 print*,'entr dt > m ',entr(ig,k)*ptimestep,masse(ig,k)56 abort_message = ' '55 print*,'entr*dt>m,1',k,entr(ig,k)*ptimestep,masse(ig,k) 56 abort_message = 'entr dt > m, 1st' 57 57 CALL abort_gcm (modname,abort_message,1) 58 58 endif … … 188 188 cfl=max(cfl,fm(ig,k)/zzm) 189 189 if (entr(ig,k).gt.zzm) then 190 print*,'entr dt > m ',entr(ig,k)*ptimestep,masse(ig,k)191 abort_message = ' '190 print*,'entr*dt>m,2',k,entr(ig,k)*ptimestep,masse(ig,k) 191 abort_message = 'entr dt > m, 2nd' 192 192 CALL abort_gcm (modname,abort_message,1) 193 193 endif -
LMDZ5/trunk/libf/phylmd/thermcell_main.F90
r1907 r1943 83 83 REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1) 84 84 real pphi(ngrid,nlay) 85 LOGICAL debut 85 86 86 87 ! local: … … 178 179 !FH/IM save f0 179 180 real zlevinter(klon) 180 logical debut181 181 real seuil 182 182 real csc(klon,klev) … … 256 256 ! 257 257 258 seuil=0.25 259 260 if (debut) then 261 ! call getin('dvdq',dvdq) 262 ! call getin('dqimpl',dqimpl) 263 264 if (iflag_thermals==15.or.iflag_thermals==16) then 265 dvdq=0 266 dqimpl=-1 267 else 268 dvdq=1 269 dqimpl=1 270 endif 271 272 fm0=0. 273 entr0=0. 274 detr0=0. 275 276 277 #undef wrgrads_thermcell 278 #ifdef wrgrads_thermcell 279 ! Initialisation des sorties grads pour les thermiques. 280 ! Pour l'instant en 1D sur le point igout. 281 ! Utilise par thermcell_out3d.h 282 str10='therm' 283 call inigrads(1,1,rlond(igout),1.,-180.,180.,jjm, & 284 & rlatd(igout),-90.,90.,1.,llm,pplay(igout,:),1., & 285 & ptimestep,str10,'therm ') 286 #endif 287 288 289 290 endif 291 292 fm=0. ; entr=0. ; detr=0. 293 294 295 icount=icount+1 258 seuil=0.25 259 260 if (debut) then 261 if (iflag_thermals==15.or.iflag_thermals==16) then 262 dvdq=0 263 dqimpl=-1 264 else 265 dvdq=1 266 dqimpl=1 267 endif 268 269 fm0=0. 270 entr0=0. 271 detr0=0. 272 endif 273 fm=0. ; entr=0. ; detr=0. 274 icount=icount+1 296 275 297 276 !IM 090508 beg -
LMDZ5/trunk/libf/phylmd/thermcell_old.F
r1907 r1943 1 SUBROUTINE thermcell_2002(ngrid,nlay,ptimestep 1 SUBROUTINE thermcell_2002(ngrid,nlay,ptimestep,iflag_thermals 2 2 s ,pplay,pplev,pphi 3 3 s ,pu,pv,pt,po 4 4 s ,pduadj,pdvadj,pdtadj,pdoadj 5 s ,fm0,entr0 6 c s ,pu_therm,pv_therm 5 s ,fm0,entr0,fraca,wa_moy 7 6 s ,r_aspect,l_mix,w2di,tho) 8 7 9 8 USE dimphy 9 USE write_field_phy 10 10 IMPLICIT NONE 11 11 … … 35 35 36 36 #include "dimensions.h" 37 cccc#include "dimphy.h"38 37 #include "YOMCST.h" 39 38 … … 41 40 c ---------- 42 41 43 INTEGER ngrid,nlay,w2di 42 INTEGER ngrid,nlay,w2di,iflag_thermals 44 43 REAL tho 45 44 real ptimestep,l_mix,r_aspect … … 50 49 REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1) 51 50 real pphi(ngrid,nlay) 52 53 integer idetr 54 save idetr 55 data idetr/3/ 56 c$OMP THREADPRIVATE(idetr) 51 real fraca(ngrid,nlay+1),zw2(ngrid,nlay+1) 52 53 integer,save :: idetr=3,lev_out=1 54 c$OMP THREADPRIVATE(idetr,lev_out) 57 55 58 56 c local: 59 57 c ------ 60 58 59 INTEGER, SAVE :: dvdq=0,flagdq=0,dqimpl=1 60 LOGICAL, SAVE :: debut=.true. 61 !$OMP THREADPRIVATE(dvdq,flagdq,debut,dqimpl) 62 61 63 INTEGER ig,k,l,lmax(klon,klev+1),lmaxa(klon),lmix(klon) 62 real zmax(klon),zw,zz,z w2(klon,klev+1),ztva(klon,klev),zzz64 real zmax(klon),zw,zz,ztva(klon,klev),zzz 63 65 64 66 real zlev(klon,klev+1),zlay(klon,klev) … … 79 81 real zha(klon,klev) 80 82 real wa_moy(klon,klev+1) 81 real fraca(klon,klev+1)82 83 real fracc(klon,klev+1) 83 84 real zf,zf2 … … 86 87 87 88 real count_time 88 integer isplit,nsplit,ialt89 parameter (nsplit=10)90 data isplit/0/91 save isplit92 c$OMP THREADPRIVATE(isplit)93 89 94 90 logical sorties … … 142 138 c --------------------------------------------------- 143 139 144 140 ! print*,'0 OK convect8' 145 141 146 142 DO 1010 l=1,nlay … … 155 151 1010 CONTINUE 156 152 157 cprint*,'1 OK convect8'153 ! print*,'1 OK convect8' 158 154 c -------------------- 159 155 c … … 177 173 c----------------------------------------------------------------------- 178 174 175 if (debut) then 176 flagdq=(iflag_thermals-1000)/100 177 dvdq=(iflag_thermals-(1000+flagdq*100))/10 178 if (flagdq==2) dqimpl=-1 179 if (flagdq==3) dqimpl=1 180 debut=.false. 181 endif 182 print*,'TH flag th ',iflag_thermals,flagdq,dvdq,dqimpl 183 179 184 do l=2,nlay 180 185 do ig=1,ngrid … … 192 197 enddo 193 198 194 cprint*,'2 OK convect8'199 ! print*,'2 OK convect8' 195 200 c----------------------------------------------------------------------- 196 201 c Calcul des densites … … 217 222 enddo 218 223 219 cprint*,'3 OK convect8'224 ! print*,'3 OK convect8' 220 225 c------------------------------------------------------------------ 221 226 c Calcul de w2, quarre de w a partir de la cape … … 274 279 enddo 275 280 276 cprint*,'4 OK convect8'281 ! print*,'4 OK convect8' 277 282 c Calcul de la couche correspondant a la hauteur du thermique 278 283 do k=1,nlay-1 … … 287 292 enddo 288 293 289 cprint*,'5 OK convect8'294 ! print*,'5 OK convect8' 290 295 c Calcule du w max du thermique 291 296 do k=1,nlay … … 315 320 enddo 316 321 317 cprint*,'6 OK convect8'322 ! print*,'6 OK convect8' 318 323 c Longueur caracteristique correspondant a la hauteur des thermiques. 319 324 do ig=1,ngrid … … 330 335 c call dump2d(iim,jjm-1,zmax(2:ngrid-1),'ZMAX ') 331 336 337 ! print*,'OKl336' 332 338 c Calcul de l'entrainement. 333 339 c Le rapport d'aspect relie la largeur de l'ascendance a l'epaisseur … … 348 354 enddo 349 355 350 c print*,'7 OK convect8' 356 357 ! print*,'7 OK convect8' 351 358 do k=1,klev+1 352 359 do ig=1,ngrid … … 359 366 enddo 360 367 361 cprint*,'8 OK convect8'368 ! print*,'8 OK convect8' 362 369 do ig=1,ngrid 363 370 lmaxa(ig)=1 … … 367 374 368 375 376 ! print*,'OKl372' 369 377 do l=1,nlay-2 370 378 do ig=1,ngrid … … 422 430 enddo 423 431 424 cprint*,'9 OK convect8'432 ! print*,'9 OK convect8' 425 433 c print*,'WA1 ',wa_moy 426 434 … … 433 441 c de la vitesse d'entrainement horizontal dans la couche alimentante. 434 442 443 ! print*,'OKl439' 435 444 do l=2,nlay 436 445 do ig=1,ngrid … … 463 472 enddo 464 473 465 cprint*,'10 OK convect8'474 ! print*,'10 OK convect8' 466 475 c print*,'WA2 ',wa_moy 467 476 c calcul de la fraction de la maille concernée par l'ascendance en tenant … … 500 509 enddo 501 510 502 cprint*,'11 OK convect8'511 ! print*,'11 OK convect8' 503 512 c print*,'Ea3 ',wa_moy 504 513 c------------------------------------------------------------------ … … 535 544 enddo 536 545 537 cprint*,'12 OK convect8'546 ! print*,'12 OK convect8' 538 547 c print*,'WA4 ',wa_moy 539 548 cc------------------------------------------------------------------ … … 589 598 590 599 4444 continue 600 ! print*,'OK 444 ' 591 601 592 602 if (w2di.eq.1) then … … 598 608 endif 599 609 600 if ( 1.eq.1) then610 if (flagdq==0) then 601 611 call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse 602 612 . ,zh,zdhadj,zha) 603 613 call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse 604 614 . ,zo,pdoadj,zoa) 605 else 615 print*,'THERMALS OPT 1' 616 else if (flagdq==1) then 606 617 call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca 607 618 . ,zh,zdhadj,zha) 608 619 call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca 609 620 . ,zo,pdoadj,zoa) 621 print*,'THERMALS OPT 2' 622 else 623 call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, 624 . zh,zdhadj,zha,lev_out) 625 call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, 626 . zo,pdoadj,zoa,lev_out) 627 print*,'THERMALS OPT 3',dqimpl 610 628 endif 611 629 612 if (1.eq.0) then 613 call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse 614 . ,fraca,zmax 615 . ,zu,zv,pduadj,pdvadj,zua,zva) 616 else 630 print*,'TH VENT ',dvdq 631 if (dvdq==0) then 632 ! print*,'TH VENT OK ',dvdq 617 633 call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse 618 634 . ,zu,pduadj,zua) 619 635 call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse 620 636 . ,zv,pdvadj,zva) 637 else if (dvdq==1) then 638 call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse 639 . ,fraca,zmax 640 . ,zu,zv,pduadj,pdvadj,zua,zva) 641 else if (dvdq==2) then 642 call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse 643 & ,fraca,zmax 644 & ,zu,zv,pduadj,pdvadj,zua,zva,lev_out) 645 else if (dvdq==3) then 646 call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse 647 & ,zu,pduadj,zua,lev_out) 648 call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse 649 & ,zv,pdvadj,zva,lev_out) 621 650 endif 651 652 ! CALL writefield_phy('duadj',pduadj,klev) 622 653 623 654 do l=1,nlay … … 632 663 633 664 634 cprint*,'13 OK convect8'665 ! print*,'13 OK convect8' 635 666 c print*,'WA5 ',wa_moy 636 667 do l=1,nlay … … 654 685 c enddo 655 686 656 cprint*,'14 OK convect8'687 ! print*,'14 OK convect8' 657 688 c------------------------------------------------------------------ 658 689 c Calculs pour les sorties … … 679 710 enddo 680 711 enddo 681 682 c print*,'15 OK convect8'683 684 isplit=isplit+1685 686 687 c #define und688 goto 123689 #ifdef und690 CALL writeg1d(1,nlay,wd,'wd ','wd ')691 CALL writeg1d(1,nlay,zwa,'wa ','wa ')692 CALL writeg1d(1,nlay,fracd,'fracd ','fracd ')693 CALL writeg1d(1,nlay,fraca,'fraca ','fraca ')694 CALL writeg1d(1,nlay,wa_moy,'wam ','wam ')695 CALL writeg1d(1,nlay,zla,'la ','la ')696 CALL writeg1d(1,nlay,zld,'ld ','ld ')697 CALL writeg1d(1,nlay,pt,'pt ','pt ')698 CALL writeg1d(1,nlay,zh,'zh ','zh ')699 CALL writeg1d(1,nlay,zha,'zha ','zha ')700 CALL writeg1d(1,nlay,zu,'zu ','zu ')701 CALL writeg1d(1,nlay,zv,'zv ','zv ')702 CALL writeg1d(1,nlay,zo,'zo ','zo ')703 CALL writeg1d(1,nlay,wh,'wh ','wh ')704 CALL writeg1d(1,nlay,wu,'wu ','wu ')705 CALL writeg1d(1,nlay,wv,'wv ','wv ')706 CALL writeg1d(1,nlay,wo,'w15uo ','wXo ')707 CALL writeg1d(1,nlay,zdhadj,'zdhadj ','zdhadj ')708 CALL writeg1d(1,nlay,pduadj,'pduadj ','pduadj ')709 CALL writeg1d(1,nlay,pdvadj,'pdvadj ','pdvadj ')710 CALL writeg1d(1,nlay,pdoadj,'pdoadj ','pdoadj ')711 CALL writeg1d(1,nlay,entr ,'entr ','entr ')712 CALL writeg1d(1,nlay,detr ,'detr ','detr ')713 CALL writeg1d(1,nlay,fm ,'fm ','fm ')714 715 CALL writeg1d(1,nlay,pdtadj,'pdtadj ','pdtadj ')716 CALL writeg1d(1,nlay,pplay,'pplay ','pplay ')717 CALL writeg1d(1,nlay,pplev,'pplev ','pplev ')718 c recalcul des flux en diagnostique...719 c print*,'PAS DE TEMPS ',ptimestep720 call dt2F(pplev,pplay,pt,pdtadj,wh)721 CALL writeg1d(1,nlay,wh,'wh2 ','wh2 ')722 #endif723 123 continue724 ! #define troisD725 #ifdef troisD726 c if (sorties) then727 print*,'Debut des wrgradsfi'728 729 c print*,'16 OK convect8'730 call wrgradsfi(1,nlay,wd,'wd ','wd ')731 call wrgradsfi(1,nlay,zwa,'wa ','wa ')732 call wrgradsfi(1,nlay,fracd,'fracd ','fracd ')733 call wrgradsfi(1,nlay,fraca,'fraca ','fraca ')734 call wrgradsfi(1,nlay,xxx,'xxx ','xxx ')735 call wrgradsfi(1,nlay,wa_moy,'wam ','wam ')736 c print*,'WA6 ',wa_moy737 call wrgradsfi(1,nlay,zla,'la ','la ')738 call wrgradsfi(1,nlay,zld,'ld ','ld ')739 call wrgradsfi(1,nlay,pt,'pt ','pt ')740 call wrgradsfi(1,nlay,zh,'zh ','zh ')741 call wrgradsfi(1,nlay,zha,'zha ','zha ')742 call wrgradsfi(1,nlay,zua,'zua ','zua ')743 call wrgradsfi(1,nlay,zva,'zva ','zva ')744 call wrgradsfi(1,nlay,zu,'zu ','zu ')745 call wrgradsfi(1,nlay,zv,'zv ','zv ')746 call wrgradsfi(1,nlay,zo,'zo ','zo ')747 call wrgradsfi(1,nlay,wh,'wh ','wh ')748 call wrgradsfi(1,nlay,wu,'wu ','wu ')749 call wrgradsfi(1,nlay,wv,'wv ','wv ')750 call wrgradsfi(1,nlay,wo,'wo ','wo ')751 call wrgradsfi(1,1,zmax,'zmax ','zmax ')752 call wrgradsfi(1,nlay,zdhadj,'zdhadj ','zdhadj ')753 call wrgradsfi(1,nlay,pduadj,'pduadj ','pduadj ')754 call wrgradsfi(1,nlay,pdvadj,'pdvadj ','pdvadj ')755 call wrgradsfi(1,nlay,pdoadj,'pdoadj ','pdoadj ')756 call wrgradsfi(1,nlay,entr,'entr ','entr ')757 call wrgradsfi(1,nlay,detr,'detr ','detr ')758 call wrgradsfi(1,nlay,fm,'fm ','fm ')759 call wrgradsfi(1,nlay,fmc,'fmc ','fmc ')760 call wrgradsfi(1,nlay,zw2,'zw2 ','zw2 ')761 call wrgradsfi(1,nlay,ztva,'ztva ','ztva ')762 call wrgradsfi(1,nlay,ztv,'ztv ','ztv ')763 764 call wrgradsfi(1,nlay,zo,'zo ','zo ')765 call wrgradsfi(1,nlay,larg_cons,'Lc ','Lc ')766 call wrgradsfi(1,nlay,larg_detr,'Ldetr ','Ldetr ')767 768 769 c print*,'17 OK convect8'770 771 do k=1,klev/10772 write(str2,'(i2.2)') k773 str10='wa'//str2774 do l=1,nlay775 do ig=1,ngrid776 zsortie(ig,l)=wa(ig,k,l)777 enddo778 enddo779 CALL wrgradsfi(1,nlay,zsortie,str10,str10)780 do l=1,nlay781 do ig=1,ngrid782 zsortie(ig,l)=larg_part(ig,k,l)783 enddo784 enddo785 str10='la'//str2786 CALL wrgradsfi(1,nlay,zsortie,str10,str10)787 enddo788 789 790 c print*,'18 OK convect8'791 c endif792 print*,'Fin des wrgradsfi'793 #endif794 795 712 endif 796 713 714 ! print*,'15 OK convect8' 715 716 797 717 c if(wa_moy(1,4).gt.1.e-10) stop 798 718 799 cprint*,'19 OK convect8'719 ! print*,'19 OK convect8' 800 720 return 801 721 end … … 925 845 real ratqsdiff(klon,klev) 926 846 real count_time 927 integer isplit,nsplit,ialt 928 parameter (nsplit=10) 929 data isplit/0/ 930 save isplit 931 c$OMP THREADPRIVATE(isplit) 847 integer ialt 932 848 933 849 logical sorties … … 2197 2113 enddo 2198 2114 2199 2115 ! print*,'12 OK convect8' 2200 2116 c print*,'WA4 ',wa_moy 2201 2117 cc------------------------------------------------------------------ … … 2347 2263 c enddo 2348 2264 2349 2265 ! print*,'14 OK convect8' 2350 2266 c------------------------------------------------------------------ 2351 2267 c Calculs pour les sorties … … 2433 2349 enddo 2434 2350 enddo 2435 cdeja fait2436 c do l=1,nlay2437 c do ig=1,ngrid2438 c detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)2439 c if (detr(ig,l).lt.0.) then2440 c entr(ig,l)=entr(ig,l)-detr(ig,l)2441 c detr(ig,l)=0.2442 c print*,'WARNING !!! detrainement negatif ',ig,l2443 c endif2444 c enddo2445 c enddo2446 2447 c print*,'15 OK convect8'2448 2449 isplit=isplit+12450 2451 2452 c #define und2453 goto 1232454 #ifdef und2455 CALL writeg1d(1,nlay,wd,'wd ','wd ')2456 CALL writeg1d(1,nlay,zwa,'wa ','wa ')2457 CALL writeg1d(1,nlay,fracd,'fracd ','fracd ')2458 CALL writeg1d(1,nlay,fraca,'fraca ','fraca ')2459 CALL writeg1d(1,nlay,wa_moy,'wam ','wam ')2460 CALL writeg1d(1,nlay,zla,'la ','la ')2461 CALL writeg1d(1,nlay,zld,'ld ','ld ')2462 CALL writeg1d(1,nlay,pt,'pt ','pt ')2463 CALL writeg1d(1,nlay,zh,'zh ','zh ')2464 CALL writeg1d(1,nlay,zha,'zha ','zha ')2465 CALL writeg1d(1,nlay,zu,'zu ','zu ')2466 CALL writeg1d(1,nlay,zv,'zv ','zv ')2467 CALL writeg1d(1,nlay,zo,'zo ','zo ')2468 CALL writeg1d(1,nlay,wh,'wh ','wh ')2469 CALL writeg1d(1,nlay,wu,'wu ','wu ')2470 CALL writeg1d(1,nlay,wv,'wv ','wv ')2471 CALL writeg1d(1,nlay,wo,'w15uo ','wXo ')2472 CALL writeg1d(1,nlay,zdhadj,'zdhadj ','zdhadj ')2473 CALL writeg1d(1,nlay,pduadj,'pduadj ','pduadj ')2474 CALL writeg1d(1,nlay,pdvadj,'pdvadj ','pdvadj ')2475 CALL writeg1d(1,nlay,pdoadj,'pdoadj ','pdoadj ')2476 CALL writeg1d(1,nlay,entr ,'entr ','entr ')2477 CALL writeg1d(1,nlay,detr ,'detr ','detr ')2478 CALL writeg1d(1,nlay,fm ,'fm ','fm ')2479 2480 CALL writeg1d(1,nlay,pdtadj,'pdtadj ','pdtadj ')2481 CALL writeg1d(1,nlay,pplay,'pplay ','pplay ')2482 CALL writeg1d(1,nlay,pplev,'pplev ','pplev ')2483 2484 c recalcul des flux en diagnostique...2485 c print*,'PAS DE TEMPS ',ptimestep2486 call dt2F(pplev,pplay,pt,pdtadj,wh)2487 CALL writeg1d(1,nlay,wh,'wh2 ','wh2 ')2488 #endif2489 123 continue2490 ! #define troisD2491 #ifdef troisD2492 c if (sorties) then2493 print*,'Debut des wrgradsfi'2494 2495 c print*,'16 OK convect8'2496 c call wrgradsfi(1,nlay,wd,'wd ','wd ')2497 c call wrgradsfi(1,nlay,zwa,'wa ','wa ')2498 call wrgradsfi(1,nlay,fracd,'fracd ','fracd ')2499 call wrgradsfi(1,nlay,fraca,'fraca ','fraca ')2500 c call wrgradsfi(1,nlay,xxx,'xxx ','xxx ')2501 c call wrgradsfi(1,nlay,wa_moy,'wam ','wam ')2502 c print*,'WA6 ',wa_moy2503 c call wrgradsfi(1,nlay,zla,'la ','la ')2504 c call wrgradsfi(1,nlay,zld,'ld ','ld ')2505 call wrgradsfi(1,nlay,pt,'pt ','pt ')2506 call wrgradsfi(1,nlay,zh,'zh ','zh ')2507 call wrgradsfi(1,nlay,zha,'zha ','zha ')2508 call wrgradsfi(1,nlay,zua,'zua ','zua ')2509 call wrgradsfi(1,nlay,zva,'zva ','zva ')2510 call wrgradsfi(1,nlay,zu,'zu ','zu ')2511 call wrgradsfi(1,nlay,zv,'zv ','zv ')2512 call wrgradsfi(1,nlay,zo,'zo ','zo ')2513 call wrgradsfi(1,nlay,wh,'wh ','wh ')2514 call wrgradsfi(1,nlay,wu,'wu ','wu ')2515 call wrgradsfi(1,nlay,wv,'wv ','wv ')2516 call wrgradsfi(1,nlay,wo,'wo ','wo ')2517 call wrgradsfi(1,1,zmax,'zmax ','zmax ')2518 call wrgradsfi(1,nlay,zdhadj,'zdhadj ','zdhadj ')2519 call wrgradsfi(1,nlay,pduadj,'pduadj ','pduadj ')2520 call wrgradsfi(1,nlay,pdvadj,'pdvadj ','pdvadj ')2521 call wrgradsfi(1,nlay,pdoadj,'pdoadj ','pdoadj ')2522 call wrgradsfi(1,nlay,entr,'entr ','entr ')2523 call wrgradsfi(1,nlay,detr,'detr ','detr ')2524 call wrgradsfi(1,nlay,fm,'fm ','fm ')2525 call wrgradsfi(1,nlay,fmc,'fmc ','fmc ')2526 call wrgradsfi(1,nlay,zw2,'zw2 ','zw2 ')2527 call wrgradsfi(1,nlay,w_est,'w_est ','w_est ')2528 con sort les moments2529 call wrgradsfi(1,nlay,thetath2,'zh2 ','zh2 ')2530 call wrgradsfi(1,nlay,wth2,'w2 ','w2 ')2531 call wrgradsfi(1,nlay,wth3,'w3 ','w3 ')2532 call wrgradsfi(1,nlay,q2,'q2 ','q2 ')2533 call wrgradsfi(1,nlay,dtheta,'dT ','dT ')2534 c2535 call wrgradsfi(1,nlay,zw_sec,'zw_sec ','zw_sec ')2536 call wrgradsfi(1,nlay,ztva,'ztva ','ztva ')2537 call wrgradsfi(1,nlay,ztv,'ztv ','ztv ')2538 call wrgradsfi(1,nlay,nu,'nu ','nu ')2539 2540 call wrgradsfi(1,nlay,zo,'zo ','zo ')2541 call wrgradsfi(1,nlay,zoa,'zoa ','zoa ')2542 c call wrgradsfi(1,nlay,larg_cons,'Lc ','Lc ')2543 c call wrgradsfi(1,nlay,larg_detr,'Ldetr ','Ldetr ')2544 2545 cAM:nouveaux diagnostiques2546 call wrgradsfi(1,nlay,zthl,'zthl ','zthl ')2547 call wrgradsfi(1,nlay,zta,'zta ','zta ')2548 call wrgradsfi(1,nlay,zl,'zl ','zl ')2549 call wrgradsfi(1,nlay,zdthladj,'zdthladj ',2550 s 'zdthladj ')2551 call wrgradsfi(1,nlay,ztla,'ztla ','ztla ')2552 call wrgradsfi(1,nlay,zqta,'zqta ','zqta ')2553 call wrgradsfi(1,nlay,zqla,'zqla ','zqla ')2554 call wrgradsfi(1,nlay,deltaz,'deltaz ','deltaz ')2555 cCR:nouveaux diagnostiques2556 call wrgradsfi(1,nlay,entr_star ,'entr_star ','entr_star ')2557 call wrgradsfi(1,nlay,detr_star ,'detr_star ','detr_star ')2558 call wrgradsfi(1,nlay,f_star ,'f_star ','f_star ')2559 call wrgradsfi(1,nlay,zqsat ,'zqsat ','zqsat ')2560 call wrgradsfi(1,nlay,zqsatth ,'qsath ','qsath ')2561 call wrgradsfi(1,nlay,alim_star ,'alim_star ','alim_star ')2562 call wrgradsfi(1,nlay,alim ,'alim ','alim ')2563 call wrgradsfi(1,1,f,'f ','f ')2564 call wrgradsfi(1,1,alim_star_tot,'a_s_t ','a_s_t ')2565 call wrgradsfi(1,1,alim_star2,'a_2 ','a_2 ')2566 call wrgradsfi(1,1,zmax,'zmax ','zmax ')2567 call wrgradsfi(1,1,zmax_sec,'z_sec ','z_sec ')2568 c call wrgradsfi(1,1,zmax_sec2,'zz_se ','zz_se ')2569 call wrgradsfi(1,1,zmix,'zmix ','zmix ')2570 call wrgradsfi(1,1,nivcon,'nivcon ','nivcon ')2571 call wrgradsfi(1,1,zcon,'zcon ','zcon ')2572 zsortie1d(:)=lmax(:)2573 call wrgradsfi(1,1,zsortie1d,'lmax ','lmax ')2574 call wrgradsfi(1,1,wmax,'wmax ','wmax ')2575 call wrgradsfi(1,1,wmax_sec,'w_sec ','w_sec ')2576 zsortie1d(:)=lmix(:)2577 call wrgradsfi(1,1,zsortie1d,'lmix ','lmix ')2578 zsortie1d(:)=lentr(:)2579 call wrgradsfi(1,1,zsortie1d,'lentr ','lentr ')2580 2581 c print*,'17 OK convect8'2582 2583 do k=1,klev/102584 write(str2,'(i2.2)') k2585 str10='wa'//str22586 do l=1,nlay2587 do ig=1,ngrid2588 zsortie(ig,l)=wa(ig,k,l)2589 enddo2590 enddo2591 CALL wrgradsfi(1,nlay,zsortie,str10,str10)2592 do l=1,nlay2593 do ig=1,ngrid2594 zsortie(ig,l)=larg_part(ig,k,l)2595 enddo2596 enddo2597 str10='la'//str22598 CALL wrgradsfi(1,nlay,zsortie,str10,str10)2599 enddo2600 2601 2602 c print*,'18 OK convect8'2603 c endif2604 print*,'Fin des wrgradsfi'2605 #endif2606 2351 2607 2352 endif 2608 2353 2609 c if(wa_moy(1,4).gt.1.e-10) stop 2610 2611 print*,'19 OK convect8' 2354 ! print*,'19 OK convect8' 2612 2355 return 2613 2356 end 2357 2614 2358 SUBROUTINE thermcell_eau(ngrid,nlay,ptimestep 2615 2359 s ,pplay,pplev,pphi … … 2710 2454 2711 2455 real count_time 2712 integer isplit,nsplit,ialt 2713 parameter (nsplit=10) 2714 data isplit/0/ 2715 save isplit 2716 c$OMP THREADPRIVATE(isplit) 2456 integer ialt 2717 2457 2718 2458 logical sorties … … 2861 2601 c --------------------------------------------------- 2862 2602 2863 2603 ! print*,'0 OK convect8' 2864 2604 2865 2605 DO 1010 l=1,nlay … … 3568 3308 c------------------------------------------------------------------ 3569 3309 3570 if(sorties) then 3310 return 3311 end 3312 3313 SUBROUTINE thermcell(ngrid,nlay,ptimestep 3314 s ,pplay,pplev,pphi 3315 s ,pu,pv,pt,po 3316 s ,pduadj,pdvadj,pdtadj,pdoadj 3317 s ,fm0,entr0 3318 c s ,pu_therm,pv_therm 3319 s ,r_aspect,l_mix,w2di,tho) 3320 3321 USE dimphy 3322 IMPLICIT NONE 3323 3324 c======================================================================= 3325 c 3326 c Calcul du transport verticale dans la couche limite en presence 3327 c de "thermiques" explicitement representes 3328 c 3329 c Réécriture à partir d'un listing papier à Habas, le 14/02/00 3330 c 3331 c le thermique est supposé homogène et dissipé par mélange avec 3332 c son environnement. la longueur l_mix contrôle l'efficacité du 3333 c mélange 3334 c 3335 c Le calcul du transport des différentes espèces se fait en prenant 3336 c en compte: 3337 c 1. un flux de masse montant 3338 c 2. un flux de masse descendant 3339 c 3. un entrainement 3340 c 4. un detrainement 3341 c 3342 c======================================================================= 3343 3344 c----------------------------------------------------------------------- 3345 c declarations: 3346 c ------------- 3347 3348 #include "dimensions.h" 3349 cccc#include "dimphy.h" 3350 #include "YOMCST.h" 3351 3352 c arguments: 3353 c ---------- 3354 3355 INTEGER ngrid,nlay,w2di 3356 REAL tho 3357 real ptimestep,l_mix,r_aspect 3358 REAL pt(ngrid,nlay),pdtadj(ngrid,nlay) 3359 REAL pu(ngrid,nlay),pduadj(ngrid,nlay) 3360 REAL pv(ngrid,nlay),pdvadj(ngrid,nlay) 3361 REAL po(ngrid,nlay),pdoadj(ngrid,nlay) 3362 REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1) 3363 real pphi(ngrid,nlay) 3364 3365 integer idetr 3366 save idetr 3367 data idetr/3/ 3368 c$OMP THREADPRIVATE(idetr) 3369 3370 c local: 3371 c ------ 3372 3373 INTEGER ig,k,l,lmaxa(klon),lmix(klon) 3374 real zsortie1d(klon) 3375 c CR: on remplace lmax(klon,klev+1) 3376 INTEGER lmax(klon),lmin(klon),lentr(klon) 3377 real linter(klon) 3378 real zmix(klon), fracazmix(klon) 3379 c RC 3380 real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz 3381 3382 real zlev(klon,klev+1),zlay(klon,klev) 3383 REAL zh(klon,klev),zdhadj(klon,klev) 3384 REAL ztv(klon,klev) 3385 real zu(klon,klev),zv(klon,klev),zo(klon,klev) 3386 REAL wh(klon,klev+1) 3387 real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1) 3388 real zla(klon,klev+1) 3389 real zwa(klon,klev+1) 3390 real zld(klon,klev+1) 3391 real zwd(klon,klev+1) 3392 real zsortie(klon,klev) 3393 real zva(klon,klev) 3394 real zua(klon,klev) 3395 real zoa(klon,klev) 3396 3397 real zha(klon,klev) 3398 real wa_moy(klon,klev+1) 3399 real fraca(klon,klev+1) 3400 real fracc(klon,klev+1) 3401 real zf,zf2 3402 real thetath2(klon,klev),wth2(klon,klev) 3403 ! common/comtherm/thetath2,wth2 3404 3405 real count_time 3406 integer ialt 3407 3408 logical sorties 3409 real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev) 3410 real zpspsk(klon,klev) 3411 3412 c real wmax(klon,klev),wmaxa(klon) 3413 real wmax(klon),wmaxa(klon) 3414 real wa(klon,klev,klev+1) 3415 real wd(klon,klev+1) 3416 real larg_part(klon,klev,klev+1) 3417 real fracd(klon,klev+1) 3418 real xxx(klon,klev+1) 3419 real larg_cons(klon,klev+1) 3420 real larg_detr(klon,klev+1) 3421 real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev) 3422 real pu_therm(klon,klev),pv_therm(klon,klev) 3423 real fm(klon,klev+1),entr(klon,klev) 3424 real fmc(klon,klev+1) 3425 3426 cCR:nouvelles variables 3427 real f_star(klon,klev+1),entr_star(klon,klev) 3428 real entr_star_tot(klon),entr_star2(klon) 3429 real f(klon), f0(klon) 3430 real zlevinter(klon) 3431 logical first 3432 data first /.false./ 3433 save first 3434 c$OMP THREADPRIVATE(first) 3435 cRC 3436 3437 character*2 str2 3438 character*10 str10 3439 3440 character (len=20) :: modname='thermcell' 3441 character (len=80) :: abort_message 3442 3443 LOGICAL vtest(klon),down 3444 3445 EXTERNAL SCOPY 3446 3447 integer ncorrec,ll 3448 save ncorrec 3449 data ncorrec/0/ 3450 c$OMP THREADPRIVATE(ncorrec) 3451 3452 c 3453 c----------------------------------------------------------------------- 3454 c initialisation: 3455 c --------------- 3456 c 3457 sorties=.true. 3458 IF(ngrid.NE.klon) THEN 3459 PRINT* 3460 PRINT*,'STOP dans convadj' 3461 PRINT*,'ngrid =',ngrid 3462 PRINT*,'klon =',klon 3463 ENDIF 3464 c 3465 c----------------------------------------------------------------------- 3466 c incrementation eventuelle de tendances precedentes: 3467 c --------------------------------------------------- 3468 3469 ! print*,'0 OK convect8' 3470 3471 DO 1010 l=1,nlay 3472 DO 1015 ig=1,ngrid 3473 zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA 3474 zh(ig,l)=pt(ig,l)/zpspsk(ig,l) 3475 zu(ig,l)=pu(ig,l) 3476 zv(ig,l)=pv(ig,l) 3477 zo(ig,l)=po(ig,l) 3478 ztv(ig,l)=zh(ig,l)*(1.+0.61*zo(ig,l)) 3479 1015 CONTINUE 3480 1010 CONTINUE 3481 3482 ! print*,'1 OK convect8' 3483 c -------------------- 3484 c 3485 c 3486 c + + + + + + + + + + + 3487 c 3488 c 3489 c wa, fraca, wd, fracd -------------------- zlev(2), rhobarz 3490 c wh,wt,wo ... 3491 c 3492 c + + + + + + + + + + + zh,zu,zv,zo,rho 3493 c 3494 c 3495 c -------------------- zlev(1) 3496 c \\\\\\\\\\\\\\\\\\\\ 3497 c 3498 c 3499 3500 c----------------------------------------------------------------------- 3501 c Calcul des altitudes des couches 3502 c----------------------------------------------------------------------- 3503 3504 do l=2,nlay 3505 do ig=1,ngrid 3506 zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG 3507 enddo 3508 enddo 3509 do ig=1,ngrid 3510 zlev(ig,1)=0. 3511 zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG 3512 enddo 3571 3513 do l=1,nlay 3572 3514 do ig=1,ngrid 3573 zla(ig,l)=(1.-fracd(ig,l))*zmax(ig) 3574 zld(ig,l)=fracd(ig,l)*zmax(ig) 3575 if(1.-fracd(ig,l).gt.1.e-10) 3576 s zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l)) 3577 enddo 3578 enddo 3515 zlay(ig,l)=pphi(ig,l)/RG 3516 enddo 3517 enddo 3518 3519 c print*,'2 OK convect8' 3520 c----------------------------------------------------------------------- 3521 c Calcul des densites 3522 c----------------------------------------------------------------------- 3579 3523 3580 3524 do l=1,nlay 3525 do ig=1,ngrid 3526 rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)*RD*zh(ig,l)) 3527 enddo 3528 enddo 3529 3530 do l=2,nlay 3531 do ig=1,ngrid 3532 rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1)) 3533 enddo 3534 enddo 3535 3536 do k=1,nlay 3537 do l=1,nlay+1 3538 do ig=1,ngrid 3539 wa(ig,k,l)=0. 3540 enddo 3541 enddo 3542 enddo 3543 3544 c print*,'3 OK convect8' 3545 c------------------------------------------------------------------ 3546 c Calcul de w2, quarre de w a partir de la cape 3547 c a partir de w2, on calcule wa, vitesse de l'ascendance 3548 c 3549 c ATTENTION: Dans cette version, pour cause d'economie de memoire, 3550 c w2 est stoke dans wa 3551 c 3552 c ATTENTION: dans convect8, on n'utilise le calcule des wa 3553 c independants par couches que pour calculer l'entrainement 3554 c a la base et la hauteur max de l'ascendance. 3555 c 3556 c Indicages: 3557 c l'ascendance provenant du niveau k traverse l'interface l avec 3558 c une vitesse wa(k,l). 3559 c 3560 c -------------------- 3561 c 3562 c + + + + + + + + + + 3563 c 3564 c wa(k,l) ---- -------------------- l 3565 c /\ 3566 c /||\ + + + + + + + + + + 3567 c || 3568 c || -------------------- 3569 c || 3570 c || + + + + + + + + + + 3571 c || 3572 c || -------------------- 3573 c ||__ 3574 c |___ + + + + + + + + + + k 3575 c 3576 c -------------------- 3577 c 3578 c 3579 c 3580 c------------------------------------------------------------------ 3581 3582 cCR: ponderation entrainement des couches instables 3583 cdef des entr_star tels que entr=f*entr_star 3584 do l=1,klev 3585 do ig=1,ngrid 3586 entr_star(ig,l)=0. 3587 enddo 3588 enddo 3589 c determination de la longueur de la couche d entrainement 3590 do ig=1,ngrid 3591 lentr(ig)=1 3592 enddo 3593 3594 con ne considere que les premieres couches instables 3595 do k=nlay-2,1,-1 3596 do ig=1,ngrid 3597 if (ztv(ig,k).gt.ztv(ig,k+1).and. 3598 s ztv(ig,k+1).le.ztv(ig,k+2)) then 3599 lentr(ig)=k 3600 endif 3601 enddo 3602 enddo 3603 3604 c determination du lmin: couche d ou provient le thermique 3605 do ig=1,ngrid 3606 lmin(ig)=1 3607 enddo 3608 do ig=1,ngrid 3609 do l=nlay,2,-1 3610 if (ztv(ig,l-1).gt.ztv(ig,l)) then 3611 lmin(ig)=l-1 3612 endif 3613 enddo 3614 enddo 3615 c 3616 c definition de l'entrainement des couches 3617 do l=1,klev-1 3618 do ig=1,ngrid 3619 if (ztv(ig,l).gt.ztv(ig,l+1).and. 3620 s l.ge.lmin(ig).and.l.le.lentr(ig)) then 3621 entr_star(ig,l)=(ztv(ig,l)-ztv(ig,l+1))* 3622 s (zlev(ig,l+1)-zlev(ig,l)) 3623 endif 3624 enddo 3625 enddo 3626 c pas de thermique si couches 1->5 stables 3627 do ig=1,ngrid 3628 if (lmin(ig).gt.5) then 3629 do l=1,klev 3630 entr_star(ig,l)=0. 3631 enddo 3632 endif 3633 enddo 3634 c calcul de l entrainement total 3635 do ig=1,ngrid 3636 entr_star_tot(ig)=0. 3637 enddo 3638 do ig=1,ngrid 3639 do k=1,klev 3640 entr_star_tot(ig)=entr_star_tot(ig)+entr_star(ig,k) 3641 enddo 3642 enddo 3643 c 3644 print*,'fin calcul entr_star' 3645 do k=1,klev 3646 do ig=1,ngrid 3647 ztva(ig,k)=ztv(ig,k) 3648 enddo 3649 enddo 3650 cRC 3651 c print*,'7 OK convect8' 3652 do k=1,klev+1 3653 do ig=1,ngrid 3654 zw2(ig,k)=0. 3655 fmc(ig,k)=0. 3656 cCR 3657 f_star(ig,k)=0. 3658 cRC 3659 larg_cons(ig,k)=0. 3660 larg_detr(ig,k)=0. 3661 wa_moy(ig,k)=0. 3662 enddo 3663 enddo 3664 3665 c print*,'8 OK convect8' 3666 do ig=1,ngrid 3667 linter(ig)=1. 3668 lmaxa(ig)=1 3669 lmix(ig)=1 3670 wmaxa(ig)=0. 3671 enddo 3672 3673 cCR: 3674 do l=1,nlay-2 3675 do ig=1,ngrid 3676 if (ztv(ig,l).gt.ztv(ig,l+1) 3677 s .and.entr_star(ig,l).gt.1.e-10 3678 s .and.zw2(ig,l).lt.1e-10) then 3679 f_star(ig,l+1)=entr_star(ig,l) 3680 ctest:calcul de dteta 3681 zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1) 3682 s *(zlev(ig,l+1)-zlev(ig,l)) 3683 s *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l)) 3684 larg_detr(ig,l)=0. 3685 else if ((zw2(ig,l).ge.1e-10).and. 3686 s (f_star(ig,l)+entr_star(ig,l).gt.1.e-10)) then 3687 f_star(ig,l+1)=f_star(ig,l)+entr_star(ig,l) 3688 ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l) 3689 s *ztv(ig,l))/f_star(ig,l+1) 3690 zw2(ig,l+1)=zw2(ig,l)*(f_star(ig,l)/f_star(ig,l+1))**2+ 3691 s 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) 3692 s *(zlev(ig,l+1)-zlev(ig,l)) 3693 endif 3694 c determination de zmax continu par interpolation lineaire 3695 if (zw2(ig,l+1).lt.0.) then 3696 ctest 3697 if (abs(zw2(ig,l+1)-zw2(ig,l)).lt.1e-10) then 3698 print*,'pb linter' 3699 endif 3700 linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) 3701 s -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l)) 3702 zw2(ig,l+1)=0. 3703 lmaxa(ig)=l 3704 else 3705 if (zw2(ig,l+1).lt.0.) then 3706 print*,'pb1 zw2<0' 3707 endif 3708 wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) 3709 endif 3710 if (wa_moy(ig,l+1).gt.wmaxa(ig)) then 3711 c lmix est le niveau de la couche ou w (wa_moy) est maximum 3712 lmix(ig)=l+1 3713 wmaxa(ig)=wa_moy(ig,l+1) 3714 endif 3715 enddo 3716 enddo 3717 print*,'fin calcul zw2' 3718 c 3719 c Calcul de la couche correspondant a la hauteur du thermique 3720 do ig=1,ngrid 3721 lmax(ig)=lentr(ig) 3722 enddo 3723 do ig=1,ngrid 3724 do l=nlay,lentr(ig)+1,-1 3725 if (zw2(ig,l).le.1.e-10) then 3726 lmax(ig)=l-1 3727 endif 3728 enddo 3729 enddo 3730 c pas de thermique si couches 1->5 stables 3731 do ig=1,ngrid 3732 if (lmin(ig).gt.5) then 3733 lmax(ig)=1 3734 lmin(ig)=1 3735 endif 3736 enddo 3737 c 3738 c Determination de zw2 max 3739 do ig=1,ngrid 3740 wmax(ig)=0. 3741 enddo 3742 3743 do l=1,nlay 3744 do ig=1,ngrid 3745 if (l.le.lmax(ig)) then 3746 if (zw2(ig,l).lt.0.)then 3747 print*,'pb2 zw2<0' 3748 endif 3749 zw2(ig,l)=sqrt(zw2(ig,l)) 3750 wmax(ig)=max(wmax(ig),zw2(ig,l)) 3751 else 3752 zw2(ig,l)=0. 3753 endif 3754 enddo 3755 enddo 3756 3757 c Longueur caracteristique correspondant a la hauteur des thermiques. 3758 do ig=1,ngrid 3759 zmax(ig)=0. 3760 zlevinter(ig)=zlev(ig,1) 3761 enddo 3762 do ig=1,ngrid 3763 c calcul de zlevinter 3764 zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))* 3765 s linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1) 3766 s -zlev(ig,lmax(ig))) 3767 zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig))) 3768 enddo 3769 3770 print*,'avant fermeture' 3771 c Fermeture,determination de f 3772 do ig=1,ngrid 3773 entr_star2(ig)=0. 3774 enddo 3775 do ig=1,ngrid 3776 if (entr_star_tot(ig).LT.1.e-10) then 3777 f(ig)=0. 3778 else 3779 do k=lmin(ig),lentr(ig) 3780 entr_star2(ig)=entr_star2(ig)+entr_star(ig,k)**2 3781 s /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))) 3782 enddo 3783 c Nouvelle fermeture 3784 f(ig)=wmax(ig)/(max(500.,zmax(ig))*r_aspect 3785 s *entr_star2(ig))*entr_star_tot(ig) 3786 ctest 3787 c if (first) then 3788 c f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig) 3789 c s *wmax(ig)) 3790 c endif 3791 endif 3792 c f0(ig)=f(ig) 3793 c first=.true. 3794 enddo 3795 print*,'apres fermeture' 3796 3797 c Calcul de l'entrainement 3798 do k=1,klev 3799 do ig=1,ngrid 3800 entr(ig,k)=f(ig)*entr_star(ig,k) 3801 enddo 3802 enddo 3803 c Calcul des flux 3804 do ig=1,ngrid 3805 do l=1,lmax(ig)-1 3806 fmc(ig,l+1)=fmc(ig,l)+entr(ig,l) 3807 enddo 3808 enddo 3809 3810 cRC 3811 3812 3813 c print*,'9 OK convect8' 3814 c print*,'WA1 ',wa_moy 3815 3816 c determination de l'indice du debut de la mixed layer ou w decroit 3817 3818 c calcul de la largeur de chaque ascendance dans le cas conservatif. 3819 c dans ce cas simple, on suppose que la largeur de l'ascendance provenant 3820 c d'une couche est égale à la hauteur de la couche alimentante. 3821 c La vitesse maximale dans l'ascendance est aussi prise comme estimation 3822 c de la vitesse d'entrainement horizontal dans la couche alimentante. 3823 3824 do l=2,nlay 3825 do ig=1,ngrid 3826 if (l.le.lmaxa(ig)) then 3827 zw=max(wa_moy(ig,l),1.e-10) 3828 larg_cons(ig,l)=zmax(ig)*r_aspect 3829 s *fmc(ig,l)/(rhobarz(ig,l)*zw) 3830 endif 3831 enddo 3832 enddo 3833 3834 do l=2,nlay 3835 do ig=1,ngrid 3836 if (l.le.lmaxa(ig)) then 3837 c if (idetr.eq.0) then 3838 c cette option est finalement en dur. 3839 if ((l_mix*zlev(ig,l)).lt.0.)then 3840 print*,'pb l_mix*zlev<0' 3841 endif 3842 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l)) 3843 c else if (idetr.eq.1) then 3844 c larg_detr(ig,l)=larg_cons(ig,l) 3845 c s *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig)) 3846 c else if (idetr.eq.2) then 3847 c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l)) 3848 c s *sqrt(wa_moy(ig,l)) 3849 c else if (idetr.eq.4) then 3850 c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l)) 3851 c s *wa_moy(ig,l) 3852 c endif 3853 endif 3854 enddo 3855 enddo 3856 3857 c print*,'10 OK convect8' 3858 c print*,'WA2 ',wa_moy 3859 c calcul de la fraction de la maille concernée par l'ascendance en tenant 3860 c compte de l'epluchage du thermique. 3861 c 3862 cCR def de zmix continu (profil parabolique des vitesses) 3863 do ig=1,ngrid 3864 if (lmix(ig).gt.1.) then 3865 c test 3866 if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig))) 3867 s *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1))) 3868 s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1)) 3869 s *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10) 3870 s then 3871 c 3872 zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig))) 3873 s *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2) 3874 s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1)) 3875 s *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2)) 3876 s /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig))) 3877 s *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1))) 3878 s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1)) 3879 s *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig)))))) 3880 else 3881 zmix(ig)=zlev(ig,lmix(ig)) 3882 print*,'pb zmix' 3883 endif 3884 else 3885 zmix(ig)=0. 3886 endif 3887 ctest 3888 if ((zmax(ig)-zmix(ig)).lt.0.) then 3889 zmix(ig)=0.99*zmax(ig) 3890 c print*,'pb zmix>zmax' 3891 endif 3892 enddo 3893 c 3894 c calcul du nouveau lmix correspondant 3895 do ig=1,ngrid 3896 do l=1,klev 3897 if (zmix(ig).ge.zlev(ig,l).and. 3898 s zmix(ig).lt.zlev(ig,l+1)) then 3899 lmix(ig)=l 3900 endif 3901 enddo 3902 enddo 3903 c 3904 do l=2,nlay 3905 do ig=1,ngrid 3906 if(larg_cons(ig,l).gt.1.) then 3907 c print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK' 3908 fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l)) 3909 s /(r_aspect*zmax(ig)) 3910 c test 3911 fraca(ig,l)=max(fraca(ig,l),0.) 3912 fraca(ig,l)=min(fraca(ig,l),0.5) 3913 fracd(ig,l)=1.-fraca(ig,l) 3914 fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig)) 3915 else 3916 c wa_moy(ig,l)=0. 3917 fraca(ig,l)=0. 3918 fracc(ig,l)=0. 3919 fracd(ig,l)=1. 3920 endif 3921 enddo 3922 enddo 3923 cCR: calcul de fracazmix 3924 do ig=1,ngrid 3925 fracazmix(ig)=(fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/ 3926 s (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig) 3927 s +fraca(ig,lmix(ig))-zlev(ig,lmix(ig))*(fraca(ig,lmix(ig)+1) 3928 s -fraca(ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig))) 3929 enddo 3930 c 3931 do l=2,nlay 3932 do ig=1,ngrid 3933 if(larg_cons(ig,l).gt.1.) then 3934 if (l.gt.lmix(ig)) then 3935 ctest 3936 if (zmax(ig)-zmix(ig).lt.1.e-10) then 3937 c print*,'pb xxx' 3938 xxx(ig,l)=(lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig)) 3939 else 3940 xxx(ig,l)=(zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig)) 3941 endif 3942 if (idetr.eq.0) then 3943 fraca(ig,l)=fracazmix(ig) 3944 else if (idetr.eq.1) then 3945 fraca(ig,l)=fracazmix(ig)*xxx(ig,l) 3946 else if (idetr.eq.2) then 3947 fraca(ig,l)=fracazmix(ig)*(1.-(1.-xxx(ig,l))**2) 3948 else 3949 fraca(ig,l)=fracazmix(ig)*xxx(ig,l)**2 3950 endif 3951 c print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL' 3952 fraca(ig,l)=max(fraca(ig,l),0.) 3953 fraca(ig,l)=min(fraca(ig,l),0.5) 3954 fracd(ig,l)=1.-fraca(ig,l) 3955 fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig)) 3956 endif 3957 endif 3958 enddo 3959 enddo 3960 3961 print*,'fin calcul fraca' 3962 c print*,'11 OK convect8' 3963 c print*,'Ea3 ',wa_moy 3964 c------------------------------------------------------------------ 3965 c Calcul de fracd, wd 3966 c somme wa - wd = 0 3967 c------------------------------------------------------------------ 3968 3969 3970 do ig=1,ngrid 3971 fm(ig,1)=0. 3972 fm(ig,nlay+1)=0. 3973 enddo 3974 3975 do l=2,nlay 3976 do ig=1,ngrid 3977 fm(ig,l)=fraca(ig,l)*wa_moy(ig,l)*rhobarz(ig,l) 3978 cCR:test 3979 if (entr(ig,l-1).lt.1e-10.and.fm(ig,l).gt.fm(ig,l-1) 3980 s .and.l.gt.lmix(ig)) then 3981 fm(ig,l)=fm(ig,l-1) 3982 c write(1,*)'ajustement fm, l',l 3983 endif 3984 c write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l) 3985 cRC 3986 enddo 3987 do ig=1,ngrid 3988 if(fracd(ig,l).lt.0.1) then 3989 abort_message = 'fracd trop petit' 3990 CALL abort_gcm (modname,abort_message,1) 3991 else 3992 c vitesse descendante "diagnostique" 3993 wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l)) 3994 endif 3995 enddo 3996 enddo 3997 3998 do l=1,nlay 3999 do ig=1,ngrid 4000 c masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) 4001 masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG 4002 enddo 4003 enddo 4004 4005 ! print*,'12 OK convect8' 4006 c print*,'WA4 ',wa_moy 4007 cc------------------------------------------------------------------ 4008 c calcul du transport vertical 4009 c------------------------------------------------------------------ 4010 4011 go to 4444 4012 c print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep 4013 do l=2,nlay-1 4014 do ig=1,ngrid 4015 if(fm(ig,l+1)*ptimestep.gt.masse(ig,l) 4016 s .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then 4017 c print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM=' 4018 c s ,fm(ig,l+1)*ptimestep 4019 c s ,' M=',masse(ig,l),masse(ig,l+1) 4020 endif 4021 enddo 4022 enddo 4023 4024 do l=1,nlay 4025 do ig=1,ngrid 4026 if(entr(ig,l)*ptimestep.gt.masse(ig,l)) then 4027 c print*,'WARN!!! E>M ig=',ig,' l=',l,' E==' 4028 c s ,entr(ig,l)*ptimestep 4029 c s ,' M=',masse(ig,l) 4030 endif 4031 enddo 4032 enddo 4033 4034 do l=1,nlay 4035 do ig=1,ngrid 4036 if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then 4037 c print*,'WARN!!! fm exagere ig=',ig,' l=',l 4038 c s ,' FM=',fm(ig,l) 4039 endif 4040 if(.not.masse(ig,l).ge.1.e-10 4041 s .or..not.masse(ig,l).le.1.e4) then 4042 c print*,'WARN!!! masse exagere ig=',ig,' l=',l 4043 c s ,' M=',masse(ig,l) 4044 c print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)', 4045 c s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l) 4046 c print*,'zlev(ig,l+1),zlev(ig,l)' 4047 c s ,zlev(ig,l+1),zlev(ig,l) 4048 c print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)' 4049 c s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1) 4050 endif 4051 if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then 4052 c print*,'WARN!!! entr exagere ig=',ig,' l=',l 4053 c s ,' E=',entr(ig,l) 4054 endif 4055 enddo 4056 enddo 4057 4058 4444 continue 4059 4060 cCR:redefinition du entr 4061 do l=1,nlay 3581 4062 do ig=1,ngrid 3582 4063 detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1) … … 3588 4069 enddo 3589 4070 enddo 4071 cRC 4072 if (w2di.eq.1) then 4073 fm0=fm0+ptimestep*(fm-fm0)/tho 4074 entr0=entr0+ptimestep*(entr-entr0)/tho 4075 else 4076 fm0=fm 4077 entr0=entr 4078 endif 4079 4080 if (1.eq.1) then 4081 call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse 4082 . ,zh,zdhadj,zha) 4083 call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse 4084 . ,zo,pdoadj,zoa) 4085 else 4086 call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca 4087 . ,zh,zdhadj,zha) 4088 call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca 4089 . ,zo,pdoadj,zoa) 4090 endif 4091 4092 if (1.eq.0) then 4093 call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse 4094 . ,fraca,zmax 4095 . ,zu,zv,pduadj,pdvadj,zua,zva) 4096 else 4097 call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse 4098 . ,zu,pduadj,zua) 4099 call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse 4100 . ,zv,pdvadj,zva) 4101 endif 4102 4103 do l=1,nlay 4104 do ig=1,ngrid 4105 zf=0.5*(fracc(ig,l)+fracc(ig,l+1)) 4106 zf2=zf/(1.-zf) 4107 thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l))**2 4108 wth2(ig,l)=zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2 4109 enddo 4110 enddo 4111 4112 4113 4114 c print*,'13 OK convect8' 4115 c print*,'WA5 ',wa_moy 4116 do l=1,nlay 4117 do ig=1,ngrid 4118 pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l) 4119 enddo 4120 enddo 4121 4122 4123 c do l=1,nlay 4124 c do ig=1,ngrid 4125 c if(abs(pdtadj(ig,l))*86400..gt.500.) then 4126 c print*,'WARN!!! ig=',ig,' l=',l 4127 c s ,' pdtadj=',pdtadj(ig,l) 4128 c endif 4129 c if(abs(pdoadj(ig,l))*86400..gt.1.) then 4130 c print*,'WARN!!! ig=',ig,' l=',l 4131 c s ,' pdoadj=',pdoadj(ig,l) 4132 c endif 4133 c enddo 4134 c enddo 4135 4136 ! print*,'14 OK convect8' 4137 c------------------------------------------------------------------ 4138 c Calculs pour les sorties 4139 c------------------------------------------------------------------ 4140 4141 if(sorties) then 4142 do l=1,nlay 4143 do ig=1,ngrid 4144 zla(ig,l)=(1.-fracd(ig,l))*zmax(ig) 4145 zld(ig,l)=fracd(ig,l)*zmax(ig) 4146 if(1.-fracd(ig,l).gt.1.e-10) 4147 s zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l)) 4148 enddo 4149 enddo 4150 4151 cdeja fait 4152 c do l=1,nlay 4153 c do ig=1,ngrid 4154 c detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1) 4155 c if (detr(ig,l).lt.0.) then 4156 c entr(ig,l)=entr(ig,l)-detr(ig,l) 4157 c detr(ig,l)=0. 4158 c print*,'WARNING !!! detrainement negatif ',ig,l 4159 c endif 4160 c enddo 4161 c enddo 3590 4162 3591 4163 c print*,'15 OK convect8' 3592 4164 3593 isplit=isplit+13594 3595 4165 3596 4166 c #define und 3597 4167 goto 123 3598 4168 #ifdef und 3599 4169 CALL writeg1d(1,nlay,wd,'wd ','wd ') … … 3632 4202 #endif 3633 4203 123 continue 3634 !#define troisD4204 #define troisD 3635 4205 #ifdef troisD 3636 4206 c if (sorties) then … … 3676 4246 call wrgradsfi(1,nlay,larg_detr,'Ldetr ','Ldetr ') 3677 4247 3678 cAM:nouveaux diagnostiques3679 call wrgradsfi(1,nlay,zthl,'zthl ','zthl ')3680 call wrgradsfi(1,nlay,zta,'zta ','zta ')3681 call wrgradsfi(1,nlay,zl,'zl ','zl ')3682 call wrgradsfi(1,nlay,zdthladj,'zdthladj ',3683 s 'zdthladj ')3684 call wrgradsfi(1,nlay,ztla,'ztla ','ztla ')3685 call wrgradsfi(1,nlay,zqta,'zqta ','zqta ')3686 call wrgradsfi(1,nlay,zqla,'zqla ','zqla ')3687 4248 cCR:nouveaux diagnostiques 3688 4249 call wrgradsfi(1,nlay,entr_star ,'entr_star ','entr_star ') … … 3728 4289 c if(wa_moy(1,4).gt.1.e-10) stop 3729 4290 3730 c print*,'19 OK convect8' 3731 return 3732 end 3733 3734 SUBROUTINE thermcell(ngrid,nlay,ptimestep 3735 s ,pplay,pplev,pphi 3736 s ,pu,pv,pt,po 3737 s ,pduadj,pdvadj,pdtadj,pdoadj 3738 s ,fm0,entr0 3739 c s ,pu_therm,pv_therm 3740 s ,r_aspect,l_mix,w2di,tho) 3741 3742 USE dimphy 3743 IMPLICIT NONE 3744 3745 c======================================================================= 3746 c 3747 c Calcul du transport verticale dans la couche limite en presence 3748 c de "thermiques" explicitement representes 3749 c 3750 c Réécriture à partir d'un listing papier à Habas, le 14/02/00 3751 c 3752 c le thermique est supposé homogène et dissipé par mélange avec 3753 c son environnement. la longueur l_mix contrôle l'efficacité du 3754 c mélange 3755 c 3756 c Le calcul du transport des différentes espèces se fait en prenant 3757 c en compte: 3758 c 1. un flux de masse montant 3759 c 2. un flux de masse descendant 3760 c 3. un entrainement 3761 c 4. un detrainement 3762 c 3763 c======================================================================= 3764 3765 c----------------------------------------------------------------------- 3766 c declarations: 3767 c ------------- 3768 3769 #include "dimensions.h" 3770 cccc#include "dimphy.h" 3771 #include "YOMCST.h" 3772 3773 c arguments: 3774 c ---------- 3775 3776 INTEGER ngrid,nlay,w2di 3777 REAL tho 3778 real ptimestep,l_mix,r_aspect 3779 REAL pt(ngrid,nlay),pdtadj(ngrid,nlay) 3780 REAL pu(ngrid,nlay),pduadj(ngrid,nlay) 3781 REAL pv(ngrid,nlay),pdvadj(ngrid,nlay) 3782 REAL po(ngrid,nlay),pdoadj(ngrid,nlay) 3783 REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1) 3784 real pphi(ngrid,nlay) 3785 3786 integer idetr 3787 save idetr 3788 data idetr/3/ 3789 c$OMP THREADPRIVATE(idetr) 3790 3791 c local: 3792 c ------ 3793 3794 INTEGER ig,k,l,lmaxa(klon),lmix(klon) 3795 real zsortie1d(klon) 3796 c CR: on remplace lmax(klon,klev+1) 3797 INTEGER lmax(klon),lmin(klon),lentr(klon) 3798 real linter(klon) 3799 real zmix(klon), fracazmix(klon) 3800 c RC 3801 real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz 3802 3803 real zlev(klon,klev+1),zlay(klon,klev) 3804 REAL zh(klon,klev),zdhadj(klon,klev) 3805 REAL ztv(klon,klev) 3806 real zu(klon,klev),zv(klon,klev),zo(klon,klev) 3807 REAL wh(klon,klev+1) 3808 real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1) 3809 real zla(klon,klev+1) 3810 real zwa(klon,klev+1) 3811 real zld(klon,klev+1) 3812 real zwd(klon,klev+1) 3813 real zsortie(klon,klev) 3814 real zva(klon,klev) 3815 real zua(klon,klev) 3816 real zoa(klon,klev) 3817 3818 real zha(klon,klev) 3819 real wa_moy(klon,klev+1) 3820 real fraca(klon,klev+1) 3821 real fracc(klon,klev+1) 3822 real zf,zf2 3823 real thetath2(klon,klev),wth2(klon,klev) 3824 ! common/comtherm/thetath2,wth2 3825 3826 real count_time 3827 integer isplit,nsplit,ialt 3828 parameter (nsplit=10) 3829 data isplit/0/ 3830 save isplit 3831 c$OMP THREADPRIVATE(isplit) 3832 3833 logical sorties 3834 real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev) 3835 real zpspsk(klon,klev) 3836 3837 c real wmax(klon,klev),wmaxa(klon) 3838 real wmax(klon),wmaxa(klon) 3839 real wa(klon,klev,klev+1) 3840 real wd(klon,klev+1) 3841 real larg_part(klon,klev,klev+1) 3842 real fracd(klon,klev+1) 3843 real xxx(klon,klev+1) 3844 real larg_cons(klon,klev+1) 3845 real larg_detr(klon,klev+1) 3846 real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev) 3847 real pu_therm(klon,klev),pv_therm(klon,klev) 3848 real fm(klon,klev+1),entr(klon,klev) 3849 real fmc(klon,klev+1) 3850 3851 cCR:nouvelles variables 3852 real f_star(klon,klev+1),entr_star(klon,klev) 3853 real entr_star_tot(klon),entr_star2(klon) 3854 real f(klon), f0(klon) 3855 real zlevinter(klon) 3856 logical first 3857 data first /.false./ 3858 save first 3859 c$OMP THREADPRIVATE(first) 3860 cRC 3861 3862 character*2 str2 3863 character*10 str10 3864 3865 character (len=20) :: modname='thermcell' 3866 character (len=80) :: abort_message 3867 3868 LOGICAL vtest(klon),down 3869 3870 EXTERNAL SCOPY 3871 3872 integer ncorrec,ll 3873 save ncorrec 3874 data ncorrec/0/ 3875 c$OMP THREADPRIVATE(ncorrec) 3876 3877 c 3878 c----------------------------------------------------------------------- 3879 c initialisation: 3880 c --------------- 3881 c 3882 sorties=.true. 3883 IF(ngrid.NE.klon) THEN 3884 PRINT* 3885 PRINT*,'STOP dans convadj' 3886 PRINT*,'ngrid =',ngrid 3887 PRINT*,'klon =',klon 3888 ENDIF 3889 c 3890 c----------------------------------------------------------------------- 3891 c incrementation eventuelle de tendances precedentes: 3892 c --------------------------------------------------- 3893 3894 print*,'0 OK convect8' 3895 3896 DO 1010 l=1,nlay 3897 DO 1015 ig=1,ngrid 3898 zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA 3899 zh(ig,l)=pt(ig,l)/zpspsk(ig,l) 3900 zu(ig,l)=pu(ig,l) 3901 zv(ig,l)=pv(ig,l) 3902 zo(ig,l)=po(ig,l) 3903 ztv(ig,l)=zh(ig,l)*(1.+0.61*zo(ig,l)) 3904 1015 CONTINUE 3905 1010 CONTINUE 3906 3907 print*,'1 OK convect8' 3908 c -------------------- 3909 c 3910 c 3911 c + + + + + + + + + + + 3912 c 3913 c 3914 c wa, fraca, wd, fracd -------------------- zlev(2), rhobarz 3915 c wh,wt,wo ... 3916 c 3917 c + + + + + + + + + + + zh,zu,zv,zo,rho 3918 c 3919 c 3920 c -------------------- zlev(1) 3921 c \\\\\\\\\\\\\\\\\\\\ 3922 c 3923 c 3924 3925 c----------------------------------------------------------------------- 3926 c Calcul des altitudes des couches 3927 c----------------------------------------------------------------------- 3928 3929 do l=2,nlay 3930 do ig=1,ngrid 3931 zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG 3932 enddo 3933 enddo 3934 do ig=1,ngrid 3935 zlev(ig,1)=0. 3936 zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG 3937 enddo 3938 do l=1,nlay 3939 do ig=1,ngrid 3940 zlay(ig,l)=pphi(ig,l)/RG 3941 enddo 3942 enddo 3943 3944 c print*,'2 OK convect8' 3945 c----------------------------------------------------------------------- 3946 c Calcul des densites 3947 c----------------------------------------------------------------------- 3948 3949 do l=1,nlay 3950 do ig=1,ngrid 3951 rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)*RD*zh(ig,l)) 3952 enddo 3953 enddo 3954 3955 do l=2,nlay 3956 do ig=1,ngrid 3957 rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1)) 3958 enddo 3959 enddo 3960 3961 do k=1,nlay 3962 do l=1,nlay+1 3963 do ig=1,ngrid 3964 wa(ig,k,l)=0. 3965 enddo 3966 enddo 3967 enddo 3968 3969 c print*,'3 OK convect8' 3970 c------------------------------------------------------------------ 3971 c Calcul de w2, quarre de w a partir de la cape 3972 c a partir de w2, on calcule wa, vitesse de l'ascendance 3973 c 3974 c ATTENTION: Dans cette version, pour cause d'economie de memoire, 3975 c w2 est stoke dans wa 3976 c 3977 c ATTENTION: dans convect8, on n'utilise le calcule des wa 3978 c independants par couches que pour calculer l'entrainement 3979 c a la base et la hauteur max de l'ascendance. 3980 c 3981 c Indicages: 3982 c l'ascendance provenant du niveau k traverse l'interface l avec 3983 c une vitesse wa(k,l). 3984 c 3985 c -------------------- 3986 c 3987 c + + + + + + + + + + 3988 c 3989 c wa(k,l) ---- -------------------- l 3990 c /\ 3991 c /||\ + + + + + + + + + + 3992 c || 3993 c || -------------------- 3994 c || 3995 c || + + + + + + + + + + 3996 c || 3997 c || -------------------- 3998 c ||__ 3999 c |___ + + + + + + + + + + k 4000 c 4001 c -------------------- 4002 c 4003 c 4004 c 4005 c------------------------------------------------------------------ 4006 4007 cCR: ponderation entrainement des couches instables 4008 cdef des entr_star tels que entr=f*entr_star 4009 do l=1,klev 4010 do ig=1,ngrid 4011 entr_star(ig,l)=0. 4012 enddo 4013 enddo 4014 c determination de la longueur de la couche d entrainement 4015 do ig=1,ngrid 4016 lentr(ig)=1 4017 enddo 4018 4019 con ne considere que les premieres couches instables 4020 do k=nlay-2,1,-1 4021 do ig=1,ngrid 4022 if (ztv(ig,k).gt.ztv(ig,k+1).and. 4023 s ztv(ig,k+1).le.ztv(ig,k+2)) then 4024 lentr(ig)=k 4025 endif 4026 enddo 4027 enddo 4028 4029 c determination du lmin: couche d ou provient le thermique 4030 do ig=1,ngrid 4031 lmin(ig)=1 4032 enddo 4033 do ig=1,ngrid 4034 do l=nlay,2,-1 4035 if (ztv(ig,l-1).gt.ztv(ig,l)) then 4036 lmin(ig)=l-1 4037 endif 4038 enddo 4039 enddo 4040 c 4041 c definition de l'entrainement des couches 4042 do l=1,klev-1 4043 do ig=1,ngrid 4044 if (ztv(ig,l).gt.ztv(ig,l+1).and. 4045 s l.ge.lmin(ig).and.l.le.lentr(ig)) then 4046 entr_star(ig,l)=(ztv(ig,l)-ztv(ig,l+1))* 4047 s (zlev(ig,l+1)-zlev(ig,l)) 4048 endif 4049 enddo 4050 enddo 4051 c pas de thermique si couches 1->5 stables 4052 do ig=1,ngrid 4053 if (lmin(ig).gt.5) then 4054 do l=1,klev 4055 entr_star(ig,l)=0. 4056 enddo 4057 endif 4058 enddo 4059 c calcul de l entrainement total 4060 do ig=1,ngrid 4061 entr_star_tot(ig)=0. 4062 enddo 4063 do ig=1,ngrid 4064 do k=1,klev 4065 entr_star_tot(ig)=entr_star_tot(ig)+entr_star(ig,k) 4066 enddo 4067 enddo 4068 c 4069 print*,'fin calcul entr_star' 4070 do k=1,klev 4071 do ig=1,ngrid 4072 ztva(ig,k)=ztv(ig,k) 4073 enddo 4074 enddo 4075 cRC 4076 c print*,'7 OK convect8' 4077 do k=1,klev+1 4078 do ig=1,ngrid 4079 zw2(ig,k)=0. 4080 fmc(ig,k)=0. 4081 cCR 4082 f_star(ig,k)=0. 4083 cRC 4084 larg_cons(ig,k)=0. 4085 larg_detr(ig,k)=0. 4086 wa_moy(ig,k)=0. 4087 enddo 4088 enddo 4089 4090 c print*,'8 OK convect8' 4091 do ig=1,ngrid 4092 linter(ig)=1. 4093 lmaxa(ig)=1 4094 lmix(ig)=1 4095 wmaxa(ig)=0. 4096 enddo 4097 4098 cCR: 4099 do l=1,nlay-2 4100 do ig=1,ngrid 4101 if (ztv(ig,l).gt.ztv(ig,l+1) 4102 s .and.entr_star(ig,l).gt.1.e-10 4103 s .and.zw2(ig,l).lt.1e-10) then 4104 f_star(ig,l+1)=entr_star(ig,l) 4105 ctest:calcul de dteta 4106 zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1) 4107 s *(zlev(ig,l+1)-zlev(ig,l)) 4108 s *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l)) 4109 larg_detr(ig,l)=0. 4110 else if ((zw2(ig,l).ge.1e-10).and. 4111 s (f_star(ig,l)+entr_star(ig,l).gt.1.e-10)) then 4112 f_star(ig,l+1)=f_star(ig,l)+entr_star(ig,l) 4113 ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l) 4114 s *ztv(ig,l))/f_star(ig,l+1) 4115 zw2(ig,l+1)=zw2(ig,l)*(f_star(ig,l)/f_star(ig,l+1))**2+ 4116 s 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) 4117 s *(zlev(ig,l+1)-zlev(ig,l)) 4118 endif 4119 c determination de zmax continu par interpolation lineaire 4120 if (zw2(ig,l+1).lt.0.) then 4121 ctest 4122 if (abs(zw2(ig,l+1)-zw2(ig,l)).lt.1e-10) then 4123 print*,'pb linter' 4124 endif 4125 linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) 4126 s -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l)) 4127 zw2(ig,l+1)=0. 4128 lmaxa(ig)=l 4129 else 4130 if (zw2(ig,l+1).lt.0.) then 4131 print*,'pb1 zw2<0' 4132 endif 4133 wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) 4134 endif 4135 if (wa_moy(ig,l+1).gt.wmaxa(ig)) then 4136 c lmix est le niveau de la couche ou w (wa_moy) est maximum 4137 lmix(ig)=l+1 4138 wmaxa(ig)=wa_moy(ig,l+1) 4139 endif 4140 enddo 4141 enddo 4142 print*,'fin calcul zw2' 4143 c 4144 c Calcul de la couche correspondant a la hauteur du thermique 4145 do ig=1,ngrid 4146 lmax(ig)=lentr(ig) 4147 enddo 4148 do ig=1,ngrid 4149 do l=nlay,lentr(ig)+1,-1 4150 if (zw2(ig,l).le.1.e-10) then 4151 lmax(ig)=l-1 4152 endif 4153 enddo 4154 enddo 4155 c pas de thermique si couches 1->5 stables 4156 do ig=1,ngrid 4157 if (lmin(ig).gt.5) then 4158 lmax(ig)=1 4159 lmin(ig)=1 4160 endif 4161 enddo 4162 c 4163 c Determination de zw2 max 4164 do ig=1,ngrid 4165 wmax(ig)=0. 4166 enddo 4167 4168 do l=1,nlay 4169 do ig=1,ngrid 4170 if (l.le.lmax(ig)) then 4171 if (zw2(ig,l).lt.0.)then 4172 print*,'pb2 zw2<0' 4173 endif 4174 zw2(ig,l)=sqrt(zw2(ig,l)) 4175 wmax(ig)=max(wmax(ig),zw2(ig,l)) 4176 else 4177 zw2(ig,l)=0. 4178 endif 4179 enddo 4180 enddo 4181 4182 c Longueur caracteristique correspondant a la hauteur des thermiques. 4183 do ig=1,ngrid 4184 zmax(ig)=0. 4185 zlevinter(ig)=zlev(ig,1) 4186 enddo 4187 do ig=1,ngrid 4188 c calcul de zlevinter 4189 zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))* 4190 s linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1) 4191 s -zlev(ig,lmax(ig))) 4192 zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig))) 4193 enddo 4194 4195 print*,'avant fermeture' 4196 c Fermeture,determination de f 4197 do ig=1,ngrid 4198 entr_star2(ig)=0. 4199 enddo 4200 do ig=1,ngrid 4201 if (entr_star_tot(ig).LT.1.e-10) then 4202 f(ig)=0. 4203 else 4204 do k=lmin(ig),lentr(ig) 4205 entr_star2(ig)=entr_star2(ig)+entr_star(ig,k)**2 4206 s /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))) 4207 enddo 4208 c Nouvelle fermeture 4209 f(ig)=wmax(ig)/(max(500.,zmax(ig))*r_aspect 4210 s *entr_star2(ig))*entr_star_tot(ig) 4211 ctest 4212 c if (first) then 4213 c f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig) 4214 c s *wmax(ig)) 4215 c endif 4216 endif 4217 c f0(ig)=f(ig) 4218 c first=.true. 4219 enddo 4220 print*,'apres fermeture' 4221 4222 c Calcul de l'entrainement 4223 do k=1,klev 4224 do ig=1,ngrid 4225 entr(ig,k)=f(ig)*entr_star(ig,k) 4226 enddo 4227 enddo 4228 c Calcul des flux 4229 do ig=1,ngrid 4230 do l=1,lmax(ig)-1 4231 fmc(ig,l+1)=fmc(ig,l)+entr(ig,l) 4232 enddo 4233 enddo 4234 4235 cRC 4236 4237 4238 c print*,'9 OK convect8' 4239 c print*,'WA1 ',wa_moy 4240 4241 c determination de l'indice du debut de la mixed layer ou w decroit 4242 4243 c calcul de la largeur de chaque ascendance dans le cas conservatif. 4244 c dans ce cas simple, on suppose que la largeur de l'ascendance provenant 4245 c d'une couche est égale à la hauteur de la couche alimentante. 4246 c La vitesse maximale dans l'ascendance est aussi prise comme estimation 4247 c de la vitesse d'entrainement horizontal dans la couche alimentante. 4248 4249 do l=2,nlay 4250 do ig=1,ngrid 4251 if (l.le.lmaxa(ig)) then 4252 zw=max(wa_moy(ig,l),1.e-10) 4253 larg_cons(ig,l)=zmax(ig)*r_aspect 4254 s *fmc(ig,l)/(rhobarz(ig,l)*zw) 4255 endif 4256 enddo 4257 enddo 4258 4259 do l=2,nlay 4260 do ig=1,ngrid 4261 if (l.le.lmaxa(ig)) then 4262 c if (idetr.eq.0) then 4263 c cette option est finalement en dur. 4264 if ((l_mix*zlev(ig,l)).lt.0.)then 4265 print*,'pb l_mix*zlev<0' 4266 endif 4267 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l)) 4268 c else if (idetr.eq.1) then 4269 c larg_detr(ig,l)=larg_cons(ig,l) 4270 c s *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig)) 4271 c else if (idetr.eq.2) then 4272 c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l)) 4273 c s *sqrt(wa_moy(ig,l)) 4274 c else if (idetr.eq.4) then 4275 c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l)) 4276 c s *wa_moy(ig,l) 4277 c endif 4278 endif 4279 enddo 4280 enddo 4281 4282 c print*,'10 OK convect8' 4283 c print*,'WA2 ',wa_moy 4284 c calcul de la fraction de la maille concernée par l'ascendance en tenant 4285 c compte de l'epluchage du thermique. 4286 c 4287 cCR def de zmix continu (profil parabolique des vitesses) 4288 do ig=1,ngrid 4289 if (lmix(ig).gt.1.) then 4290 c test 4291 if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig))) 4292 s *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1))) 4293 s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1)) 4294 s *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10) 4295 s then 4296 c 4297 zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig))) 4298 s *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2) 4299 s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1)) 4300 s *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2)) 4301 s /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig))) 4302 s *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1))) 4303 s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1)) 4304 s *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig)))))) 4305 else 4306 zmix(ig)=zlev(ig,lmix(ig)) 4307 print*,'pb zmix' 4308 endif 4309 else 4310 zmix(ig)=0. 4311 endif 4312 ctest 4313 if ((zmax(ig)-zmix(ig)).lt.0.) then 4314 zmix(ig)=0.99*zmax(ig) 4315 c print*,'pb zmix>zmax' 4316 endif 4317 enddo 4318 c 4319 c calcul du nouveau lmix correspondant 4320 do ig=1,ngrid 4321 do l=1,klev 4322 if (zmix(ig).ge.zlev(ig,l).and. 4323 s zmix(ig).lt.zlev(ig,l+1)) then 4324 lmix(ig)=l 4325 endif 4326 enddo 4327 enddo 4328 c 4329 do l=2,nlay 4330 do ig=1,ngrid 4331 if(larg_cons(ig,l).gt.1.) then 4332 c print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK' 4333 fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l)) 4334 s /(r_aspect*zmax(ig)) 4335 c test 4336 fraca(ig,l)=max(fraca(ig,l),0.) 4337 fraca(ig,l)=min(fraca(ig,l),0.5) 4338 fracd(ig,l)=1.-fraca(ig,l) 4339 fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig)) 4340 else 4341 c wa_moy(ig,l)=0. 4342 fraca(ig,l)=0. 4343 fracc(ig,l)=0. 4344 fracd(ig,l)=1. 4345 endif 4346 enddo 4347 enddo 4348 cCR: calcul de fracazmix 4349 do ig=1,ngrid 4350 fracazmix(ig)=(fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/ 4351 s (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig) 4352 s +fraca(ig,lmix(ig))-zlev(ig,lmix(ig))*(fraca(ig,lmix(ig)+1) 4353 s -fraca(ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig))) 4354 enddo 4355 c 4356 do l=2,nlay 4357 do ig=1,ngrid 4358 if(larg_cons(ig,l).gt.1.) then 4359 if (l.gt.lmix(ig)) then 4360 ctest 4361 if (zmax(ig)-zmix(ig).lt.1.e-10) then 4362 c print*,'pb xxx' 4363 xxx(ig,l)=(lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig)) 4364 else 4365 xxx(ig,l)=(zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig)) 4366 endif 4367 if (idetr.eq.0) then 4368 fraca(ig,l)=fracazmix(ig) 4369 else if (idetr.eq.1) then 4370 fraca(ig,l)=fracazmix(ig)*xxx(ig,l) 4371 else if (idetr.eq.2) then 4372 fraca(ig,l)=fracazmix(ig)*(1.-(1.-xxx(ig,l))**2) 4373 else 4374 fraca(ig,l)=fracazmix(ig)*xxx(ig,l)**2 4375 endif 4376 c print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL' 4377 fraca(ig,l)=max(fraca(ig,l),0.) 4378 fraca(ig,l)=min(fraca(ig,l),0.5) 4379 fracd(ig,l)=1.-fraca(ig,l) 4380 fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig)) 4381 endif 4382 endif 4383 enddo 4384 enddo 4385 4386 print*,'fin calcul fraca' 4387 c print*,'11 OK convect8' 4388 c print*,'Ea3 ',wa_moy 4389 c------------------------------------------------------------------ 4390 c Calcul de fracd, wd 4391 c somme wa - wd = 0 4392 c------------------------------------------------------------------ 4393 4394 4395 do ig=1,ngrid 4396 fm(ig,1)=0. 4397 fm(ig,nlay+1)=0. 4398 enddo 4399 4400 do l=2,nlay 4401 do ig=1,ngrid 4402 fm(ig,l)=fraca(ig,l)*wa_moy(ig,l)*rhobarz(ig,l) 4403 cCR:test 4404 if (entr(ig,l-1).lt.1e-10.and.fm(ig,l).gt.fm(ig,l-1) 4405 s .and.l.gt.lmix(ig)) then 4406 fm(ig,l)=fm(ig,l-1) 4407 c write(1,*)'ajustement fm, l',l 4408 endif 4409 c write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l) 4410 cRC 4411 enddo 4412 do ig=1,ngrid 4413 if(fracd(ig,l).lt.0.1) then 4414 abort_message = 'fracd trop petit' 4415 CALL abort_gcm (modname,abort_message,1) 4416 else 4417 c vitesse descendante "diagnostique" 4418 wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l)) 4419 endif 4420 enddo 4421 enddo 4422 4423 do l=1,nlay 4424 do ig=1,ngrid 4425 c masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) 4426 masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG 4427 enddo 4428 enddo 4429 4430 print*,'12 OK convect8' 4431 c print*,'WA4 ',wa_moy 4432 cc------------------------------------------------------------------ 4433 c calcul du transport vertical 4434 c------------------------------------------------------------------ 4435 4436 go to 4444 4437 c print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep 4438 do l=2,nlay-1 4439 do ig=1,ngrid 4440 if(fm(ig,l+1)*ptimestep.gt.masse(ig,l) 4441 s .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then 4442 c print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM=' 4443 c s ,fm(ig,l+1)*ptimestep 4444 c s ,' M=',masse(ig,l),masse(ig,l+1) 4445 endif 4446 enddo 4447 enddo 4448 4449 do l=1,nlay 4450 do ig=1,ngrid 4451 if(entr(ig,l)*ptimestep.gt.masse(ig,l)) then 4452 c print*,'WARN!!! E>M ig=',ig,' l=',l,' E==' 4453 c s ,entr(ig,l)*ptimestep 4454 c s ,' M=',masse(ig,l) 4455 endif 4456 enddo 4457 enddo 4458 4459 do l=1,nlay 4460 do ig=1,ngrid 4461 if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then 4462 c print*,'WARN!!! fm exagere ig=',ig,' l=',l 4463 c s ,' FM=',fm(ig,l) 4464 endif 4465 if(.not.masse(ig,l).ge.1.e-10 4466 s .or..not.masse(ig,l).le.1.e4) then 4467 c print*,'WARN!!! masse exagere ig=',ig,' l=',l 4468 c s ,' M=',masse(ig,l) 4469 c print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)', 4470 c s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l) 4471 c print*,'zlev(ig,l+1),zlev(ig,l)' 4472 c s ,zlev(ig,l+1),zlev(ig,l) 4473 c print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)' 4474 c s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1) 4475 endif 4476 if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then 4477 c print*,'WARN!!! entr exagere ig=',ig,' l=',l 4478 c s ,' E=',entr(ig,l) 4479 endif 4480 enddo 4481 enddo 4482 4483 4444 continue 4484 4485 cCR:redefinition du entr 4486 do l=1,nlay 4487 do ig=1,ngrid 4488 detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1) 4489 if (detr(ig,l).lt.0.) then 4490 entr(ig,l)=entr(ig,l)-detr(ig,l) 4491 detr(ig,l)=0. 4492 c print*,'WARNING !!! detrainement negatif ',ig,l 4493 endif 4494 enddo 4495 enddo 4496 cRC 4497 if (w2di.eq.1) then 4498 fm0=fm0+ptimestep*(fm-fm0)/tho 4499 entr0=entr0+ptimestep*(entr-entr0)/tho 4500 else 4501 fm0=fm 4502 entr0=entr 4503 endif 4504 4505 if (1.eq.1) then 4506 call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse 4507 . ,zh,zdhadj,zha) 4508 call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse 4509 . ,zo,pdoadj,zoa) 4510 else 4511 call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca 4512 . ,zh,zdhadj,zha) 4513 call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca 4514 . ,zo,pdoadj,zoa) 4515 endif 4516 4517 if (1.eq.0) then 4518 call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse 4519 . ,fraca,zmax 4520 . ,zu,zv,pduadj,pdvadj,zua,zva) 4521 else 4522 call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse 4523 . ,zu,pduadj,zua) 4524 call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse 4525 . ,zv,pdvadj,zva) 4526 endif 4527 4528 do l=1,nlay 4529 do ig=1,ngrid 4530 zf=0.5*(fracc(ig,l)+fracc(ig,l+1)) 4531 zf2=zf/(1.-zf) 4532 thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l))**2 4533 wth2(ig,l)=zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2 4534 enddo 4535 enddo 4536 4537 4538 4539 c print*,'13 OK convect8' 4540 c print*,'WA5 ',wa_moy 4541 do l=1,nlay 4542 do ig=1,ngrid 4543 pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l) 4544 enddo 4545 enddo 4546 4547 4548 c do l=1,nlay 4549 c do ig=1,ngrid 4550 c if(abs(pdtadj(ig,l))*86400..gt.500.) then 4551 c print*,'WARN!!! ig=',ig,' l=',l 4552 c s ,' pdtadj=',pdtadj(ig,l) 4553 c endif 4554 c if(abs(pdoadj(ig,l))*86400..gt.1.) then 4555 c print*,'WARN!!! ig=',ig,' l=',l 4556 c s ,' pdoadj=',pdoadj(ig,l) 4557 c endif 4558 c enddo 4559 c enddo 4560 4561 print*,'14 OK convect8' 4562 c------------------------------------------------------------------ 4563 c Calculs pour les sorties 4564 c------------------------------------------------------------------ 4565 4566 if(sorties) then 4567 do l=1,nlay 4568 do ig=1,ngrid 4569 zla(ig,l)=(1.-fracd(ig,l))*zmax(ig) 4570 zld(ig,l)=fracd(ig,l)*zmax(ig) 4571 if(1.-fracd(ig,l).gt.1.e-10) 4572 s zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l)) 4573 enddo 4574 enddo 4575 4576 cdeja fait 4577 c do l=1,nlay 4578 c do ig=1,ngrid 4579 c detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1) 4580 c if (detr(ig,l).lt.0.) then 4581 c entr(ig,l)=entr(ig,l)-detr(ig,l) 4582 c detr(ig,l)=0. 4583 c print*,'WARNING !!! detrainement negatif ',ig,l 4584 c endif 4585 c enddo 4586 c enddo 4587 4588 c print*,'15 OK convect8' 4589 4590 isplit=isplit+1 4591 4592 4593 c #define und 4594 goto 123 4595 #ifdef und 4596 CALL writeg1d(1,nlay,wd,'wd ','wd ') 4597 CALL writeg1d(1,nlay,zwa,'wa ','wa ') 4598 CALL writeg1d(1,nlay,fracd,'fracd ','fracd ') 4599 CALL writeg1d(1,nlay,fraca,'fraca ','fraca ') 4600 CALL writeg1d(1,nlay,wa_moy,'wam ','wam ') 4601 CALL writeg1d(1,nlay,zla,'la ','la ') 4602 CALL writeg1d(1,nlay,zld,'ld ','ld ') 4603 CALL writeg1d(1,nlay,pt,'pt ','pt ') 4604 CALL writeg1d(1,nlay,zh,'zh ','zh ') 4605 CALL writeg1d(1,nlay,zha,'zha ','zha ') 4606 CALL writeg1d(1,nlay,zu,'zu ','zu ') 4607 CALL writeg1d(1,nlay,zv,'zv ','zv ') 4608 CALL writeg1d(1,nlay,zo,'zo ','zo ') 4609 CALL writeg1d(1,nlay,wh,'wh ','wh ') 4610 CALL writeg1d(1,nlay,wu,'wu ','wu ') 4611 CALL writeg1d(1,nlay,wv,'wv ','wv ') 4612 CALL writeg1d(1,nlay,wo,'w15uo ','wXo ') 4613 CALL writeg1d(1,nlay,zdhadj,'zdhadj ','zdhadj ') 4614 CALL writeg1d(1,nlay,pduadj,'pduadj ','pduadj ') 4615 CALL writeg1d(1,nlay,pdvadj,'pdvadj ','pdvadj ') 4616 CALL writeg1d(1,nlay,pdoadj,'pdoadj ','pdoadj ') 4617 CALL writeg1d(1,nlay,entr ,'entr ','entr ') 4618 CALL writeg1d(1,nlay,detr ,'detr ','detr ') 4619 CALL writeg1d(1,nlay,fm ,'fm ','fm ') 4620 4621 CALL writeg1d(1,nlay,pdtadj,'pdtadj ','pdtadj ') 4622 CALL writeg1d(1,nlay,pplay,'pplay ','pplay ') 4623 CALL writeg1d(1,nlay,pplev,'pplev ','pplev ') 4624 4625 c recalcul des flux en diagnostique... 4626 c print*,'PAS DE TEMPS ',ptimestep 4627 call dt2F(pplev,pplay,pt,pdtadj,wh) 4628 CALL writeg1d(1,nlay,wh,'wh2 ','wh2 ') 4629 #endif 4630 123 continue 4631 #define troisD 4632 #ifdef troisD 4633 c if (sorties) then 4634 print*,'Debut des wrgradsfi' 4635 4636 c print*,'16 OK convect8' 4637 call wrgradsfi(1,nlay,wd,'wd ','wd ') 4638 call wrgradsfi(1,nlay,zwa,'wa ','wa ') 4639 call wrgradsfi(1,nlay,fracd,'fracd ','fracd ') 4640 call wrgradsfi(1,nlay,fraca,'fraca ','fraca ') 4641 call wrgradsfi(1,nlay,xxx,'xxx ','xxx ') 4642 call wrgradsfi(1,nlay,wa_moy,'wam ','wam ') 4643 c print*,'WA6 ',wa_moy 4644 call wrgradsfi(1,nlay,zla,'la ','la ') 4645 call wrgradsfi(1,nlay,zld,'ld ','ld ') 4646 call wrgradsfi(1,nlay,pt,'pt ','pt ') 4647 call wrgradsfi(1,nlay,zh,'zh ','zh ') 4648 call wrgradsfi(1,nlay,zha,'zha ','zha ') 4649 call wrgradsfi(1,nlay,zua,'zua ','zua ') 4650 call wrgradsfi(1,nlay,zva,'zva ','zva ') 4651 call wrgradsfi(1,nlay,zu,'zu ','zu ') 4652 call wrgradsfi(1,nlay,zv,'zv ','zv ') 4653 call wrgradsfi(1,nlay,zo,'zo ','zo ') 4654 call wrgradsfi(1,nlay,wh,'wh ','wh ') 4655 call wrgradsfi(1,nlay,wu,'wu ','wu ') 4656 call wrgradsfi(1,nlay,wv,'wv ','wv ') 4657 call wrgradsfi(1,nlay,wo,'wo ','wo ') 4658 call wrgradsfi(1,1,zmax,'zmax ','zmax ') 4659 call wrgradsfi(1,nlay,zdhadj,'zdhadj ','zdhadj ') 4660 call wrgradsfi(1,nlay,pduadj,'pduadj ','pduadj ') 4661 call wrgradsfi(1,nlay,pdvadj,'pdvadj ','pdvadj ') 4662 call wrgradsfi(1,nlay,pdoadj,'pdoadj ','pdoadj ') 4663 call wrgradsfi(1,nlay,entr,'entr ','entr ') 4664 call wrgradsfi(1,nlay,detr,'detr ','detr ') 4665 call wrgradsfi(1,nlay,fm,'fm ','fm ') 4666 call wrgradsfi(1,nlay,fmc,'fmc ','fmc ') 4667 call wrgradsfi(1,nlay,zw2,'zw2 ','zw2 ') 4668 call wrgradsfi(1,nlay,ztva,'ztva ','ztva ') 4669 call wrgradsfi(1,nlay,ztv,'ztv ','ztv ') 4670 4671 call wrgradsfi(1,nlay,zo,'zo ','zo ') 4672 call wrgradsfi(1,nlay,larg_cons,'Lc ','Lc ') 4673 call wrgradsfi(1,nlay,larg_detr,'Ldetr ','Ldetr ') 4674 4675 cCR:nouveaux diagnostiques 4676 call wrgradsfi(1,nlay,entr_star ,'entr_star ','entr_star ') 4677 call wrgradsfi(1,nlay,f_star ,'f_star ','f_star ') 4678 call wrgradsfi(1,1,zmax,'zmax ','zmax ') 4679 call wrgradsfi(1,1,zmix,'zmix ','zmix ') 4680 zsortie1d(:)=lmax(:) 4681 call wrgradsfi(1,1,zsortie1d,'lmax ','lmax ') 4682 call wrgradsfi(1,1,wmax,'wmax ','wmax ') 4683 zsortie1d(:)=lmix(:) 4684 call wrgradsfi(1,1,zsortie1d,'lmix ','lmix ') 4685 zsortie1d(:)=lentr(:) 4686 call wrgradsfi(1,1,zsortie1d,'lentr ','lentr ') 4687 4688 c print*,'17 OK convect8' 4689 4690 do k=1,klev/10 4691 write(str2,'(i2.2)') k 4692 str10='wa'//str2 4693 do l=1,nlay 4694 do ig=1,ngrid 4695 zsortie(ig,l)=wa(ig,k,l) 4696 enddo 4697 enddo 4698 CALL wrgradsfi(1,nlay,zsortie,str10,str10) 4699 do l=1,nlay 4700 do ig=1,ngrid 4701 zsortie(ig,l)=larg_part(ig,k,l) 4702 enddo 4703 enddo 4704 str10='la'//str2 4705 CALL wrgradsfi(1,nlay,zsortie,str10,str10) 4706 enddo 4707 4708 4709 c print*,'18 OK convect8' 4710 c endif 4711 print*,'Fin des wrgradsfi' 4712 #endif 4713 4714 endif 4715 4716 c if(wa_moy(1,4).gt.1.e-10) stop 4717 4718 print*,'19 OK convect8' 4291 ! print*,'19 OK convect8' 4719 4292 return 4720 4293 end … … 5239 4812 5240 4813 real count_time 5241 integer isplit,nsplit,ialt 5242 parameter (nsplit=10) 5243 data isplit/0/ 5244 save isplit 5245 c$OMP THREADPRIVATE(isplit) 4814 integer ialt 5246 4815 5247 4816 logical sorties … … 6007 5576 enddo 6008 5577 6009 cdeja fait 6010 c do l=1,nlay 6011 c do ig=1,ngrid 6012 c detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1) 6013 c if (detr(ig,l).lt.0.) then 6014 c entr(ig,l)=entr(ig,l)-detr(ig,l) 6015 c detr(ig,l)=0. 6016 c print*,'WARNING !!! detrainement negatif ',ig,l 6017 c endif 6018 c enddo 6019 c enddo 6020 6021 c print*,'15 OK convect8' 6022 6023 isplit=isplit+1 6024 6025 6026 c #define und 6027 goto 123 6028 #ifdef und 6029 CALL writeg1d(1,nlay,wd,'wd ','wd ') 6030 CALL writeg1d(1,nlay,zwa,'wa ','wa ') 6031 CALL writeg1d(1,nlay,fracd,'fracd ','fracd ') 6032 CALL writeg1d(1,nlay,fraca,'fraca ','fraca ') 6033 CALL writeg1d(1,nlay,wa_moy,'wam ','wam ') 6034 CALL writeg1d(1,nlay,zla,'la ','la ') 6035 CALL writeg1d(1,nlay,zld,'ld ','ld ') 6036 CALL writeg1d(1,nlay,pt,'pt ','pt ') 6037 CALL writeg1d(1,nlay,zh,'zh ','zh ') 6038 CALL writeg1d(1,nlay,zha,'zha ','zha ') 6039 CALL writeg1d(1,nlay,zu,'zu ','zu ') 6040 CALL writeg1d(1,nlay,zv,'zv ','zv ') 6041 CALL writeg1d(1,nlay,zo,'zo ','zo ') 6042 CALL writeg1d(1,nlay,wh,'wh ','wh ') 6043 CALL writeg1d(1,nlay,wu,'wu ','wu ') 6044 CALL writeg1d(1,nlay,wv,'wv ','wv ') 6045 CALL writeg1d(1,nlay,wo,'w15uo ','wXo ') 6046 CALL writeg1d(1,nlay,zdhadj,'zdhadj ','zdhadj ') 6047 CALL writeg1d(1,nlay,pduadj,'pduadj ','pduadj ') 6048 CALL writeg1d(1,nlay,pdvadj,'pdvadj ','pdvadj ') 6049 CALL writeg1d(1,nlay,pdoadj,'pdoadj ','pdoadj ') 6050 CALL writeg1d(1,nlay,entr ,'entr ','entr ') 6051 CALL writeg1d(1,nlay,detr ,'detr ','detr ') 6052 CALL writeg1d(1,nlay,fm ,'fm ','fm ') 6053 6054 CALL writeg1d(1,nlay,pdtadj,'pdtadj ','pdtadj ') 6055 CALL writeg1d(1,nlay,pplay,'pplay ','pplay ') 6056 CALL writeg1d(1,nlay,pplev,'pplev ','pplev ') 6057 6058 c recalcul des flux en diagnostique... 6059 c print*,'PAS DE TEMPS ',ptimestep 6060 call dt2F(pplev,pplay,pt,pdtadj,wh) 6061 CALL writeg1d(1,nlay,wh,'wh2 ','wh2 ') 6062 #endif 6063 123 continue 6064 ! #define troisD 6065 #ifdef troisD 6066 c if (sorties) then 6067 print*,'Debut des wrgradsfi' 6068 6069 c print*,'16 OK convect8' 6070 call wrgradsfi(1,nlay,wd,'wd ','wd ') 6071 call wrgradsfi(1,nlay,zwa,'wa ','wa ') 6072 call wrgradsfi(1,nlay,fracd,'fracd ','fracd ') 6073 call wrgradsfi(1,nlay,fraca,'fraca ','fraca ') 6074 call wrgradsfi(1,nlay,xxx,'xxx ','xxx ') 6075 call wrgradsfi(1,nlay,wa_moy,'wam ','wam ') 6076 c print*,'WA6 ',wa_moy 6077 call wrgradsfi(1,nlay,zla,'la ','la ') 6078 call wrgradsfi(1,nlay,zld,'ld ','ld ') 6079 call wrgradsfi(1,nlay,pt,'pt ','pt ') 6080 call wrgradsfi(1,nlay,zh,'zh ','zh ') 6081 call wrgradsfi(1,nlay,zha,'zha ','zha ') 6082 call wrgradsfi(1,nlay,zua,'zua ','zua ') 6083 call wrgradsfi(1,nlay,zva,'zva ','zva ') 6084 call wrgradsfi(1,nlay,zu,'zu ','zu ') 6085 call wrgradsfi(1,nlay,zv,'zv ','zv ') 6086 call wrgradsfi(1,nlay,zo,'zo ','zo ') 6087 call wrgradsfi(1,nlay,wh,'wh ','wh ') 6088 call wrgradsfi(1,nlay,wu,'wu ','wu ') 6089 call wrgradsfi(1,nlay,wv,'wv ','wv ') 6090 call wrgradsfi(1,nlay,wo,'wo ','wo ') 6091 call wrgradsfi(1,1,zmax,'zmax ','zmax ') 6092 call wrgradsfi(1,nlay,zdhadj,'zdhadj ','zdhadj ') 6093 call wrgradsfi(1,nlay,pduadj,'pduadj ','pduadj ') 6094 call wrgradsfi(1,nlay,pdvadj,'pdvadj ','pdvadj ') 6095 call wrgradsfi(1,nlay,pdoadj,'pdoadj ','pdoadj ') 6096 call wrgradsfi(1,nlay,entr,'entr ','entr ') 6097 call wrgradsfi(1,nlay,detr,'detr ','detr ') 6098 call wrgradsfi(1,nlay,fm,'fm ','fm ') 6099 call wrgradsfi(1,nlay,fmc,'fmc ','fmc ') 6100 call wrgradsfi(1,nlay,zw2,'zw2 ','zw2 ') 6101 call wrgradsfi(1,nlay,ztva,'ztva ','ztva ') 6102 call wrgradsfi(1,nlay,ztv,'ztv ','ztv ') 6103 6104 call wrgradsfi(1,nlay,zo,'zo ','zo ') 6105 call wrgradsfi(1,nlay,larg_cons,'Lc ','Lc ') 6106 call wrgradsfi(1,nlay,larg_detr,'Ldetr ','Ldetr ') 6107 6108 cCR:nouveaux diagnostiques 6109 call wrgradsfi(1,nlay,entr_star ,'entr_star ','entr_star ') 6110 call wrgradsfi(1,nlay,f_star ,'f_star ','f_star ') 6111 call wrgradsfi(1,1,zmax,'zmax ','zmax ') 6112 call wrgradsfi(1,1,zmix,'zmix ','zmix ') 6113 zsortie1d(:)=lmax(:) 6114 call wrgradsfi(1,1,zsortie1d,'lmax ','lmax ') 6115 call wrgradsfi(1,1,wmax,'wmax ','wmax ') 6116 zsortie1d(:)=lmix(:) 6117 call wrgradsfi(1,1,zsortie1d,'lmix ','lmix ') 6118 zsortie1d(:)=lentr(:) 6119 call wrgradsfi(1,1,zsortie1d,'lentr ','lentr ') 6120 6121 c print*,'17 OK convect8' 5578 6122 5579 6123 5580 do k=1,klev/10 … … 6139 5596 enddo 6140 5597 6141 6142 c print*,'18 OK convect8'6143 c endif6144 print*,'Fin des wrgradsfi'6145 #endif6146 6147 5598 endif 6148 5599 6149 c if(wa_moy(1,4).gt.1.e-10) stop 6150 6151 c print*,'19 OK convect8' 5600 5601 6152 5602 return 6153 5603 end
Note: See TracChangeset
for help on using the changeset viewer.