Index: trunk/LMDZ.MARS/libf/phymars/aeropacity.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/aeropacity.F (revision 517) +++ trunk/LMDZ.MARS/libf/phymars/aeropacity.F (revision 520) @@ -1,4 +1,4 @@ SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls, - & pq,tauscaling,tauref,tau,aerosol,reffrad,nueffrad, + & pq,tauscaling,tauref,tau,taucloudtes,aerosol,reffrad,nueffrad, & QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d) @@ -386,19 +386,19 @@ ENDDO ENDDO -c 3. Outputs - IF (ngrid.NE.1) THEN - CALL WRITEDIAGFI(ngridmx,'tauVIS','tauext VIS refwvl', - & ' ',2,taucloudvis) - CALL WRITEDIAGFI(ngridmx,'tauTES','tauabs IR refwvl', - & ' ',2,taucloudtes) - IF (callstats) THEN - CALL wstats(ngridmx,'tauVIS','tauext VIS refwvl', - & ' ',2,taucloudvis) - CALL wstats(ngridmx,'tauTES','tauabs IR refwvl', - & ' ',2,taucloudtes) - ENDIF - ELSE -c CALL writeg1d(ngrid,1,taucloudtes,'tautes','NU') - ENDIF +c 3. Outputs -- Now done in physiq.F +! IF (ngrid.NE.1) THEN +! CALL WRITEDIAGFI(ngridmx,'tauVIS','tauext VIS refwvl', +! & ' ',2,taucloudvis) +! CALL WRITEDIAGFI(ngridmx,'tauTES','tauabs IR refwvl', +! & ' ',2,taucloudtes) +! IF (callstats) THEN +! CALL wstats(ngridmx,'tauVIS','tauext VIS refwvl', +! & ' ',2,taucloudvis) +! CALL wstats(ngridmx,'tauTES','tauabs IR refwvl', +! & ' ',2,taucloudtes) +! ENDIF +! ELSE +! CALL writeg1d(ngrid,1,taucloudtes,'tautes','NU') +! ENDIF c================================================================== END SELECT aerkind Index: trunk/LMDZ.MARS/libf/phymars/callradite.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/callradite.F (revision 517) +++ trunk/LMDZ.MARS/libf/phymars/callradite.F (revision 520) @@ -2,5 +2,6 @@ $ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, $ dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw, - & tauref,tau,aerosol,tauscaling,rdust,rice,nuice,co2ice) + & tauref,tau,aerosol,tauscaling,taucloudtes,rdust,rice, + & nuice,co2ice) IMPLICIT NONE @@ -174,4 +175,8 @@ REAL tauref(ngrid), tau(ngrid,naerkind) + REAL taucloudtes(ngridmx)! Cloud opacity at infrared + ! reference wavelength using + ! Qabs instead of Qext + ! (direct comparison with TES) REAL aerosol(ngrid,nlayer,naerkind) REAL rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m) @@ -395,6 +400,6 @@ c Computing aerosol optical depth in each layer: CALL aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls, - & pq,tauscaling,tauref,tau,aerosol,reffrad,nueffrad, - & QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d) + & pq,tauscaling,tauref,tau,taucloudtes,aerosol,reffrad, + & nueffrad,QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d) c Starting loop on sub-domain Index: trunk/LMDZ.MARS/libf/phymars/callsedim.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/callsedim.F (revision 517) +++ trunk/LMDZ.MARS/libf/phymars/callsedim.F (revision 520) @@ -4,4 +4,6 @@ & pq, pdqfi, pdqsed,pdqs_sed,nq, & tau, tauscaling) +! to use 'getin' + USE ioipsl_getincom IMPLICIT NONE @@ -21,5 +23,5 @@ c declarations: c ------------- - + #include "dimensions.h" #include "dimphys.h" @@ -69,4 +71,6 @@ c Sedimentation radius of water ice real rsedcloud(ngridmx,nlayermx) + real beta ! correction for the shape of the ice particles (cf. newsedim) + save beta c Cloud density (kg.m-3) real rhocloud(ngridmx,nlayermx) @@ -118,4 +122,5 @@ SAVE idust_mass,idust_number SAVE iccn_mass,iccn_number + c Firstcall: @@ -201,10 +206,15 @@ ENDIF !of if (scavenging) - IF (water) THEN - write(*,*) "water_param nueff Sedimentation:", nuice_sed - IF (activice) THEN - write(*,*) "water_param nueff Radiative:", nuice_ref - ENDIF - ENDIF +! IF (water) THEN +! write(*,*) "correction for the shape of the ice particles ?" +! beta=0.75 ! default value +! call getin("ice_shape",beta) +! write(*,*) " ice_shape = ",beta +! +! write(*,*) "water_param nueff Sedimentation:", nuice_sed +! IF (activice) THEN +! write(*,*) "water_param nueff Radiative:", nuice_ref +! ENDIF +! ENDIF firstcall=.false. @@ -269,13 +279,13 @@ if ((doubleq.and. & ((iq.eq.idust_mass).or. - & (iq.eq.idust_number))).or. - & (scavenging.and. - & ((iq.eq.iccn_mass).or. - & (iq.eq.iccn_number)))) then + & (iq.eq.idust_number)))) then !.or. +! & (scavenging.and. +! & ((iq.eq.iccn_mass).or. +! & (iq.eq.iccn_number)))) then c Computing size distribution: c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - if ((iq.eq.idust_mass).or.(iq.eq.idust_number)) then +c if ((iq.eq.idust_mass).or.(iq.eq.idust_number)) then do l=1,nlay do ig=1, ngrid @@ -284,12 +294,12 @@ end do sigma0 = varian - else - do l=1,nlay - do ig=1, ngrid - r0(ig,l)=r0ccn(ig,l) - end do - end do - sigma0 = sqrt(log(1.+nuice_sed)) - endif +c else ! ccn +c do l=1,nlay +c do ig=1, ngrid +c r0(ig,l)=r0ccn(ig,l) +c end do +c end do +c sigma0 = sqrt(log(1.+nuice_sed)) +c endif c Computing mass mixing ratio for each particle size @@ -297,5 +307,5 @@ IF ((iq.EQ.idust_mass).or.(iq.EQ.iccn_mass)) then radpower = 2 - ELSE + ELSE ! number radpower = -1 ENDIF @@ -348,13 +358,15 @@ do ir=1,nr - IF ((iq.EQ.idust_mass).or.(iq.EQ.idust_number)) then +! IF ((iq.EQ.idust_mass).or.(iq.EQ.idust_number)) then + ! Dust sedimentation call newsedim(ngrid,nlay,1,1,ptimestep, & pplev,masse,epaisseur,pt,rd(ir),rho_dust,qr(1,1,ir), & wq,0.5) - ELSE - call newsedim(ngrid,nlay,1,ngrid*nlay,ptimestep, - & pplev,masse,epaisseur,pt,rd(ir),rhocloud,qr(1,1,ir), - & wq,1.0) - ENDIF +! ELSE +! ! CCN sedimentation +! call newsedim(ngrid,nlay,1,ngrid*nlay,ptimestep, +! & pplev,masse,epaisseur,pt,rd(ir),rhocloud,qr(1,1,ir), +! & wq,beta) +! ENDIF c Tendencies @@ -371,25 +383,29 @@ enddo ! of do ir=1,nr c ----------------------------------------------------------------- -c WATER CYCLE CASE -c ----------------------------------------------------------------- - else if (water.and.(iq.eq.igcm_h2o_ice)) then - if (microphys) then - call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, - & ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rhocloud, - & zqi(1,1,iq),wq,1.0) - else - call newsedim(ngrid,nlay,ngrid*nlay,1, - & ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rho_q(iq), - & zqi(1,1,iq),wq,1.0) - endif ! of if (microphys) +c WATER CYCLE CASE --- NOW DONE IN WATERCLOUD !! +c ----------------------------------------------------------------- +! else if (water.and.(iq.eq.igcm_h2o_ice)) then +! if (microphys) then +! ! water ice sedimentation +! call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay, +! & ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rhocloud, +! & zqi(1,1,iq),wq,beta) +! else +! ! water ice sedimentation +! call newsedim(ngrid,nlay,ngrid*nlay,1, +! & ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rho_q(iq), +! & zqi(1,1,iq),wq,beta) +! endif ! of if (microphys) c Tendencies c ~~~~~~~~~~ - do ig=1,ngrid - pdqs_sed(ig,iq)=wq(ig,1)/ptimestep - end do +! do ig=1,ngrid +! pdqs_sed(ig,iq)=wq(ig,1)/ptimestep +! end do c ----------------------------------------------------------------- c GENERAL CASE c ----------------------------------------------------------------- - else +c else + else if ((iq .ne. iccn_mass) .and. (iq .ne. iccn_number) + & .and. (iq .ne. igcm_h2o_ice)) then ! because it is done in watercloud call newsedim(ngrid,nlay,1,1,ptimestep, & pplev,masse,epaisseur,pt,radius(iq),rho_q(iq), @@ -418,6 +434,5 @@ c Update the dust particle size "rdust" c ------------------------------------- - if (doubleq) then - DO l = 1, nlay + DO l = 1, nlay DO ig=1,ngrid rdust(ig,l)= @@ -426,42 +441,39 @@ rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6) ENDDO - ENDDO - endif !of if (doubleq) + ENDDO c Update the ice particle size "rice" c ------------------------------------- - IF(scavenging) THEN - DO l = 1, nlay - DO ig=1,ngrid - Mo = zqi(ig,l,igcm_h2o_ice) + - & zqi(ig,l,iccn_mass)* tauscaling(ig) + 1.e-30 - No = zqi(ig,l,iccn_number)* tauscaling(ig)+ 1.e-30 - rhocloud(ig,l) = zqi(ig,l,igcm_h2o_ice) / Mo * rho_ice - & +zqi(ig,l,iccn_mass)* tauscaling(ig) - & / Mo * rho_dust - rhocloud(ig,l) = - & min(max(rhocloud(ig,l),rho_ice),rho_dust) - rice(ig,l) = - & ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.) - if ((Mo.lt.1.e-20) .or. (No.le.1)) rice(ig,l) = 1.e-8 +! IF(scavenging) THEN +! DO l = 1, nlay +! DO ig=1,ngrid +! Mo = zqi(ig,l,igcm_h2o_ice) + +! & zqi(ig,l,iccn_mass)* tauscaling(ig) + 1.e-30 +! No = zqi(ig,l,iccn_number)* tauscaling(ig)+ 1.e-30 +! rhocloud(ig,l) = zqi(ig,l,igcm_h2o_ice) / Mo * rho_ice +! & +zqi(ig,l,iccn_mass)* tauscaling(ig) +! & / Mo * rho_dust +! rhocloud(ig,l) = +! & min(max(rhocloud(ig,l),rho_ice),rho_dust) +! rice(ig,l) = +! & ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.) +! if ((Mo.lt.1.e-15) .or. (No.le.50)) rice(ig,l) = 1.e-8 c print*, "Mice,Mo, No",zqi(ig,l,igcm_h2o_ice),Mo, No c print*, "rice, rho apres", rice(ig,l), rhocloud(ig,l) - - ENDDO - ENDDO - ELSE - if (water) then - DO l = 1, nlay - DO ig=1,ngrid - ccntyp = - & 1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-zlay(ig,l)/10000.) - ccntyp = ccntyp /ccn_factor - rice(ig,l)=max( CBRT ( (zqi(ig,l,igcm_h2o_ice)/rho_ice - & +ccntyp*(4./3.)*pi*rdust(ig,l)**3.) - & /(ccntyp*4./3.*pi) ), rdust(ig,l)) - ENDDO - ENDDO - endif ! of if (water) - ENDIF ! of IF (scavenging) +! +! ENDDO +! ENDDO +! ELSE +! DO l = 1, nlay +! DO ig=1,ngrid +! ccntyp = +! & 1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-zlay(ig,l)/10000.) +! ccntyp = ccntyp /ccn_factor +! rice(ig,l)=max( CBRT ( (zqi(ig,l,igcm_h2o_ice)/rho_ice +! & +ccntyp*(4./3.)*pi*rdust(ig,l)**3.) +! & /(ccntyp*4./3.*pi) ), rdust(ig,l)) +! ENDDO +! ENDDO +! ENDIF RETURN Index: trunk/LMDZ.MARS/libf/phymars/improvedclouds.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/improvedclouds.F (revision 517) +++ trunk/LMDZ.MARS/libf/phymars/improvedclouds.F (revision 520) @@ -1,7 +1,9 @@ subroutine improvedclouds(ngrid,nlay,ptimestep, - & pplev,pplay,pt,pdt, - & pq,pdq,pdqcloud,pdqscloud,pdtcloud, + & pplev,pplay,zlev,pt,pdt, + & pq,pdq,pdqcloud,pdtcloud, & nq,tauscaling,rdust,rice,nuice, & rsedcloud,rhocloud) +! to use 'getin' + USE ioipsl_getincom implicit none c------------------------------------------------------------------ @@ -22,4 +24,5 @@ c Authors: J.-B. Madeleine, based on the work by Franck Montmessin c (October 2011) +c T. Navarro, debug & correction (October-December 2011) c------------------------------------------------------------------ #include "dimensions.h" @@ -37,5 +40,8 @@ integer nq ! nombre de traceurs REAL ptimestep ! pas de temps physique (s) - REAL pplev(ngrid,nlay+1),pplay(ngrid,nlay) ! pression au milieu des couches (Pa) + REAL pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa) + REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa) + REAL zlev(ngrid,nlay+1) ! altitude at layer boundaries + REAL pt(ngrid,nlay) ! temperature at the middle of the ! layers (K) @@ -57,5 +63,4 @@ REAL pdqcloud(ngrid,nlay,nq) ! tendance de la condensation ! H2O(kg/kg.s-1) - REAL pdqscloud(ngrid,nq) ! flux en surface (kg.m-2.s-1) REAL pdtcloud(ngrid,nlay) ! tendance temperature due ! a la chaleur latente @@ -75,4 +80,5 @@ INTEGER ig,l,i + REAL zq(ngridmx,nlayermx,nqmx) ! local value of tracers @@ -106,10 +112,10 @@ DOUBLE PRECISION rb_cld(nbin_cld+1)! boundary values of each rad_cld bin (m) SAVE rb_cld - DOUBLE PRECISION dr_cld(nbin_cld)! width of each rad_cld bin (m) - DOUBLE PRECISION vol_cld(nbin_cld) ! particle volume for each bin (m3) + DOUBLE PRECISION dr_cld(nbin_cld) ! width of each rad_cld bin (m) + DOUBLE PRECISION vol_cld(nbin_cld) ! particle volume for each bin (m3) REAL sigma_ice ! Variance of the ice and CCN distributions SAVE sigma_ice - + c---------------------------------- c some outputs for 1D -- TEMPORARY @@ -126,8 +132,13 @@ REAL Mccn_col(ngridmx) ! total column ccn mass REAL Nccn_col(ngridmx) ! total column ccn mass + REAL dMice_col(ngridmx) ! total column ice mass change + REAL drice_col(ngridmx) ! total column ice radius change + REAL icetot(ngridmx) ! total column ice mass + REAL rice_out(ngridmx,nlayermx) ! ice radius change REAL rate_out(ngridmx,nlayermx) ! nucleation rate INTEGER count LOGICAL output_sca ! scavenging outputs flag for tests + output_sca = .false. c---------------------------------- @@ -202,5 +213,4 @@ c----------------------------------------------------------------------- c Initialize the tendencies - pdqscloud(1:ngrid,1:nq)=0 pdqcloud(1:ngrid,1:nlay,1:nq)=0 pdtcloud(1:ngrid,1:nlay)=0 @@ -249,5 +259,4 @@ enddo - c------------------------------------------------------------------ c Cloud microphysical scheme @@ -269,6 +278,6 @@ IF ((satu .ge. 1)! ) THEN ! if we have condensation - & .or. ( zq(ig,l,igcm_ccn_number) - & .ge. 10) ) THEN ! or sublimation + & .or. ( zq(ig,l,igcm_ccn_number)*tauscaling(ig) + & .ge. 2) ) THEN ! or sublimation @@ -287,5 +296,5 @@ & -derf( dlog(rb_cld(i) * Rm) * dev2 ) ) enddo - + ! sumcheck = 0 ! do i = 1, nbin_cld @@ -320,16 +329,11 @@ rate_out(ig,l) = 0 do i = 1, nbin_cld + ! schema simple n_aer(i) = n_aer(i) / ( 1. + rate(i)*ptimestep ) m_aer(i) = m_aer(i) / ( 1. + rate(i)*ptimestep ) dN = dN + n_aer(i) * rate(i) * ptimestep dM = dM + m_aer(i) * rate(i) * ptimestep - rate_out(ig,l)=rate_out(ig,l)+rate(i) + !rate_out(ig,l)=rate_out(ig,l)+rate(i) enddo - -! dN = min( max(dN,-zq(ig,l,igcm_ccn_number) ), -! & zq(ig,l,igcm_dust_number) ) -! -! dM = min( max(dM,-zq(ig,l,igcm_ccn_mass) ), -! & zq(ig,l,igcm_dust_mass) ) Mo = zq0(ig,l,igcm_h2o_ice) + @@ -382,5 +386,5 @@ zq(ig,l,igcm_h2o_vap) = zq0(ig,l,igcm_h2o_vap) - & dMice - + c If all the ice particles sublimate, all the condensation c nuclei are released: @@ -411,5 +415,6 @@ zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + dM zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + dN - + + pdqcloud(ig,l,igcm_dust_mass)=(zq(ig,l,igcm_dust_mass) & -zq0(ig,l,igcm_dust_mass))/ptimestep @@ -425,4 +430,5 @@ & -zq0(ig,l,igcm_h2o_ice))/ptimestep + count = count +1 ELSE ! for coherence (rdust, rice computations etc ..) @@ -447,6 +453,7 @@ lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)**2)*1.e+3 pdtcloud(ig,l)=-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp +c pdtcloud(ig,l)=pdt(ig,l)-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp -c----- update rice & rhocloud AFTER scavenging +c----- update rice & rhocloud AFTER microphysic Mo = zq(ig,l,igcm_h2o_ice) + & zq(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30 @@ -456,7 +463,10 @@ & / Mo * rho_dust rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust) + + rice_out(ig,l)=rice(ig,l) rice(ig,l) = & ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.) if ((Mo.lt.1.e-20) .or. (No.le.1)) rice(ig,l) = 1.e-8 + rice_out(ig,l)=rice(ig,l)-rice_out(ig,l) nuice(ig,l)=nuice_ref ! used for rad. transfer calculations @@ -474,5 +484,7 @@ ENDDO ENDDO - + +! print*, 'SATU' +! print*, satu_out(1,:) c------------------------------------------------------------------ @@ -494,4 +506,7 @@ Mccn_col(:) = 0 Nccn_col(:) = 0 + dMice_col(:) = 0 + drice_col(:) = 0 + icetot(:) = 0 do l=1, nlay do ig=1,ngrid @@ -512,6 +527,17 @@ & zq(ig,l,igcm_ccn_number)*tauscaling(ig) & *(pplev(ig,l) - pplev(ig,l+1)) / g + dMice_col(ig) = dMice_col(ig) + + & Mcon_out(ig,l) + & *(pplev(ig,l) - pplev(ig,l+1)) / g + drice_col(ig) = drice_col(ig) + + & rice_out(ig,l)*zq(ig,l,igcm_h2o_ice) + & *(pplev(ig,l) - pplev(ig,l+1)) / g + icetot(ig) = icetot(ig) + + & zq(ig,l,igcm_h2o_ice) + & *(pplev(ig,l) - pplev(ig,l+1)) / g enddo ! of do ig=1,ngrid enddo ! of do l=1,nlay + + drice_col(ig) = drice_col(ig)/icetot(ig) @@ -549,4 +575,8 @@ call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1, & dN_out) + call WRITEDIAGFI(ngrid,"dicetot","ice col var","kg/m2",0, + & dMice_col) + call WRITEDIAGFI(ngrid,"dricetot","ice col var","m",0, + & drice_col) call WRITEDIAGFI(ngrid,"mcond","h2o condensed mass","kg",1, & Mcon_out) Index: trunk/LMDZ.MARS/libf/phymars/microphys.h =================================================================== --- trunk/LMDZ.MARS/libf/phymars/microphys.h (revision 517) +++ trunk/LMDZ.MARS/libf/phymars/microphys.h (revision 520) @@ -5,5 +5,5 @@ ! Number of bins - INTEGER, PARAMETER :: nbin_cld = 50 + INTEGER, PARAMETER :: nbin_cld = 5 ! Reference temperature, T=273.15 K Index: trunk/LMDZ.MARS/libf/phymars/nuclea.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/nuclea.F (revision 517) +++ trunk/LMDZ.MARS/libf/phymars/nuclea.F (revision 520) @@ -44,8 +44,41 @@ integer i + + LOGICAL firstcall + DATA firstcall/.true./ + SAVE firstcall c ************************************************* mtetalocal = mteta + +cccccccccccccccccccccccccccccccccccccccccccccccccc +ccccccccccc ESSAIS TN MTETA = F (T) cccccccccccccc +c if (temp .gt. 200) then +c mtetalocal = mtetalocal +c else if (temp .lt. 190) then +c mtetalocal = mtetalocal-0.05 +c else +c mtetalocal = mtetalocal - (190-temp)*0.005 +c endif +c----------------exp law, see Trainer 2008, J. Phys. Chem. C 2009, 113, 2036\u20132040 + !mtetalocal = max(mtetalocal - 6005*exp(-0.065*temp),0.1) + !mtetalocal = max(mtetalocal - 6005*exp(-0.068*temp),0.1) + !print*, mtetalocal, temp +cccccccccccccccccccccccccccccccccccccccccccccccccc +cccccccccccccccccccccccccccccccccccccccccccccccccc + IF (firstcall) THEN + print*, ' ' + print*, 'dear user, please keep in mind that' + print*, 'contact parameter IS constant' + !print*, 'contact parameter IS NOT constant:' + !print*, 'max(mteta - 6005*exp(-0.065*temp),0.1)' + !print*, 'max(mteta - 6005*exp(-0.068*temp),0.1)' + print*, ' ' + firstcall=.false. + END IF +cccccccccccccccccccccccccccccccccccccccccccccccccc +cccccccccccccccccccccccccccccccccccccccccccccccccc + if (sat .gt. 1.) then ! minimum condition to activate nucleation Index: trunk/LMDZ.MARS/libf/phymars/physiq.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/physiq.F (revision 517) +++ trunk/LMDZ.MARS/libf/phymars/physiq.F (revision 520) @@ -201,6 +201,6 @@ real rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius real rhocloud(ngridmx,nlayermx) ! Cloud density (kg.m-3) - REAL surfdust(ngridmx,nlayermx) ! dust surface area (m2/m3, if photochemistry) - REAL surfice(ngridmx,nlayermx) ! ice surface area (m2/m3, if photochemistry) + REAL surfdust(ngridmx,nlayermx) ! dust surface area (micron2/cm3, if photochemistry) + REAL surfice(ngridmx,nlayermx) ! ice surface area (micron2/cm3, if photochemistry) c Variables used by the slope model @@ -308,14 +308,22 @@ real rho(ngridmx,nlayermx) ! density real vmr(ngridmx,nlayermx) ! volume mixing ratio - REAL colden(ngridmx,nqmx) ! vertical column of tracers + !real colden(ngridmx,nqmx) ! vertical column !FL REAL mtot(ngridmx) ! Total mass of water vapor (kg/m2) REAL icetot(ngridmx) ! Total mass of water ice (kg/m2) + REAL ccntot(ngridmx) ! Total number of ccn (nbr/m2) REAL rave(ngridmx) ! Mean water ice effective radius (m) REAL opTES(ngridmx,nlayermx)! abs optical depth at 825 cm-1 REAL tauTES(ngridmx) ! column optical depth at 825 cm-1 REAL Qabsice ! Water ice absorption coefficient + REAL taucloudtes(ngridmx)! Cloud opacity at infrared + ! reference wavelength using + ! Qabs instead of Qext + ! (direct comparison with TES) + c Test 1d/3d scavenging - real h2o_tot + real h2otot(ngridmx) + REAL satu(ngridmx,nlayermx) ! satu ratio for output + REAL zqsat(ngridmx,nlayermx) ! saturation REAL time_phys @@ -332,6 +340,6 @@ REAL lmax_th_out(ngridmx),zmax_th(ngridmx) - REAL, SAVE :: wstar(ngridmx) - REAL, SAVE :: hfmax_th(ngridmx) + REAL, SAVE :: wmax_th(ngridmx) + REAL hfmax_th(ngridmx) REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx) REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nqmx) @@ -342,5 +350,5 @@ REAL z_out ! height of interpolation between z0 and z1 [meters] REAL ustar(ngridmx),tstar(ngridmx) ! friction velocity and friction potential temp - REAL L_mo(ngridmx),wstarpbl(ngridmx),vhf(ngridmx),vvv(ngridmx) + REAL L_mo(ngridmx) REAL zu2(ngridmx) c======================================================================= @@ -359,5 +367,5 @@ fluxrad(:)=0 - wstar(:)=0. + wmax_th(:)=0. c read startfi @@ -569,6 +577,4 @@ enddo endif - else - zdtnlte(:,:)=0. endif @@ -589,5 +595,6 @@ $ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, $ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw, - & tauref,tau,aerosol,tauscaling,rdust,rice,nuice,co2ice) + $ tauref,tau,aerosol,tauscaling,taucloudtes,rdust,rice, + $ nuice,co2ice) c Outputs for basic check (middle of domain) @@ -747,10 +754,6 @@ DO ig=1, ngridmx IF (zh(ig,1) .lt. tsurf(ig)) THEN - wstar(ig)=1. - hfmax_th(ig)=0.2 - ELSE - wstar(ig)=0. - hfmax_th(ig)=0. - ENDIF + wmax_th(ig)=1. + ENDIF ENDDO ENDIF @@ -768,5 +771,5 @@ $ zdum1,zdum2,zdh,pdq,zflubid, $ zdudif,zdvdif,zdhdif,zdtsdif,q2, - & zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th) + & zdqdif,zdqsdif,wmax_th,zcdv,zcdh) #ifdef MESOSCALE @@ -829,5 +832,5 @@ $ pplay,pplev,pphi,zpopsk, $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th, - $ dtke_th,hfmax_th,wstar) + $ dtke_th,hfmax_th,wmax_th) DO l=1,nlayer @@ -858,6 +861,5 @@ else !of if calltherm lmax_th(:)=0 - wstar(:)=0. - hfmax_th(:)=0. + wmax_th(:)=0. lmax_th_out(:)=0. end if @@ -911,10 +913,10 @@ c ------------------------------------------- - IF (tituscap) THEN - !!! get the actual co2 seasonal cap from Titus observations - CALL geticecover( ngrid, 180.*zls/pi, +#ifdef MESOSCALE + !!! get the actual co2 seasonal cap from Titus observations + CALL geticecover( ngrid, 180.*zls/pi, . 180.*long/pi, 180.*lati/pi, co2ice ) - co2ice = co2ice * 10000. - ENDIF + co2ice = co2ice * 10000. +#endif IF (callcond) THEN @@ -962,4 +964,5 @@ c ---------------------------------------- IF (water) THEN + call watercloud(ngrid,nlayer,ptimestep, @@ -976,5 +979,5 @@ ENDDO endif - + ! increment water vapour and ice atmospheric tracers tendencies IF (water) THEN @@ -987,4 +990,10 @@ & pdq(ig,l,igcm_h2o_ice)+ & zdqcloud(ig,l,igcm_h2o_ice) + if (((pq(ig,l,igcm_h2o_ice) + + & pdq(ig,l,igcm_h2o_ice)*ptimestep)) .le. 0) + & then + pdq(ig,l,igcm_h2o_ice) = + & - pq(ig,l,igcm_h2o_ice)/ptimestep + 1e-30 + endif IF (scavenging) THEN pdq(ig,l,igcm_ccn_mass)= @@ -994,4 +1003,10 @@ & pdq(ig,l,igcm_ccn_number)+ & zdqcloud(ig,l,igcm_ccn_number) + pdq(ig,l,igcm_dust_mass)= + & pdq(ig,l,igcm_dust_mass)+ + & zdqcloud(ig,l,igcm_dust_mass) + pdq(ig,l,igcm_dust_number)= + & pdq(ig,l,igcm_dust_number)+ + & zdqcloud(ig,l,igcm_dust_number) !!!!!!!!!!!!!!!!!!!!! We need to check that we have Nccn & Ndust > 0 !!!!!!!!!!!!!!!!!!!!! This is due to single precision rounding problems @@ -1001,4 +1016,6 @@ pdq(ig,l,igcm_ccn_number) = & - pq(ig,l,igcm_ccn_number)/ptimestep + 1e-30 + pdq(ig,l,igcm_ccn_mass) = + & - pq(ig,l,igcm_ccn_mass)/ptimestep + 1e-30 endif if (((pq(ig,l,igcm_dust_number) + @@ -1007,13 +1024,9 @@ pdq(ig,l,igcm_dust_number) = & - pq(ig,l,igcm_dust_number)/ptimestep + 1e-30 + pdq(ig,l,igcm_dust_mass) = + & - pq(ig,l,igcm_dust_mass)/ptimestep + 1e-30 endif !!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!! - pdq(ig,l,igcm_dust_mass)= - & pdq(ig,l,igcm_dust_mass)+ - & zdqcloud(ig,l,igcm_dust_mass) - pdq(ig,l,igcm_dust_number)= - & pdq(ig,l,igcm_dust_number)+ - & zdqcloud(ig,l,igcm_dust_number) ENDIF ENDDO @@ -1029,5 +1042,58 @@ END IF ! of IF (water) -c 7b. Aerosol particles + +c 7b. Chemical species +c ------------------ + +#ifndef MESOSCALE +c -------------- +c photochemistry : +c -------------- + IF (photochem .or. thermochem) then +!NB: Photochemistry includes condensation of H2O2 + PRINT*, 'SURFDUST,SURFICE TO BE IMPLEMENTED. YAAAAAARG.' + call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0, + $ zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim, + $ zdqcloud,zdqscloud,tauref,co2ice, + $ pu,pdu,pv,pdv,surfdust,surfice) + + ! increment values of tracers: + DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry + ! tracers is zero anyways + DO l=1,nlayer + DO ig=1,ngrid + pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq) + ENDDO + ENDDO + ENDDO ! of DO iq=1,nq + ! add condensation tendency for H2O2 + if (igcm_h2o2.ne.0) then + DO l=1,nlayer + DO ig=1,ngrid + pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2) + & +zdqcloud(ig,l,igcm_h2o2) + ENDDO + ENDDO + endif + + ! increment surface values of tracers: + DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry + ! tracers is zero anyways + DO ig=1,ngrid + dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq) + ENDDO + ENDDO ! of DO iq=1,nq + ! add condensation tendency for H2O2 + if (igcm_h2o2.ne.0) then + DO ig=1,ngrid + dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2) + & +zdqscloud(ig,igcm_h2o2) + ENDDO + endif + + END IF ! of IF (photochem.or.thermochem) +#endif + +c 7c. Aerosol particles c ------------------- @@ -1070,4 +1136,7 @@ & pq, pdq, zdqsed, zdqssed,nq, & tau,tauscaling) + + + !print*, 'h2o_ice zdqsed ds physiq', zdqsed(1,:,igcm_h2o_ice) DO iq=1, nq @@ -1084,63 +1153,5 @@ ENDDO END IF ! of IF (sedimentation) - -c 7c. Chemical species -c ------------------ - -#ifndef MESOSCALE -c -------------- -c photochemistry : -c -------------- - IF (photochem .or. thermochem) then - -! dust and ice surface area - call surfacearea(ngrid, nlayer, ptimestep, pplay, zzlay, - $ pt, pq, pdq, nq, - $ rdust, rice, tau, tauscaling, - $ surfdust, surfice) -! call photochemistry - call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0, - $ zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim, - $ zdqcloud,zdqscloud,tauref,co2ice, - $ pu,pdu,pv,pdv,surfdust,surfice) - - ! increment values of tracers: - DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry - ! tracers is zero anyways - DO l=1,nlayer - DO ig=1,ngrid - pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq) - ENDDO - ENDDO - ENDDO ! of DO iq=1,nq - - ! add condensation tendency for H2O2 - if (igcm_h2o2.ne.0) then - DO l=1,nlayer - DO ig=1,ngrid - pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2) - & +zdqcloud(ig,l,igcm_h2o2) - ENDDO - ENDDO - endif - - ! increment surface values of tracers: - DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry - ! tracers is zero anyways - DO ig=1,ngrid - dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq) - ENDDO - ENDDO ! of DO iq=1,nq - - ! add condensation tendency for H2O2 - if (igcm_h2o2.ne.0) then - DO ig=1,ngrid - dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2) - & +zdqscloud(ig,igcm_h2o2) - ENDDO - endif - - END IF ! of IF (photochem.or.thermochem) -#endif + @@ -1158,4 +1169,5 @@ ENDDO ! (ig) ENDDO ! (iq) + endif ! of if (tracer) @@ -1232,5 +1244,5 @@ ENDIF ! of IF (tracer.AND.water.AND.(ngridmx.NE.1)) -c + c 9.2 Compute soil temperatures and subsurface heat flux: c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1291,5 +1303,5 @@ DO ig=1,ngridmx DO l=1,nlayermx - zh(ig,l) = zt(ig,l)*(zplev(ig,1)/zplay(ig,l))**rcp + zh(ig,l) = zt(ig,l)*(zplay(ig,l)/zplev(ig,1))**rcp ENDDO ENDDO @@ -1374,4 +1386,15 @@ if (tracer) then if (water) then + + if (scavenging) then + ccntot(:)= 0 + do ig=1,ngrid + do l=1,nlayermx + ccntot(ig) = ccntot(ig) + + & zq(ig,l,igcm_ccn_number)*tauscaling(ig) + & *(pplev(ig,l) - pplev(ig,l+1)) / g + enddo + enddo + endif mtot(:)=0 @@ -1420,26 +1443,26 @@ call wstats(ngrid,"ps","Surface pressure","Pa",2,ps) call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf) - call wstats(ngrid,"co2ice","CO2 ice cover", - & "kg.m-2",2,co2ice) - call wstats(ngrid,"fluxsurf_lw", - & "Thermal IR radiative flux to surface","W.m-2",2, - & fluxsurf_lw) - call wstats(ngrid,"fluxsurf_sw", - & "Solar radiative flux to surface","W.m-2",2, - & fluxsurf_sw_tot) - call wstats(ngrid,"fluxtop_lw", - & "Thermal IR radiative flux to space","W.m-2",2, - & fluxtop_lw) - call wstats(ngrid,"fluxtop_sw", - & "Solar radiative flux to space","W.m-2",2, - & fluxtop_sw_tot) - call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt) - call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu) - call wstats(ngrid,"v","Meridional (North-South) wind", - & "m.s-1",3,zv) - call wstats(ngrid,"w","Vertical (down-up) wind", - & "m.s-1",3,pw) - call wstats(ngrid,"rho","Atmospheric density","none",3,rho) - call wstats(ngrid,"pressure","Pressure","Pa",3,pplay) +c call wstats(ngrid,"co2ice","CO2 ice cover", +c & "kg.m-2",2,co2ice) +c call wstats(ngrid,"fluxsurf_lw", +c & "Thermal IR radiative flux to surface","W.m-2",2, +c & fluxsurf_lw) +c call wstats(ngrid,"fluxsurf_sw", +c & "Solar radiative flux to surface","W.m-2",2, +c & fluxsurf_sw_tot) +c call wstats(ngrid,"fluxtop_lw", +c & "Thermal IR radiative flux to space","W.m-2",2, +c & fluxtop_lw) +c call wstats(ngrid,"fluxtop_sw", +c & "Solar radiative flux to space","W.m-2",2, +c & fluxtop_sw_tot) +c call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt) +c call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu) +c call wstats(ngrid,"v","Meridional (North-South) wind", +c & "m.s-1",3,zv) +c call wstats(ngrid,"w","Vertical (down-up) wind", +c & "m.s-1",3,pw) +c call wstats(ngrid,"rho","Atmospheric density","none",3,rho) +c call wstats(ngrid,"pressure","Pressure","Pa",3,pplay) c call wstats(ngrid,"q2", c & "Boundary layer eddy kinetic energy", @@ -1456,20 +1479,20 @@ if (tracer) then if (water) then -c vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap) -c & *mugaz/mmol(igcm_h2o_vap) -c call wstats(ngrid,"vmr_h2ovapor", -c & "H2O vapor volume mixing ratio","mol/mol", -c & 3,vmr) -c vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice) -c & *mugaz/mmol(igcm_h2o_ice) -c call wstats(ngrid,"vmr_h2oice", -c & "H2O ice volume mixing ratio","mol/mol", -c & 3,vmr) + vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap) + & *mugaz/mmol(igcm_h2o_vap) + call wstats(ngrid,"vmr_h2ovapor", + & "H2O vapor volume mixing ratio","mol/mol", + & 3,vmr) + vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice) + & *mugaz/mmol(igcm_h2o_ice) + call wstats(ngrid,"vmr_h2oice", + & "H2O ice volume mixing ratio","mol/mol", + & 3,vmr) call wstats(ngrid,"h2o_ice_s", & "surface h2o_ice","kg/m2", & 2,qsurf(1,igcm_h2o_ice)) -c call wstats(ngrid,'albedo', -c & 'albedo', -c & '',2,albedo(1:ngridmx,1)) + call wstats(ngrid,'albedo', + & 'albedo', + & '',2,albedo(1:ngridmx,1)) call wstats(ngrid,"mtot", & "total mass of water vapor","kg/m2", @@ -1478,15 +1501,21 @@ & "total mass of water ice","kg/m2", & 2,icetot) -c call wstats(ngrid,"reffice", -c & "Mean reff","m", -c & 2,rave) -c call wstats(ngrid,"rice", -c & "Ice particle size","m", -c & 3,rice) -c If activice is true, tauTES is computed in aeropacity.F; + call wstats(ngrid,"reffice", + & "Mean reff","m", + & 2,rave) + call wstats(ngrid,"ccntot", + & "condensation nuclei","Nbr/m2", + & 2,ccntot) + call wstats(ngrid,"rice", + & "Ice particle size","m", + & 3,rice) if (.not.activice) then call wstats(ngrid,"tauTESap", & "tau abs 825 cm-1","", & 2,tauTES) + else + call wstats(ngridmx,'tauTES', + & 'tau abs 825 cm-1', + & '',2,taucloudtes) endif @@ -1495,35 +1524,28 @@ if (thermochem.or.photochem) then do iq=1,nq - if (noms(iq) .ne. "dust_mass" .and. - $ noms(iq) .ne. "dust_number" .and. - $ noms(iq) .ne. "ccn_mass" .and. - $ noms(iq) .ne. "ccn_number") then - do l=1,nlayer - do ig=1,ngrid - vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq) - end do - end do - call wstats(ngrid,"vmr_"//trim(noms(iq)), - $ "Volume mixing ratio","mol/mol",3,vmr) - if ((noms(iq).eq."o") .or. (noms(iq).eq."co2").or. - $ (noms(iq).eq."o3")) then - call writediagfi(ngrid,"vmr_"//trim(noms(iq)), - $ "Volume mixing ratio","mol/mol",3,vmr) - end if - do ig = 1,ngrid - colden(ig,iq) = 0. - end do - do l=1,nlayer - do ig=1,ngrid - colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq) - $ *(pplev(ig,l)-pplev(ig,l+1)) - $ *6.022e22/(mmol(iq)*g) - end do - end do - call wstats(ngrid,"c_"//trim(noms(iq)), - $ "column","mol cm-2",2,colden(1,iq)) - call writediagfi(ngrid,"c_"//trim(noms(iq)), - $ "column","mol cm-2",2,colden(1,iq)) - end if + if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or. + . (noms(iq).eq."co").or.(noms(iq).eq."n2").or. + . (noms(iq).eq."h2").or. + . (noms(iq).eq."o3")) then + do l=1,nlayer + do ig=1,ngrid + vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq) + end do + end do + call wstats(ngrid,"vmr_"//trim(noms(iq)), + . "Volume mixing ratio","mol/mol",3,vmr) + endif +! do ig = 1,ngrid +! colden(ig,iq) = 0. !FL +! end do +! do l=1,nlayer !FL +! do ig=1,ngrid !FL +! colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq) !FL +! $ *(pplev(ig,l)-pplev(ig,l+1)) !FL +! $ *6.022e22/(mmol(iq)*g) !FL +! end do !FL +! end do !FL +! call wstats(ngrid,"c_"//trim(noms(iq)), !FL +! $ "column","mol cm-2",2,colden(1,iq)) !FL end do end if ! of if (thermochem.or.photochem) @@ -1583,6 +1605,6 @@ c call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2, c & emis) - call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay) -! call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev) +c call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay) +c call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev) call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2, & tsurf) @@ -1601,11 +1623,11 @@ c call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2, c & fluxtop_sw_tot) - call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt) - call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu) - call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv) - call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw) - call WRITEDIAGFI(ngrid,"rho","density","none",3,rho) +c call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt) +c call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu) +c call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv) +c call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw) +c call WRITEDIAGFI(ngrid,"rho","density","none",3,rho) c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2) -! call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh) +c call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh) c call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay) c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2, @@ -1615,7 +1637,13 @@ c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate', c & 'w.m-2',3,zdtlw) -c CALL WRITEDIAGFI(ngridmx,'tauTESap', -c & 'tau abs 825 cm-1', -c & '',2,tauTES) + if (.not.activice) then + CALL WRITEDIAGFI(ngridmx,'tauTESap', + & 'tau abs 825 cm-1', + & '',2,tauTES) + else + CALL WRITEDIAGFI(ngridmx,'tauTES', + & 'tau abs 825 cm-1', + & '',2,taucloudtes) + endif #ifdef MESOINI @@ -1650,4 +1678,13 @@ call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2, & co2col) +!!!!! FL +! do iq = 1,nq +! if (noms(iq) .ne. "dust_mass" .and. +! $ noms(iq) .ne. "dust_number") then +! call writediagfi(ngrid,"c_"//trim(noms(iq)), +! $ "column","mol cm-2",2,colden(1,iq)) +! end if +! end do +!!!!! FL endif ! of if (tracer.and.(igcm_co2.ne.0)) @@ -1671,11 +1708,4 @@ & 'kg.m-2',2, & qsurf(1:ngridmx,igcm_h2o_ice)) - if (photochem) then - do iq = 1,nq - call writediagfi( ngrid,trim(noms(iq)), - $ 'mix rat','units', - $ 3,zq(1:ngridmx,1:nlayermx,iq) ) - enddo - endif #endif @@ -1697,24 +1727,21 @@ c & 'Mean reff', c & 'm',2,rave) +c CALL WRITEDIAGFI(ngrid,"ccntot", +c & "condensation nuclei","Nbr/m2", +c & 2,ccntot) c call WRITEDIAGFI(ngridmx,'rice','Ice particle size', c & 'm',3,rice) -c If activice is true, tauTES is computed in aeropacity.F; - if (.not.activice) then - CALL WRITEDIAGFI(ngridmx,'tauTESap', - & 'tau abs 825 cm-1', - & '',2,tauTES) - endif call WRITEDIAGFI(ngridmx,'h2o_ice_s', & 'surface h2o_ice', & 'kg.m-2',2,qsurf(1,igcm_h2o_ice)) - if (caps) then - do ig=1,ngridmx - if (watercaptag(ig)) watercapflag(ig) = 1 - enddo +c if (caps) then +c do ig=1,ngridmx +c if (watercaptag(ig)) watercapflag(ig) = 1 +c enddo c CALL WRITEDIAGFI(ngridmx,'watercaptag', -c & 'Ice water caps', -c & '',2,watercapflag) - endif +c & 'Ice water caps', +c & '',2,watercapflag) +cs endif c CALL WRITEDIAGFI(ngridmx,'albedo', c & 'albedo', @@ -1761,6 +1788,4 @@ c call WRITEDIAGFI(ngridmx,'reffdust','reffdust', c & 'm',3,rdust*ref_r0) -c call WRITEDIAGFI(ngridmx,'rdust','rdust', -c & 'm',3,rdust) c call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr', c & 'kg/kg',3,pq(1,1,igcm_dust_mass)) @@ -1768,8 +1793,6 @@ c & 'part/kg',3,pq(1,1,igcm_dust_number)) #ifdef MESOINI - call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr', - & 'kg/kg',3,pq(1:ngridmx,1:nlayermx,igcm_dust_mass)) call WRITEDIAGFI(ngridmx,'dustN','Dust number', - & 'part/kg',3,pq(1:ngridmx,1:nlayermx,igcm_dust_number)) + & 'part/kg',3,pq(1,1,igcm_dust_number)) #endif else @@ -1784,8 +1807,8 @@ if (scavenging) then - call WRITEDIAGFI(ngridmx,'ccnq','CCN mass mr', - & 'kg/kg',3,pq(1,1,igcm_ccn_mass)) - call WRITEDIAGFI(ngridmx,'ccnN','CCN number', - & 'part/kg',3,pq(1,1,igcm_ccn_number)) +c call WRITEDIAGFI(ngridmx,'ccnq','CCN mass mr', +c & 'kg/kg',3,pq(1,1,igcm_ccn_mass)) +c call WRITEDIAGFI(ngridmx,'ccnN','CCN number', +c & 'part/kg',3,pq(1,1,igcm_ccn_number)) endif ! (scavenging) @@ -1810,9 +1833,27 @@ z_out=0. if (calltherm .and. (z_out .gt. 0.)) then - - call pbl_parameters(ngrid,nlayer,z0, - & g,zzlay,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out, - & Teta_out,u_out,ustar,tstar,wstarpbl,L_mo,vhf,vvv) - + call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th + & ,tsurf,zh,z_out,Teta_out,u_out,ustar,tstar,L_mo) + + zu2(:)=sqrt(zu(:,1)*zu(:,1)+zv(:,1)*zv(:,1)) + call WRITEDIAGFI(ngridmx,'sqrt(zu2)', + & 'horizontal velocity norm','m/s', + & 2,zu2) + + call WRITEDIAGFI(ngridmx,'Teta_out', + & 'potential temperature at z_out','K', + & 2,Teta_out) + call WRITEDIAGFI(ngridmx,'u_out', + & 'horizontal velocity norm at z_out','m/s', + & 2,u_out) + call WRITEDIAGFI(ngridmx,'u*', + & 'friction velocity','m/s', + & 2,ustar) + call WRITEDIAGFI(ngridmx,'teta*', + & 'friction potential temperature','K', + & 2,tstar) + call WRITEDIAGFI(ngrid,'L', + & 'Monin Obukhov length','m', + & 2,L_mo) else if((.not. calltherm).and.(z_out .gt. 0.)) then @@ -1851,7 +1892,7 @@ & 'maximum TH heat flux','K.m/s', & 2,hfmax_th) - call WRITEDIAGFI(ngridmx,'wstar', + call WRITEDIAGFI(ngridmx,'wmax_th', & 'maximum TH vertical velocity','m/s', - & 2,wstar) + & 2,wmax_th) endif @@ -1903,9 +1944,24 @@ z_out=0. if (calltherm .and. (z_out .gt. 0.)) then - - call pbl_parameters(ngrid,nlayer,z0, - & g,zzlay,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out, - & Teta_out,u_out,ustar,tstar,wstarpbl,L_mo,vhf,vvv) - + call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th + & ,tsurf,zh,z_out,Teta_out,u_out,ustar,tstar,L_mo) + + zu2(:)=sqrt(zu(:,1)*zu(:,1)+zv(:,1)*zv(:,1)) + call WRITEDIAGFI(ngridmx,'sqrt(zu2)', + & 'horizontal velocity norm','m/s', + & 0,zu2) + + call WRITEDIAGFI(ngridmx,'Teta_out', + & 'potential temperature at z_out','K', + & 0,Teta_out) + call WRITEDIAGFI(ngridmx,'u_out', + & 'horizontal velocity norm at z_out','m/s', + & 0,u_out) + call WRITEDIAGFI(ngridmx,'u*', + & 'friction velocity','m/s', + & 0,ustar) + call WRITEDIAGFI(ngridmx,'teta*', + & 'friction potential temperature','K', + & 0,tstar) else if((.not. calltherm).and.(z_out .gt. 0.)) then @@ -1921,13 +1977,10 @@ & 'hauteur du thermique','point', & 0,lmax_th_out) - call WRITEDIAGFI(ngridmx,'zmax_th', - & 'hauteur du thermique','m', - & 0,zmax_th) call WRITEDIAGFI(ngridmx,'hfmax_th', & 'maximum TH heat flux','K.m/s', & 0,hfmax_th) - call WRITEDIAGFI(ngridmx,'wstar', + call WRITEDIAGFI(ngridmx,'wmax_th', & 'maximum TH vertical velocity','m/s', - & 0,wstar) + & 0,wmax_th) co2col(:)=0. @@ -1973,17 +2026,63 @@ end if -cccccccccccccccccccccc scavenging outputs 1D TN ccccccccccccccccccc +cccccccccccccccccc scavenging & water outputs 1D TN ccccccccccccccc ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - IF (scavenging) THEN - CALL WRITEDIAGFI(ngridmx,'tauTESap', + IF (water) THEN + CALL WRITEDIAGFI(ngridmx,'tauTESap', & 'tau abs 825 cm-1', - & '',1,tauTES) + & '',0,tauTES) - h2o_tot = qsurf(1,igcm_h2o_ice) + CALL WRITEDIAGFI(ngridmx,'tauTES', + & 'tau abs 825 cm-1', + & '',0,taucloudtes) + + mtot = 0 + icetot = 0 + h2otot = qsurf(1,igcm_h2o_ice) + rave = 0 do l=1,nlayer - h2o_tot = h2o_tot + - & (zq(1,l,igcm_h2o_vap) + zq(1,l,igcm_h2o_ice)) - & * (pplev(1,l) - pplev(1,l+1)) / g + mtot = mtot + zq(1,l,igcm_h2o_vap) + & * (pplev(1,l) - pplev(1,l+1)) / g + icetot = icetot + zq(1,l,igcm_h2o_ice) + & * (pplev(1,l) - pplev(1,l+1)) / g + rave = rave + zq(1,l,igcm_h2o_ice) + & * (pplev(1,l) - pplev(1,l+1)) / g + & * rice(1,l) * (1.+nuice_ref) end do + rave=rave/max(icetot,1.e-30) + h2otot = h2otot+mtot+icetot + + + if (scavenging) then + ccntot= 0 + call watersat(ngridmx*nlayermx,zt,pplay,zqsat) + do l=1,nlayermx + ccntot = ccntot + + & zq(1,l,igcm_ccn_number)*tauscaling(1) + & *(pplev(1,l) - pplev(1,l+1)) / g + satu(1,l) = zq(1,l,igcm_h2o_vap)/zqsat(1,l) + satu(1,l) = (max(satu(1,l),float(1))-1) +! & * zq(1,l,igcm_h2o_vap) * +! & (pplev(1,l) - pplev(1,l+1)) / g + enddo + + CALL WRITEDIAGFI(ngridmx,'ccntot', + & 'ccntot', + & 'nbr/m2',0,ccntot) + endif + + + CALL WRITEDIAGFI(ngridmx,'h2otot', + & 'h2otot', + & 'kg/m2',0,h2otot) + CALL WRITEDIAGFI(ngridmx,'mtot', + & 'mtot', + & 'kg/m2',0,mtot) + CALL WRITEDIAGFI(ngridmx,'icetot', + & 'icetot', + & 'kg/m2',0,icetot) + CALL WRITEDIAGFI(ngridmx,'reffice', + & 'reffice', + & 'm',0,rave) @@ -1994,17 +2093,35 @@ - call WRITEDIAGFI(ngridmx,'dqssed q','sedimentation q', - & 'kg.m-2.s-1',0,zdqssed(:,igcm_dust_number)) - call WRITEDIAGFI(ngridmx,'dqsdif q','diffusion q', - & 'kg.m-2.s-1',0,zdqsdif(:,igcm_dust_number)) + call WRITEDIAGFI(ngridmx,'zdqsed_dustq','sedimentation q', + & 'kg.m-2.s-1',1,zdqsed(1,:,igcm_dust_mass)) + call WRITEDIAGFI(ngridmx,'zdqsed_dustN','sedimentation N', + & 'Nbr.m-2.s-1',1,zdqsed(1,:,igcm_dust_number)) - call WRITEDIAGFI(ngridmx,'dqssed N','sedimentation N', - & 'kg.m-2.s-1',0,zdqssed(:,igcm_dust_number)) - call WRITEDIAGFI(ngridmx,'dqsdif N','diffusion N', - & 'kg.m-2.s-1',0,zdqsdif(:,igcm_dust_number)) + call WRITEDIAGFI(ngridmx,'zdqcloud*dt ice','cloud ice', + & 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice)*ptimestep) + call WRITEDIAGFI(ngridmx,'zdqcloud*dt vap','cloud vap', + & 'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_vap)*ptimestep) + call WRITEDIAGFI(ngridmx,'zdqdif*dt ice','dif ice', + & 'kg.m-2.s-1',1,zdqdif(1,:,igcm_h2o_ice)*ptimestep) + call WRITEDIAGFI(ngridmx,'zdqdif*dt vap','dif vap', + & 'kg.m-2.s-1',1,zdqdif(1,:,igcm_h2o_vap)*ptimestep) + call WRITEDIAGFI(ngridmx,'zdqdif*dt vap 0','dif vap', + & 'kg.m-2.s-1',0,zdqdif(1,1,igcm_h2o_vap)*ptimestep) + +! if (scavenging) then +! call WRITEDIAGFI(ngridmx,'zdqsed_ccnq','sedimentation q', +! & 'kg.m-2.s-1',1,zdqsed(1,:,igcm_ccn_mass)) +! call WRITEDIAGFI(ngridmx,'zdqsed_ccnN','sedimentation N', +! & 'Nbr.m-2.s-1',1,zdqsed(1,:,igcm_ccn_number)) +! endif +! call WRITEDIAGFI(ngridmx,'zdqsed_ice','sedimentation q', +! & 'kg.m-2.s-1',1,zdqsed(1,:,igcm_h2o_ice)) +! - call WRITEDIAGFI(ngrid,"rice","ice radius","m",1, + call WRITEDIAGFI(ngrid,"rice","ice radius","m",1, & rice) - ENDIF + call WRITEDIAGFI(ngrid,"satu","vap in satu","kg/kg",1, + & satu) + ENDIF ! of IF (water) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc Index: trunk/LMDZ.MARS/libf/phymars/simpleclouds.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/simpleclouds.F (revision 517) +++ trunk/LMDZ.MARS/libf/phymars/simpleclouds.F (revision 520) @@ -1,5 +1,5 @@ subroutine simpleclouds(ngrid,nlay,ptimestep, & pplev,pplay,pzlev,pzlay,pt,pdt, - & pq,pdq,pdqcloud,pdqscloud,pdtcloud, + & pq,pdq,pdqcloud,pdtcloud, & nq,tau,rice,nuice,rsedcloud) implicit none @@ -62,5 +62,4 @@ real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation ! H2O(kg/kg.s-1) - real pdqscloud(ngrid,nq) ! flux en surface (kg.m-2.s-1) REAL pdtcloud(ngrid,nlay) ! tendance temperature due ! a la chaleur latente @@ -141,7 +140,6 @@ enddo ! of do ig=1,ngrid - enddo ! of dol=1,nlay - - pdqscloud(1:ngrid,1:nq)=0 + enddo ! of do l=1,nlay + pdqcloud(1:ngrid,1:nlay,1:nq)=0 pdtcloud(1:ngrid,1:nlay)=0 @@ -211,12 +209,12 @@ c------------------------------------------------------------------ c TEST_JBM - IF (ngrid.eq.1) THEN - call WRITEDIAGFI(ngrid,"mcond","h2o condensed mass","kg",1, - & Mcon_out) - call WRITEDIAGFI(ngrid,"rdusttyp","rdusttyp","m",1, - & rdusttyp) - call WRITEDIAGFI(ngrid,"ccntyp","ccntyp","kg-1",1, - & ccntyp) - ENDIF +! IF (ngrid.eq.1) THEN +! call WRITEDIAGFI(ngrid,"mcond","h2o condensed mass","kg",1, +! & Mcon_out) +! call WRITEDIAGFI(ngrid,"rdusttyp","rdusttyp","m",1, +! & rdusttyp) +! call WRITEDIAGFI(ngrid,"ccntyp","ccntyp","kg-1",1, +! & ccntyp) +! ENDIF c------------------------------------------------------------------ return Index: trunk/LMDZ.MARS/libf/phymars/surfini.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/surfini.F (revision 517) +++ trunk/LMDZ.MARS/libf/phymars/surfini.F (revision 520) @@ -1,3 +1,5 @@ SUBROUTINE surfini(ngrid,piceco2,qsurf,psolaralb) + ! to use 'getin' + USE ioipsl_getincom IMPLICIT NONE c======================================================================= @@ -20,7 +22,14 @@ REAL piceco2(ngrid),psolaralb(ngrid,2) REAL qsurf(ngrid,nqmx) !tracer on surface (kg/m2) + REAL icedryness ! ice dryness EXTERNAL ISMIN,ISMAX INTEGER ISMIN,ISMAX + +! Choose false to have a somewhat non resolution dependant water ice caps distribution, +! i.e based only on lat & lon values of each physical point. +! Choose true to get a carefully choosen distribution for GCM resolutions 32x24 or 64x48 +! For vizualisation : soon ... + LOGICAL,SAVE :: improvedicelocation = .true. c c======================================================================= @@ -46,8 +55,16 @@ alternate = 0 - - do ig=1,ngridmx + temptag = .true. + + write(*,*) "Ice dryness ?" + icedryness=1. ! default value + call getin("icedryness",icedryness) + write(*,*) " icedryness = ",icedryness + #ifdef MESOSCALE + + do ig=1,ngridmx + !write(*,*) "all qsurf to zero. dirty." do iq=1,nqmx @@ -66,9 +83,82 @@ watercaptag(ig) = .false. dryness(ig) = 1. - endif + endif + + endif #else - - dryness (ig) = 1 - + ! print*,'ngridmx',ngridmx + + if (improvedicelocation) then + + if (ngridmx .eq. 738) then ! hopefully 32x24 + print*,'water ice caps distribution for 32x24 resolution:' + watercaptag(1:9) = .true. ! central cap - core + watercaptag(26:33) = .true. ! central cap + + else if (ngridmx .eq. 3010) then ! hopefully 64x48 + print*,'water ice caps distribution for 64x48 resolution:' + watercaptag(1:65) = .true. ! central cap - core +! watercaptag(85) = .true. ! central cap + watercaptag(83:85) = .true. ! central cap - BIGGER + watercaptag(93:114) = .true. ! central cap +! watercaptag(254) = .true. ! Korolev crater +! watercaptag(255) = .true. ! Korolev crater --DECALE + watercaptag(242:257)= .true. ! Korolev crater -- VERY BIG +!-------------------------------------------------------------------- +! watercaptag(136) = .true. ! outlier, lat = 78.75 +! watercaptag(138) = .true. ! outlier, lat = 78.75 +! watercaptag(140) = .true. ! outlier, lat = 78.75 +! watercaptag(142) = .true. ! outlier, lat = 78.75 +! watercaptag(183) = .true. ! outlier, lat = 78.75 +! watercaptag(185) = .true. ! outlier, lat = 78.75 +! watercaptag(187) = .true. ! outlier, lat = 78.75 +! watercaptag(189) = .true. ! outlier, lat = 78.75 +! watercaptag(191) = .true. ! outlier, lat = 78.75 +! watercaptag(193) = .true. ! outlier, lat = 78.75 +!-------------------------------------------------------------------- + watercaptag(135:142)= .true. ! BIG outlier, lat = 78.75 + watercaptag(181:193)= .true. ! BIG outlier, lat = 78.75 + + else if (ngridmx .ne. 1) then + print*,'No improved ice location for this resolution!' + print*,'Please set improvedicelocation to false in surfini.' + print*,'And check lat and lon conditions for ice caps in code.' + call exit(1) + + endif + + ! print caps locations + print*,'latitude | longitude | ig' + do ig=1,ngridmx + dryness (ig) = icedryness + !!! OU alors tu mets dryness = icedrag sur outliers et 1 ailleurs + !!! OU remettre comme c'était avant Aymeric et pis c'est tout + if (watercaptag(ig)) then + print*,'ice water cap', lati(ig)*180./pi, + . long(ig)*180./pi, ig + endif + enddo + +!!------------------------ test ----------------------------- +!!----------------------------------------------------------- +!! else if (1) then ! DUMMY !!, CAREFUL !!, TOUT CA .... +!! +!! do ig=1,ngridmx +!! dryness (ig) = 1 +!! if (albedodat(ig) .ge. 0.35) then +!! watercaptag(ig) = .true. +!! print*,' albedo criteria ice water cap', lati(ig)*180./pi, +!! . long(ig)*180./pi, ig +!! endif +!! enddo +!!----------------------------------------------------------- +!!----------------------------------------------------------- + + else + + do ig=1,ngridmx + + dryness (ig) = icedryness + !!c Towards olympia planitia water caps ('relatively' low latitude ones) !!c---------------- proposition par AS -------------------- @@ -98,9 +188,9 @@ !! . ( long(ig)*180./pi .ge. 155. ) .and. !! . ( long(ig)*180./pi .le. 180. ) ) -!!c . .or. -!!c . ( ( lati(ig)*180./pi .ge. 71. ) .and. ! cap #4 (Korolev crater) -!!c . ( lati(ig)*180./pi .le. 72. ) .and. -!!c . ( long(ig)*180./pi .ge. 163. ) .and. -!!c . ( long(ig)*180./pi .le. 168. ) ) +!! . .or. +!! . ( ( lati(ig)*180./pi .ge. 71. ) .and. ! cap #4 (Korolev crater) +!! . ( lati(ig)*180./pi .le. 72. ) .and. +!! . ( long(ig)*180./pi .ge. 163. ) .and. +!! . ( long(ig)*180./pi .le. 168. ) ) !! . .or. !! . ( ( lati(ig)*180./pi .ge. 74.9 ) .and. ! cap #5 @@ -109,48 +199,69 @@ !! . ( long(ig)*180./pi .le. -120.) ) ) !! . then -!! -!! if (temptag) then -!! -!! if ((alternate .eq. 0)) then !!! 1/2 en 64x48 sinon trop large en lat -!! if (ngridmx.ne.1) watercaptag(ig)=.true. -!! write(*,*) "outliers ", lati(ig)*180./pi, -!! . long(ig)*180./pi -!! !dryness(ig) = 1. -!! alternate = 1 -!! else -!! alternate = 0 -!! endif !end if alternate = 0 -!! -!! else -!! -!! if (ngridmx.ne.1) watercaptag(ig)=.true. -!! write(*,*) "outliers ", lati(ig)*180./pi, -!! . long(ig)*180./pi -!! -!! endif ! end if temptag -!! -!! endif -!! -!! -!!c Opposite olympia planitia water cap -!!c---------------- proposition par AS -------------------- -!!c-------------------------------------------------------- -!!c if ( ( lati(ig)*180./pi .ge. 82 ) .and. -!!c . ( lati(ig)*180./pi .le. 84 ) .and. -!!c . ( ( long(ig)*180./pi .gt. -030. ) .and. -!!c . ( long(ig)*180./pi .lt. 090. ) ) ) then -!!c---------------- proposition par TN -------------------- -!!c-------------------------------------------------------- -!! if ( ( lati(ig)*180./pi .ge. 80 ) .and. -!! . ( lati(ig)*180./pi .le. 84 ) .and. -!! . ( ( long(ig)*180./pi .gt. -030. ) .and. -!! . ( long(ig)*180./pi .lt. 090. ) ) ) then -!! if (ngridmx.ne.1) then -!! watercaptag(ig)=.true. -!! write(*,*) "central cap add ", lati(ig)*180./pi, -!! . long(ig)*180./pi -!! endif -!! !dryness(ig) = 1. -!! endif + if ( ( ( lati(ig)*180./pi .ge. 77. ) .and. ! cap #2 + . ( lati(ig)*180./pi .le. 80. ) .and. + . ( long(ig)*180./pi .ge. 110. ) .and. + . ( long(ig)*180./pi .le. 181. ) ) + . .or. + + . ( ( lati(ig)*180./pi .ge. 75. ) .and. ! cap #4 (Korolev crater) + . ( lati(ig)*180./pi .le. 76. ) .and. + . ( long(ig)*180./pi .ge. 150. ) .and. + . ( long(ig)*180./pi .le. 168. ) ) + . .or. + . ( ( lati(ig)*180./pi .ge. 77 ) .and. ! cap #5 + . ( lati(ig)*180./pi .le. 80. ) .and. + . ( long(ig)*180./pi .ge. -150.) .and. + . ( long(ig)*180./pi .le. -110.) ) ) + . then + + if (temptag) then + + if ((alternate .eq. 0)) then ! 1/2 en 64x48 sinon trop large en lat + if (ngridmx.ne.1) watercaptag(ig)=.true. +! print*,'ice water cap', lati(ig)*180./pi, +! . long(ig)*180./pi, ig + !dryness(ig) = 1. + alternate = 1 + else + alternate = 0 + endif !end if alternate = 0 + + else + +! if (ngridmx.ne.1) watercaptag(ig)=.true. +! write(*,*) "outliers ", lati(ig)*180./pi, +! . long(ig)*180./pi + + endif ! end if temptag + + endif + + +c Opposite olympia planitia water cap +c---------------- proposition par AS -------------------- +c-------------------------------------------------------- +c if ( ( lati(ig)*180./pi .ge. 82 ) .and. +c . ( lati(ig)*180./pi .le. 84 ) .and. +c . ( ( long(ig)*180./pi .gt. -030. ) .and. +c . ( long(ig)*180./pi .lt. 090. ) ) ) then +c---------------- proposition par TN -------------------- +c-------------------------------------------------------- + if ( ( ( lati(ig)*180./pi .ge. 80 ) .and. + . ( lati(ig)*180./pi .le. 84 ) ) + . .and. + . ( ( long(ig)*180./pi .lt. -95. ) .or. !!! 32x24 + . ( long(ig)*180./pi .gt. 85. ) ) ) then !!! 32x24 +! . ( ( ( long(ig)*180./pi .ge. -29. ) .and. !!! 64x48 +! . ( long(ig)*180./pi .le. 90. ) ) .or. !!! 64x48 +! . ( ( long(ig)*180./pi .ge. -77. ) .and. !!! 64x48 +! . ( long(ig)*180./pi .le. -70. ) ) ) ) then !!! 64x48 + if (ngridmx.ne.1) then + watercaptag(ig)=.true. + print*,'ice water cap', lati(ig)*180./pi, + . long(ig)*180./pi, ig + endif + !dryness(ig) = 1. + endif c Central cap @@ -159,6 +270,6 @@ if (lati(ig)*180./pi.gt.84) then - PRINT*,'central cap', lati(ig)*180./pi, - . long(ig)*180./pi + print*,'ice water cap', lati(ig)*180./pi, + . long(ig)*180./pi, ig if (ngridmx.ne.1) watercaptag(ig)=.true. !dryness(ig) = 1. @@ -173,7 +284,11 @@ c endif endif ! (lati>80 deg) + + end do ! (ngridmx) + + endif ! of if (improvedicelocation) #endif - end do ! (ngridmx) - ENDIF ! (caps & water) + + ENDIF ! (caps & water) c =============================================================== @@ -210,4 +325,5 @@ PRINT*,'maximum albedo avec water caps', s psolaralb(ISMAX(ngrid,psolaralb,1),1) + ENDIF Index: trunk/LMDZ.MARS/libf/phymars/testphys1d.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/testphys1d.F (revision 517) +++ trunk/LMDZ.MARS/libf/phymars/testphys1d.F (revision 520) @@ -581,5 +581,11 @@ $ thermo,qsurf) endif - watercaptag(ngridmx)=.false. + +c Regarder si le sol est un reservoir de glace d'eau +c -------------------------------------------------- + watercaptag(ngridmx)=.false. ! Par defaut il n'y pas de glace au sol + print *,'Water ice cap on ground ?' + call getin("watercaptag",watercaptag) + write(*,*) " watercaptag = ",watercaptag Index: trunk/LMDZ.MARS/libf/phymars/vdifc.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/vdifc.F (revision 517) +++ trunk/LMDZ.MARS/libf/phymars/vdifc.F (revision 520) @@ -661,9 +661,10 @@ else if (doubleq) then do ig=1,ngrid - !!! soulevement constant + if(co2ice(ig).lt.1) then ! soulevement pas constant pdqsdif(ig,igcm_dust_mass) = & -alpha_lift(igcm_dust_mass) pdqsdif(ig,igcm_dust_number) = - & -alpha_lift(igcm_dust_number) + & -alpha_lift(igcm_dust_number) + end if end do else if (submicron) then