Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/compile =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/compile (revision 234) +++ (revision ) @@ -1,18 +1,0 @@ -#! /bin/bash - -\rm fix_no_info.inc 2> /dev/null -\rm *.o 2> /dev/null -touch fix_no_info.inc - -pgf90 -byteswapio readmeteo.F90 \ --L$NETCDF/lib -lnetcdf \ --I$NETCDF/include \ --o readmeteo.exe - -pgf90 create_readmeteo.F90 \ --L$NETCDF/lib -lnetcdf \ --I$NETCDF/include \ --o create_readmeteo.exe - -\rm fix_no_info.inc 2> /dev/null -\rm *.o 2> /dev/null Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/compile_and_exec =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/compile_and_exec (revision 234) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/compile_and_exec (revision 235) @@ -8,3 +8,2 @@ echo 1 | create_readmeteo.exe readmeteo.exe < readmeteo.def - Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/compile_ifort =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/compile_ifort (revision 234) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/compile_ifort (revision 235) @@ -18,7 +18,6 @@ \rm *.o 2> /dev/null -echo I assume input_diagfi is correctly linked -echo I assume WPSFEED folder exists or is linked - -echo 1 | create_readmeteo.exe -readmeteo.exe < readmeteo.def +#echo I assume input_diagfi is correctly linked +#echo I assume WPSFEED folder exists or is linked +#echo 1 | create_readmeteo.exe +#readmeteo.exe < readmeteo.def Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/compile_pgf =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/compile_pgf (revision 235) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/compile_pgf (revision 235) @@ -0,0 +1,18 @@ +#! /bin/bash + +\rm fix_no_info.inc 2> /dev/null +\rm *.o 2> /dev/null +touch fix_no_info.inc + +pgf90 -byteswapio readmeteo.F90 \ +-L$NETCDF/lib -lnetcdf \ +-I$NETCDF/include \ +-o readmeteo.exe + +pgf90 create_readmeteo.F90 \ +-L$NETCDF/lib -lnetcdf \ +-I$NETCDF/include \ +-o create_readmeteo.exe + +\rm fix_no_info.inc 2> /dev/null +\rm *.o 2> /dev/null Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/inifis.F =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/inifis.F (revision 235) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/inifis.F (revision 235) @@ -0,0 +1,761 @@ + SUBROUTINE inifis(ngrid,nlayer, + $ wday_ini,wdaysec, + $ wappel_phys, + $ plat,plon,parea, + $ prad,pg,pr,pcpp, + $ nq, wdt, + $ womeg,wmugaz, + $ wyear_day,wperiheli,waphelie,wperi_day,wobliquit, + $ wz0,wemin_turb,wlmixmin, + $ wemissiv,wemissiceN,wemissiceS,walbediceN,walbediceS, + $ wiceradiusN,wiceradiusS,wdtemisiceN,wdtemisiceS, + $ walbedodat,winertiedat,wphisfi, + $ wzmea,wzstd,wzsig,wzgam,wzthe, + $ wtheta,wpsi) + + + IMPLICIT NONE +c======================================================================= +c +c CAREFUL: THIS IS A VERSION TO BE USED WITH WRF !!! +c +c ... CHECK THE ****WRF lines +c +c======================================================================= +c +c subject: +c -------- +c +c Initialisation for the physical parametrisations of the LMD +c martian atmospheric general circulation modele. +c +c author: Frederic Hourdin 15 / 10 /93 +c ------- +c modified: Sebastien Lebonnois 11/06/2003 (new callphys.def) +c adapted to the WRF use - Aymeric Spiga - Jan 2007 +c +c arguments: +c ---------- +c +c input: +c ------ +c +c ngrid Size of the horizontal grid. +c All internal loops are performed on that grid. +c nlayer Number of vertical layers. +c pdayref Day of reference for the simulation +c firstcall True at the first call +c lastcall True at the last call +c pday Number of days counted from the North. Spring +c equinoxe. +c +c======================================================================= +c +c----------------------------------------------------------------------- +c declarations: +c ------------- + +#include "dimensions.h" +#include "dimphys.h" +#include "planete.h" +#include "comcstfi.h" +#include "comsaison.h" +#include "comdiurn.h" +#include "comgeomfi.h" +#include "callkeys.h" +#include "surfdat.h" +#include "slope.h" + + + INTEGER ngrid,nlayer,nq + + REAL prad,pg,pr,pcpp,pdaysec + REAL womeg,wmugaz,wdaysec + REAL wyear_day,wperiheli,waphelie,wperi_day,wobliquit + REAL wz0,wemin_turb,wlmixmin + REAL wemissiv,wemissiceN,wemissiceS,walbediceN,walbediceS + + REAL wiceradiusN,wiceradiusS,wdtemisiceN,wdtemisiceS + REAL walbedodat(ngrid),winertiedat(ngrid),wphisfi(ngrid) + REAL wzmea(ngrid),wzstd(ngrid),wzsig(ngrid) + REAL wzgam(ngrid),wzthe(ngrid) + REAL wtheta(ngrid),wpsi(ngrid) + + REAL plat(ngrid),plon(ngrid),parea(ngridmx) + integer wday_ini + REAL wdt + INTEGER ig,ierr + + INTEGER wecri_phys, wappel_phys + + EXTERNAL iniorbit,orbite + EXTERNAL SSUM + REAL SSUM + + CHARACTER ch1*12 + CHARACTER ch80*80 + + logical chem, h2o + +c ****WRF +c +c------------------------------------------------------ +c Fill some parameters in the 'include' files +c >> Do part of the job previously done by phyetat0.F +c >> Complete list of parameters is found in tabfi.F +c------------------------------------------------------ +c +c Values are defined in the module_model_constants.F WRF routine +c + ! in 'comcstfi.h' + rad=prad + cpp=pcpp + g=pg + r=pr + rcp=r/cpp + daysec=wdaysec + omeg=womeg + mugaz=wmugaz + print*,"check: rad,cpp,g,r,rcp,daysec,omeg,mugaz" + print*,rad,cpp,g,r,rcp,daysec,omeg,mugaz + + ! in 'planet.h' + year_day=wyear_day + periheli=wperiheli + aphelie=waphelie + peri_day=wperi_day + obliquit=wobliquit + z0=wz0 + emin_turb=wemin_turb + lmixmin=wlmixmin + print*,"check: year_day,periheli,aphelie,peri_day,obliquit" + print*,year_day,periheli,aphelie,peri_day,obliquit + print*,"check: z0,emin_turb,lmixmin" + print*,z0,emin_turb,lmixmin + + ! in 'surfdat.h' + emissiv=wemissiv + emisice(1)=wemissiceN + emisice(2)=wemissiceS + albedice(1)=walbediceN + albedice(2)=walbediceS + iceradius(1)=wiceradiusN + iceradius(2)=wiceradiusS + dtemisice(1)=wdtemisiceN + dtemisice(2)=wdtemisiceS + print*,"check: emissiv,emisice,albedice,iceradius,dtemisice" + print*,emissiv,emisice,albedice,iceradius,dtemisice + +c +c Values are defined in the WPS processing +c + albedodat(:)=walbedodat(:) + inertiedat(:)=winertiedat(:) + phisfi(:)=wphisfi(:) + print*,"check: albedodat(1),inertiedat(1),phisfi(1)" + print*,albedodat(1),inertiedat(1),phisfi(1) + print*,"check: albedodat(end),inertiedat(end),phisfi(end)" + print*,albedodat(ngrid),inertiedat(ngrid),phisfi(ngrid) + + ! NB: usually, gravity wave scheme is useless in mesoscale modeling + ! NB: we however keep the option for coarse grid case ... + zmea(:)=wzmea(:) + zstd(:)=wzstd(:) + zsig(:)=wzsig(:) + zgam(:)=wzgam(:) + zthe(:)=wzthe(:) + print*,"check: gw param" + print*,zmea(1),zmea(ngrid) + print*,zstd(1),zstd(ngrid) + print*,zsig(1),zsig(ngrid) + print*,zgam(1),zgam(ngrid) + print*,zthe(1),zthe(ngrid) + + ! + ! in slope.h + ! + theta_sl(:)=wtheta(:) + psi_sl(:)=wpsi(:) + print*,"check: theta_sl(1),psi_sl(1)" + print*,theta_sl(1),psi_sl(1) + print*,"check: theta_sl(end),psi_sl(end)" + print*,theta_sl(ngrid),psi_sl(ngrid) + + +c ****WRF + + + +c -------------------------------------------------------- +c The usual Tests +c -------------- + +c ****WRF +c useless here, because it was already done in the WRF driver + + IF (nlayer.NE.nlayermx) THEN + PRINT*,'STOP in inifis' + PRINT*,'Probleme de dimensions :' + PRINT*,'nlayer = ',nlayer + PRINT*,'nlayermx = ',nlayermx + STOP + ENDIF + + IF (ngrid.NE.ngridmx) THEN + PRINT*,'STOP in inifis' + PRINT*,'Probleme de dimensions :' + PRINT*,'ngrid = ',ngrid + PRINT*,'ngridmx = ',ngridmx + STOP + ENDIF + + +c -------------------------------------------------------------- +c Reading the "callphys.def" file controlling some key options +c -------------------------------------------------------------- + + callrad=.true. + calldifv=.true. + calladj=.true. + callcond=.true. + callsoil=.true. + season=.true. + diurnal=.false. + lwrite=.false. + calllott=.true. + iaervar=2 + iddist=3 + topdustref=55. + OPEN(99,file='callphys.def',status='old',form='formatted' + . ,iostat=ierr) + IF(ierr.EQ.0) THEN + !PRINT* + !PRINT* + PRINT*,'--------------------------------------------' + PRINT*,' Parametres pour la physique (callphys.def)' + PRINT*,'--------------------------------------------' + + READ(99,*) + READ(99,*) + + READ(99,fmt='(a)') ch1 + READ(99,*) tracer + WRITE(*,8000) ch1,tracer + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') diurnal + WRITE(*,8000) ch1,diurnal + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') season + WRITE(*,8000) ch1,season + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') lwrite + WRITE(*,8000) ch1,lwrite + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') callstats + WRITE(*,8000) ch1,callstats + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') calleofdump + WRITE(*,8000) ch1,calleofdump + + READ(99,*) + READ(99,*) + + READ(99,fmt='(a)') ch1 + READ(99,*,iostat=ierr) iaervar + if(ierr.ne.0) stop'no iaervar in callphys.def (old?)' +c****WRF: ligne trop longue .... + WRITE(*,8001) ch1,iaervar + + READ(99,fmt='(a)') ch1 + READ(99,*) iddist + WRITE(*,8001) ch1,iddist + + READ(99,fmt='(a)') ch1 + READ(99,*) topdustref + WRITE(*,8002) ch1,topdustref + + READ(99,*) + READ(99,*) + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') callrad + WRITE(*,8000) ch1,callrad + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') callnlte + WRITE(*,8000) ch1,callnlte + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') callnirco2 + WRITE(*,8000) ch1,callnirco2 + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') calldifv + WRITE(*,8000) ch1,calldifv + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') calladj + WRITE(*,8000) ch1,calladj + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') callcond + WRITE(*,8000) ch1,callcond + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') callsoil + WRITE(*,8000) ch1,callsoil + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') calllott + WRITE(*,8000) ch1,calllott + + READ(99,*) + READ(99,*) + + READ(99,fmt='(a)') ch1 + READ(99,*) iradia + WRITE(*,8001) ch1,iradia + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') callg2d + WRITE(*,8000) ch1,callg2d + + READ(99,fmt='(a)') ch1 + READ(99,*) rayleigh + WRITE(*,8000) ch1,rayleigh + + READ(99,*) + READ(99,*) + +c TRACERS: + + READ(99,fmt='(a)') ch1 + READ(99,*) dustbin + WRITE(*,8001) ch1,dustbin + + READ(99,fmt='(a)') ch1 + READ(99,*) active + WRITE(*,8000) ch1,active + +c Test of incompatibility: +c if active is used, then dustbin should be > 0 + + if (active.and.(dustbin.lt.1)) then + print*,'if active is used, then dustbin should > 0' + stop + endif + + READ(99,fmt='(a)') ch1 + READ(99,*) doubleq + WRITE(*,8000) ch1,doubleq + +c Test of incompatibility: +c if doubleq is used, then dustbin should be 1 + + if (doubleq.and.(dustbin.ne.1)) then + print*,'if doubleq is used, then dustbin should be 1' + stop + endif + + READ(99,fmt='(a)') ch1 + READ(99,*) lifting + WRITE(*,8000) ch1,lifting + +c Test of incompatibility: +c if lifting is used, then dustbin should be > 0 + + if (lifting.and.(dustbin.lt.1)) then + print*,'if lifting is used, then dustbin should > 0' + stop + endif + + READ(99,fmt='(a)') ch1 + READ(99,*) callddevil + WRITE(*,8000) ch1,callddevil + +c Test of incompatibility: +c if dustdevil is used, then dustbin should be > 0 + + if (callddevil.and.(dustbin.lt.1)) then + print*,'if dustdevil is used, then dustbin should > 0' + stop + endif + + READ(99,fmt='(a)') ch1 + READ(99,*) scavenging + WRITE(*,8000) ch1,scavenging + +c Test of incompatibility: +c if scavenging is used, then dustbin should be > 0 + + if (scavenging.and.(dustbin.lt.1)) then + print*,'if scavenging is used, then dustbin should > 0' + stop + endif + + READ(99,fmt='(a)') ch1 + READ(99,*) sedimentation + WRITE(*,8000) ch1,sedimentation + + READ(99,fmt='(a)') ch1 + READ(99,*) iceparty + WRITE(*,8000) ch1,iceparty + + READ(99,fmt='(a)') ch1 + READ(99,*) activice + WRITE(*,8000) ch1,activice + +c Test of incompatibility: +c if activice is used, then iceparty should be used too + + if (activice.and..not.iceparty) then + print*,'if activice is used, iceparty should be used too' + stop + endif + + READ(99,fmt='(a)') ch1 + READ(99,*) water + WRITE(*,8000) ch1,water + +c Test of incompatibility: +c if iceparty is used, then water should be used too + + if (.not.water) then + iceparty = .false. + endif + + READ(99,fmt='(a)') ch1 + READ(99,*) caps + WRITE(*,8000) ch1,caps + + READ(99,fmt='(a)') ch1 + READ(99,*) photochem + WRITE(*,8000) ch1,photochem + + READ(99,*) + READ(99,*) + +c THERMOSPHERE + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') callthermos + WRITE(*,8000) ch1,callthermos + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') thermoswater + WRITE(*,8000) ch1,thermoswater + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') callconduct + WRITE(*,8000) ch1,callconduct + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') calleuv + WRITE(*,8000) ch1,calleuv + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') callmolvis + WRITE(*,8000) ch1,callmolvis + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') callmoldiff + WRITE(*,8000) ch1,callmoldiff + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') thermochem + WRITE(*,8000) ch1,thermochem + + READ(99,fmt='(a)') ch1 + READ(99,*) solarcondate + WRITE(*,*) ch1,solarcondate + + if (.not.callthermos) then + thermoswater=.false. + callconduct=.false. + calleuv=.false. + callmolvis=.false. + callmoldiff=.false. + thermochem=.false. + end if + +c Test of incompatibility: +c if photochem is used, then water should be used too + + if (photochem.and..not.water) then + print*,'if photochem is used, water should be used too' + stop + endif + +c if callthermos is used, then thermoswater should be used too +c (if water not used already) + + if (callthermos .and. .not.water) then + if (callthermos .and. .not.thermoswater) then + print*,'if callthermos is used, water or thermoswater + & should be used too' + stop + endif + endif + + PRINT*,'--------------------------------------------' + PRINT* + PRINT* + ELSE + write(*,*) + write(*,*) 'Cannot read file callphys.def. Is it here ?' + stop + ENDIF + CLOSE(99) + +8000 FORMAT(t5,a12,l8) +8001 FORMAT(t5,a12,i8) +8002 FORMAT(t5,a12,f8.1) + + +c ****WRF +c***************************************************** +c Since it comes from WRF settings, we have to +c fill dtphys in the include file +c It must be set now, because it is used afterwards +c +c Opportunity is taken to fill ecri_phys as well +c***************************************************** + dtphys=wdt*float(wappel_phys) + print*,'Physical timestep (s) ',dtphys + ecri_phys=8.e18 !! a dummy low frequency + print*,'Physical frequency for writing ',ecri_phys +c +c ****WRF + + + PRINT* + PRINT*,'daysec',daysec + PRINT* + PRINT*,'The radiative transfer is computed:' + PRINT*,' each ',iradia,' physical time-step' + PRINT*,' or each ',iradia*dtphys,' seconds' + PRINT* + + + +c -------------------------------------------------------------- +c Managing the Longwave radiative transfer +c -------------------------------------------------------------- + +c In most cases, the run just use the following values : +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + callemis=.true. +c ilwd=10*int(daysec/dtphys) ! bug before 22/10/01 + ilwd=10*int(daysec/dtphys) + ilwn=2 + linear=.true. + ncouche=3 + alphan=0.4 + ilwb=2 + semi=0 + +c BUT people working hard on the LW may want to read them in 'radia.def' +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + OPEN(99,file='radia.def',status='old',form='formatted' + . ,iostat=ierr) + IF(ierr.EQ.0) THEN + write(*,*) 'Reading radia.def !!!' + READ(99,fmt='(a)') ch1 + READ(99,*) callemis + WRITE(*,8000) ch1,callemis + + READ(99,fmt='(a)') ch1 + READ(99,*) iradia + WRITE(*,8001) ch1,iradia + + READ(99,fmt='(a)') ch1 + READ(99,*) ilwd + WRITE(*,8001) ch1,ilwd + + READ(99,fmt='(a)') ch1 + READ(99,*) ilwn + WRITE(*,8001) ch1,ilwn + + READ(99,fmt='(a)') ch1 + READ(99,*) linear + WRITE(*,8000) ch1,linear + + READ(99,fmt='(a)') ch1 + READ(99,*) ncouche + WRITE(*,8001) ch1,ncouche + + READ(99,fmt='(a)') ch1 + READ(99,*) alphan + WRITE(*,*) ch1,alphan + + READ(99,fmt='(a)') ch1 + READ(99,*) ilwb + WRITE(*,8001) ch1,ilwb + + + READ(99,fmt='(a)') ch1 + READ(99,'(l1)') callg2d + WRITE(*,8000) ch1,callg2d + + READ(99,fmt='(a)') ch1 + READ(99,*) semi + WRITE(*,*) ch1,semi + end if + +c----------------------------------------------------------------------- +c Some more initialization: +c ------------------------ + + ! in 'comgeomfi.h' + CALL SCOPY(ngrid,plon,1,long,1) + CALL SCOPY(ngrid,plat,1,lati,1) + CALL SCOPY(ngrid,parea,1,area,1) + totarea=SSUM(ngridmx,area,1) + + ! in 'comdiurn.h' + DO ig=1,ngrid + sinlat(ig)=sin(plat(ig)) + coslat(ig)=cos(plat(ig)) + sinlon(ig)=sin(plon(ig)) + coslon(ig)=cos(plon(ig)) + ENDDO + + pi=2.*asin(1.) + +c managing the tracers, and tests: +c ------------------------------- + + if(tracer) then + +c when photochem is used, nqchem_min is the rank +c of the first chemical species + + nqchem_min = 1 + if (photochem .or. callthermos) then + chem = .true. + if (doubleq) then + nqchem_min = 3 + else + nqchem_min = dustbin+1 + end if + end if + + if (water .or. thermoswater) h2o = .true. + +c TESTS + + print*,'TRACERS:' + + if ((doubleq).and.(h2o).and. + $ (chem).and.(iceparty)) then + print*,' 1: dust ; 2: dust (doubleq)' + print*,' 3 to ',nqmx-2,': chemistry' + print*,nqmx-1,': water ice ; ',nqmx,': water vapor' + endif + + if ((doubleq).and.(h2o).and. + $ (chem).and..not.(iceparty)) then + print*,' 1: dust ; 2: dust (doubleq)' + print*,' 3 to ',nqmx-1,': chemistry' + print*,nqmx,': water vapor' + endif + + if ((doubleq).and.(h2o).and. + $ .not.(chem).and.(iceparty)) then + print*,' 1: dust ; 2: dust (doubleq)' + print*,nqmx-1,': water ice ; ',nqmx,': water vapor' + if (nqmx.ne.4) then + print*,'nqmx should be 4 with these options.' + print*,'(or check callphys.def)' + stop + endif + endif + + if ((doubleq).and.(h2o).and. + $ .not.(chem).and..not.(iceparty)) then + print*,' 1: dust ; 2: dust (doubleq)' + print*,nqmx,': water vapor' + if (nqmx.ne.3) then + print*,'nqmx should be 3 with these options...' + print*,'(or check callphys.def)' + stop + endif + endif + + if ((doubleq).and..not.(h2o)) then + print*,' 1: dust ; 2: dust (doubleq)' + if (nqmx.ne.2) then + print*,'nqmx should be 2 with these options...' + print*,'(or check callphys.def)' + stop + endif + endif + + if (.not.(doubleq).and.(h2o).and. + $ (chem).and.(iceparty)) then + if (dustbin.gt.0) then + print*,' 1 to ',dustbin,': dust bins' + endif + print*,nqchem_min,' to ',nqmx-2,': chemistry' + print*,nqmx-1,': water ice ; ',nqmx,': water vapor' + endif + if (.not.(doubleq).and.(h2o).and. + $ (chem).and..not.(iceparty)) then + if (dustbin.gt.0) then + print*,' 1 to ',dustbin,': dust bins' + endif + print*,nqchem_min,' to ',nqmx-1,': chemistry' + print*,nqmx,': water vapor' + endif + if (.not.(doubleq).and.(h2o).and. + $ .not.(chem).and.(iceparty)) then + if (dustbin.gt.0) then + print*,' 1 to ',dustbin,': dust bins' + endif + print*,nqmx-1,': water ice ; ',nqmx,': water vapor' + if (nqmx.ne.(dustbin+2)) then + print*,'nqmx should be ',(dustbin+2), + $ ' with these options...' + print*,'(or check callphys.def)' + stop + endif + endif + if (.not.(doubleq).and.(h2o).and. + $ .not.(chem).and..not.(iceparty)) then + if (dustbin.gt.0) then + print*,' 1 to ',dustbin,': dust bins' + endif + print*,nqmx,': water vapor' + if (nqmx.ne.(dustbin+1)) then + print*,'nqmx should be ',(dustbin+1), + $ ' with these options...' + print*,'(or check callphys.def)' + stop + endif + endif + if (.not.(doubleq).and..not.(h2o)) then + if (dustbin.gt.0) then + print*,' 1 to ',dustbin,': dust bins' + if (nqmx.ne.dustbin) then + print*,'nqmx should be ',dustbin, + $ ' with these options...' + print*,'(or check callphys.def)' + stop + endif + else + print*,'dustbin=',dustbin, + $ ': tracer should be F with these options...' + $ ,'UNLESS you just want to move tracers around ' + endif + endif + + endif + + RETURN + END Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/meso_inifis.F =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/meso_inifis.F (revision 234) +++ (revision ) @@ -1,761 +1,0 @@ - SUBROUTINE meso_inifis(ngrid,nlayer,nq, - $ wdt, - $ wday_ini,wdaysec, - $ wappel_phys,wecri_phys, - $ plat,plon,parea, - $ prad,pg,pr,pcpp, - $ womeg,wmugaz, - $ wyear_day,wperiheli,waphelie,wperi_day,wobliquit, - $ wz0,wemin_turb,wlmixmin, - $ wemissiv,wemissiceN,wemissiceS,walbediceN,walbediceS, - $ wiceradiusN,wiceradiusS,wdtemisiceN,wdtemisiceS, - $ walbedodat,winertiedat,wphisfi, - $ wzmea,wzstd,wzsig,wzgam,wzthe, - $ wtheta,wpsi) - - - IMPLICIT NONE -c======================================================================= -c -c CAREFUL: THIS IS A VERSION TO BE USED WITH WRF !!! -c -c ... CHECK THE ****WRF lines -c -c======================================================================= -c -c subject: -c -------- -c -c Initialisation for the physical parametrisations of the LMD -c martian atmospheric general circulation modele. -c -c author: Frederic Hourdin 15 / 10 /93 -c ------- -c modified: Sebastien Lebonnois 11/06/2003 (new callphys.def) -c adapted to the WRF use - Aymeric Spiga - Jan 2007 -c -c arguments: -c ---------- -c -c input: -c ------ -c -c ngrid Size of the horizontal grid. -c All internal loops are performed on that grid. -c nlayer Number of vertical layers. -c pdayref Day of reference for the simulation -c firstcall True at the first call -c lastcall True at the last call -c pday Number of days counted from the North. Spring -c equinoxe. -c -c======================================================================= -c -c----------------------------------------------------------------------- -c declarations: -c ------------- - -#include "dimensions.h" -#include "dimphys.h" -#include "planete.h" -#include "comcstfi.h" -#include "comsaison.h" -#include "comdiurn.h" -#include "comgeomfi.h" -#include "callkeys.h" -#include "surfdat.h" -#include "slope.h" - - - INTEGER ngrid,nlayer,nq - - REAL prad,pg,pr,pcpp,pdaysec - REAL womeg,wmugaz,wdaysec - REAL wyear_day,wperiheli,waphelie,wperi_day,wobliquit - REAL wz0,wemin_turb,wlmixmin - REAL wemissiv,wemissiceN,wemissiceS,walbediceN,walbediceS - - REAL wiceradiusN,wiceradiusS,wdtemisiceN,wdtemisiceS - REAL walbedodat(ngrid),winertiedat(ngrid),wphisfi(ngrid) - REAL wzmea(ngrid),wzstd(ngrid),wzsig(ngrid) - REAL wzgam(ngrid),wzthe(ngrid) - REAL wtheta(ngrid),wpsi(ngrid) - - REAL plat(ngrid),plon(ngrid),parea(ngridmx) - integer wday_ini - REAL wdt - INTEGER ig,ierr - - INTEGER wecri_phys, wappel_phys - - EXTERNAL iniorbit,orbite - EXTERNAL SSUM - REAL SSUM - - CHARACTER ch1*12 - CHARACTER ch80*80 - - logical chem, h2o - -c ****WRF -c -c------------------------------------------------------ -c Fill some parameters in the 'include' files -c >> Do part of the job previously done by phyetat0.F -c >> Complete list of parameters is found in tabfi.F -c------------------------------------------------------ -c -c Values are defined in the module_model_constants.F WRF routine -c - ! in 'comcstfi.h' - rad=prad - cpp=pcpp - g=pg - r=pr - rcp=r/cpp - daysec=wdaysec - omeg=womeg - mugaz=wmugaz - print*,"check: rad,cpp,g,r,rcp,daysec,omeg,mugaz" - print*,rad,cpp,g,r,rcp,daysec,omeg,mugaz - - ! in 'planet.h' - year_day=wyear_day - periheli=wperiheli - aphelie=waphelie - peri_day=wperi_day - obliquit=wobliquit - z0=wz0 - emin_turb=wemin_turb - lmixmin=wlmixmin - print*,"check: year_day,periheli,aphelie,peri_day,obliquit" - print*,year_day,periheli,aphelie,peri_day,obliquit - print*,"check: z0,emin_turb,lmixmin" - print*,z0,emin_turb,lmixmin - - ! in 'surfdat.h' - emissiv=wemissiv - emisice(1)=wemissiceN - emisice(2)=wemissiceS - albedice(1)=walbediceN - albedice(2)=walbediceS - iceradius(1)=wiceradiusN - iceradius(2)=wiceradiusS - dtemisice(1)=wdtemisiceN - dtemisice(2)=wdtemisiceS - print*,"check: emissiv,emisice,albedice,iceradius,dtemisice" - print*,emissiv,emisice,albedice,iceradius,dtemisice - -c -c Values are defined in the WPS processing -c - albedodat(:)=walbedodat(:) - inertiedat(:)=winertiedat(:) - phisfi(:)=wphisfi(:) - print*,"check: albedodat(1),inertiedat(1),phisfi(1)" - print*,albedodat(1),inertiedat(1),phisfi(1) - print*,"check: albedodat(end),inertiedat(end),phisfi(end)" - print*,albedodat(ngrid),inertiedat(ngrid),phisfi(ngrid) - - ! NB: usually, gravity wave scheme is useless in mesoscale modeling - ! NB: we however keep the option for coarse grid case ... - zmea(:)=wzmea(:) - zstd(:)=wzstd(:) - zsig(:)=wzsig(:) - zgam(:)=wzgam(:) - zthe(:)=wzthe(:) - print*,"check: gw param" - print*,zmea(1),zmea(ngrid) - print*,zstd(1),zstd(ngrid) - print*,zsig(1),zsig(ngrid) - print*,zgam(1),zgam(ngrid) - print*,zthe(1),zthe(ngrid) - - ! - ! in slope.h - ! - theta_sl(:)=wtheta(:) - psi_sl(:)=wpsi(:) - print*,"check: theta_sl(1),psi_sl(1)" - print*,theta_sl(1),psi_sl(1) - print*,"check: theta_sl(end),psi_sl(end)" - print*,theta_sl(ngrid),psi_sl(ngrid) - - -c ****WRF - - - -c -------------------------------------------------------- -c The usual Tests -c -------------- - -c ****WRF -c useless here, because it was already done in the WRF driver - - IF (nlayer.NE.nlayermx) THEN - PRINT*,'STOP in inifis' - PRINT*,'Probleme de dimensions :' - PRINT*,'nlayer = ',nlayer - PRINT*,'nlayermx = ',nlayermx - STOP - ENDIF - - IF (ngrid.NE.ngridmx) THEN - PRINT*,'STOP in inifis' - PRINT*,'Probleme de dimensions :' - PRINT*,'ngrid = ',ngrid - PRINT*,'ngridmx = ',ngridmx - STOP - ENDIF - - -c -------------------------------------------------------------- -c Reading the "callphys.def" file controlling some key options -c -------------------------------------------------------------- - - callrad=.true. - calldifv=.true. - calladj=.true. - callcond=.true. - callsoil=.true. - season=.true. - diurnal=.false. - lwrite=.false. - calllott=.true. - iaervar=2 - iddist=3 - topdustref=55. - OPEN(99,file='callphys.def',status='old',form='formatted' - . ,iostat=ierr) - IF(ierr.EQ.0) THEN - !PRINT* - !PRINT* - PRINT*,'--------------------------------------------' - PRINT*,' Parametres pour la physique (callphys.def)' - PRINT*,'--------------------------------------------' - - READ(99,*) - READ(99,*) - - READ(99,fmt='(a)') ch1 - READ(99,*) tracer - WRITE(*,8000) ch1,tracer - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') diurnal - WRITE(*,8000) ch1,diurnal - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') season - WRITE(*,8000) ch1,season - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') lwrite - WRITE(*,8000) ch1,lwrite - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') callstats - WRITE(*,8000) ch1,callstats - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') calleofdump - WRITE(*,8000) ch1,calleofdump - - READ(99,*) - READ(99,*) - - READ(99,fmt='(a)') ch1 - READ(99,*,iostat=ierr) iaervar - if(ierr.ne.0) stop'no iaervar in callphys.def (old?)' -c****WRF: ligne trop longue .... - WRITE(*,8001) ch1,iaervar - - READ(99,fmt='(a)') ch1 - READ(99,*) iddist - WRITE(*,8001) ch1,iddist - - READ(99,fmt='(a)') ch1 - READ(99,*) topdustref - WRITE(*,8002) ch1,topdustref - - READ(99,*) - READ(99,*) - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') callrad - WRITE(*,8000) ch1,callrad - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') callnlte - WRITE(*,8000) ch1,callnlte - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') callnirco2 - WRITE(*,8000) ch1,callnirco2 - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') calldifv - WRITE(*,8000) ch1,calldifv - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') calladj - WRITE(*,8000) ch1,calladj - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') callcond - WRITE(*,8000) ch1,callcond - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') callsoil - WRITE(*,8000) ch1,callsoil - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') calllott - WRITE(*,8000) ch1,calllott - - READ(99,*) - READ(99,*) - - READ(99,fmt='(a)') ch1 - READ(99,*) iradia - WRITE(*,8001) ch1,iradia - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') callg2d - WRITE(*,8000) ch1,callg2d - - READ(99,fmt='(a)') ch1 - READ(99,*) rayleigh - WRITE(*,8000) ch1,rayleigh - - READ(99,*) - READ(99,*) - -c TRACERS: - - READ(99,fmt='(a)') ch1 - READ(99,*) dustbin - WRITE(*,8001) ch1,dustbin - - READ(99,fmt='(a)') ch1 - READ(99,*) active - WRITE(*,8000) ch1,active - -c Test of incompatibility: -c if active is used, then dustbin should be > 0 - - if (active.and.(dustbin.lt.1)) then - print*,'if active is used, then dustbin should > 0' - stop - endif - - READ(99,fmt='(a)') ch1 - READ(99,*) doubleq - WRITE(*,8000) ch1,doubleq - -c Test of incompatibility: -c if doubleq is used, then dustbin should be 1 - - if (doubleq.and.(dustbin.ne.1)) then - print*,'if doubleq is used, then dustbin should be 1' - stop - endif - - READ(99,fmt='(a)') ch1 - READ(99,*) lifting - WRITE(*,8000) ch1,lifting - -c Test of incompatibility: -c if lifting is used, then dustbin should be > 0 - - if (lifting.and.(dustbin.lt.1)) then - print*,'if lifting is used, then dustbin should > 0' - stop - endif - - READ(99,fmt='(a)') ch1 - READ(99,*) callddevil - WRITE(*,8000) ch1,callddevil - -c Test of incompatibility: -c if dustdevil is used, then dustbin should be > 0 - - if (callddevil.and.(dustbin.lt.1)) then - print*,'if dustdevil is used, then dustbin should > 0' - stop - endif - - READ(99,fmt='(a)') ch1 - READ(99,*) scavenging - WRITE(*,8000) ch1,scavenging - -c Test of incompatibility: -c if scavenging is used, then dustbin should be > 0 - - if (scavenging.and.(dustbin.lt.1)) then - print*,'if scavenging is used, then dustbin should > 0' - stop - endif - - READ(99,fmt='(a)') ch1 - READ(99,*) sedimentation - WRITE(*,8000) ch1,sedimentation - - READ(99,fmt='(a)') ch1 - READ(99,*) iceparty - WRITE(*,8000) ch1,iceparty - - READ(99,fmt='(a)') ch1 - READ(99,*) activice - WRITE(*,8000) ch1,activice - -c Test of incompatibility: -c if activice is used, then iceparty should be used too - - if (activice.and..not.iceparty) then - print*,'if activice is used, iceparty should be used too' - stop - endif - - READ(99,fmt='(a)') ch1 - READ(99,*) water - WRITE(*,8000) ch1,water - -c Test of incompatibility: -c if iceparty is used, then water should be used too - - if (.not.water) then - iceparty = .false. - endif - - READ(99,fmt='(a)') ch1 - READ(99,*) caps - WRITE(*,8000) ch1,caps - - READ(99,fmt='(a)') ch1 - READ(99,*) photochem - WRITE(*,8000) ch1,photochem - - READ(99,*) - READ(99,*) - -c THERMOSPHERE - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') callthermos - WRITE(*,8000) ch1,callthermos - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') thermoswater - WRITE(*,8000) ch1,thermoswater - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') callconduct - WRITE(*,8000) ch1,callconduct - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') calleuv - WRITE(*,8000) ch1,calleuv - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') callmolvis - WRITE(*,8000) ch1,callmolvis - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') callmoldiff - WRITE(*,8000) ch1,callmoldiff - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') thermochem - WRITE(*,8000) ch1,thermochem - - READ(99,fmt='(a)') ch1 - READ(99,*) solarcondate - WRITE(*,*) ch1,solarcondate - - if (.not.callthermos) then - thermoswater=.false. - callconduct=.false. - calleuv=.false. - callmolvis=.false. - callmoldiff=.false. - thermochem=.false. - end if - -c Test of incompatibility: -c if photochem is used, then water should be used too - - if (photochem.and..not.water) then - print*,'if photochem is used, water should be used too' - stop - endif - -c if callthermos is used, then thermoswater should be used too -c (if water not used already) - - if (callthermos .and. .not.water) then - if (callthermos .and. .not.thermoswater) then - print*,'if callthermos is used, water or thermoswater - & should be used too' - stop - endif - endif - - PRINT*,'--------------------------------------------' - PRINT* - PRINT* - ELSE - write(*,*) - write(*,*) 'Cannot read file callphys.def. Is it here ?' - stop - ENDIF - CLOSE(99) - -8000 FORMAT(t5,a12,l8) -8001 FORMAT(t5,a12,i8) -8002 FORMAT(t5,a12,f8.1) - - -c ****WRF -c***************************************************** -c Since it comes from WRF settings, we have to -c fill dtphys in the include file -c It must be set now, because it is used afterwards -c -c Opportunity is taken to fill ecri_phys as well -c***************************************************** - dtphys=wdt*float(wappel_phys) - print*,'Physical timestep (s) ',dtphys - ecri_phys=wecri_phys - print*,'Physical frequency for writing ',ecri_phys -c -c ****WRF - - - PRINT* - PRINT*,'daysec',daysec - PRINT* - PRINT*,'The radiative transfer is computed:' - PRINT*,' each ',iradia,' physical time-step' - PRINT*,' or each ',iradia*dtphys,' seconds' - PRINT* - - - -c -------------------------------------------------------------- -c Managing the Longwave radiative transfer -c -------------------------------------------------------------- - -c In most cases, the run just use the following values : -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - callemis=.true. -c ilwd=10*int(daysec/dtphys) ! bug before 22/10/01 - ilwd=10*int(daysec/dtphys) - ilwn=2 - linear=.true. - ncouche=3 - alphan=0.4 - ilwb=2 - semi=0 - -c BUT people working hard on the LW may want to read them in 'radia.def' -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - OPEN(99,file='radia.def',status='old',form='formatted' - . ,iostat=ierr) - IF(ierr.EQ.0) THEN - write(*,*) 'Reading radia.def !!!' - READ(99,fmt='(a)') ch1 - READ(99,*) callemis - WRITE(*,8000) ch1,callemis - - READ(99,fmt='(a)') ch1 - READ(99,*) iradia - WRITE(*,8001) ch1,iradia - - READ(99,fmt='(a)') ch1 - READ(99,*) ilwd - WRITE(*,8001) ch1,ilwd - - READ(99,fmt='(a)') ch1 - READ(99,*) ilwn - WRITE(*,8001) ch1,ilwn - - READ(99,fmt='(a)') ch1 - READ(99,*) linear - WRITE(*,8000) ch1,linear - - READ(99,fmt='(a)') ch1 - READ(99,*) ncouche - WRITE(*,8001) ch1,ncouche - - READ(99,fmt='(a)') ch1 - READ(99,*) alphan - WRITE(*,*) ch1,alphan - - READ(99,fmt='(a)') ch1 - READ(99,*) ilwb - WRITE(*,8001) ch1,ilwb - - - READ(99,fmt='(a)') ch1 - READ(99,'(l1)') callg2d - WRITE(*,8000) ch1,callg2d - - READ(99,fmt='(a)') ch1 - READ(99,*) semi - WRITE(*,*) ch1,semi - end if - -c----------------------------------------------------------------------- -c Some more initialization: -c ------------------------ - - ! in 'comgeomfi.h' - CALL SCOPY(ngrid,plon,1,long,1) - CALL SCOPY(ngrid,plat,1,lati,1) - CALL SCOPY(ngrid,parea,1,area,1) - totarea=SSUM(ngridmx,area,1) - - ! in 'comdiurn.h' - DO ig=1,ngrid - sinlat(ig)=sin(plat(ig)) - coslat(ig)=cos(plat(ig)) - sinlon(ig)=sin(plon(ig)) - coslon(ig)=cos(plon(ig)) - ENDDO - - pi=2.*asin(1.) - -c managing the tracers, and tests: -c ------------------------------- - - if(tracer) then - -c when photochem is used, nqchem_min is the rank -c of the first chemical species - - nqchem_min = 1 - if (photochem .or. callthermos) then - chem = .true. - if (doubleq) then - nqchem_min = 3 - else - nqchem_min = dustbin+1 - end if - end if - - if (water .or. thermoswater) h2o = .true. - -c TESTS - - print*,'TRACERS:' - - if ((doubleq).and.(h2o).and. - $ (chem).and.(iceparty)) then - print*,' 1: dust ; 2: dust (doubleq)' - print*,' 3 to ',nqmx-2,': chemistry' - print*,nqmx-1,': water ice ; ',nqmx,': water vapor' - endif - - if ((doubleq).and.(h2o).and. - $ (chem).and..not.(iceparty)) then - print*,' 1: dust ; 2: dust (doubleq)' - print*,' 3 to ',nqmx-1,': chemistry' - print*,nqmx,': water vapor' - endif - - if ((doubleq).and.(h2o).and. - $ .not.(chem).and.(iceparty)) then - print*,' 1: dust ; 2: dust (doubleq)' - print*,nqmx-1,': water ice ; ',nqmx,': water vapor' - if (nqmx.ne.4) then - print*,'nqmx should be 4 with these options.' - print*,'(or check callphys.def)' - stop - endif - endif - - if ((doubleq).and.(h2o).and. - $ .not.(chem).and..not.(iceparty)) then - print*,' 1: dust ; 2: dust (doubleq)' - print*,nqmx,': water vapor' - if (nqmx.ne.3) then - print*,'nqmx should be 3 with these options...' - print*,'(or check callphys.def)' - stop - endif - endif - - if ((doubleq).and..not.(h2o)) then - print*,' 1: dust ; 2: dust (doubleq)' - if (nqmx.ne.2) then - print*,'nqmx should be 2 with these options...' - print*,'(or check callphys.def)' - stop - endif - endif - - if (.not.(doubleq).and.(h2o).and. - $ (chem).and.(iceparty)) then - if (dustbin.gt.0) then - print*,' 1 to ',dustbin,': dust bins' - endif - print*,nqchem_min,' to ',nqmx-2,': chemistry' - print*,nqmx-1,': water ice ; ',nqmx,': water vapor' - endif - if (.not.(doubleq).and.(h2o).and. - $ (chem).and..not.(iceparty)) then - if (dustbin.gt.0) then - print*,' 1 to ',dustbin,': dust bins' - endif - print*,nqchem_min,' to ',nqmx-1,': chemistry' - print*,nqmx,': water vapor' - endif - if (.not.(doubleq).and.(h2o).and. - $ .not.(chem).and.(iceparty)) then - if (dustbin.gt.0) then - print*,' 1 to ',dustbin,': dust bins' - endif - print*,nqmx-1,': water ice ; ',nqmx,': water vapor' - if (nqmx.ne.(dustbin+2)) then - print*,'nqmx should be ',(dustbin+2), - $ ' with these options...' - print*,'(or check callphys.def)' - stop - endif - endif - if (.not.(doubleq).and.(h2o).and. - $ .not.(chem).and..not.(iceparty)) then - if (dustbin.gt.0) then - print*,' 1 to ',dustbin,': dust bins' - endif - print*,nqmx,': water vapor' - if (nqmx.ne.(dustbin+1)) then - print*,'nqmx should be ',(dustbin+1), - $ ' with these options...' - print*,'(or check callphys.def)' - stop - endif - endif - if (.not.(doubleq).and..not.(h2o)) then - if (dustbin.gt.0) then - print*,' 1 to ',dustbin,': dust bins' - if (nqmx.ne.dustbin) then - print*,'nqmx should be ',dustbin, - $ ' with these options...' - print*,'(or check callphys.def)' - stop - endif - else - print*,'dustbin=',dustbin, - $ ': tracer should be F with these options...' - $ ,'UNLESS you just want to move tracers around ' - endif - endif - - endif - - RETURN - END Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/meso_physiq.F =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/meso_physiq.F (revision 234) +++ (revision ) @@ -1,1572 +1,0 @@ - SUBROUTINE meso_physiq(ngrid,nlayer,nq, - $ firstcall,lastcall, - $ wday_ini, - $ pday,ptime,ptimestep, - $ pplev,pplay,pphi, - $ pu,pv,pt,pq, - $ pw, - $ pdu,pdv,pdt,pdq,pdpsrf,tracerdyn, - $ wtsurf,wtsoil,wemis,wq2,wqsurf,wco2ice, - $ wecritphys, - $ output_tab2d, output_tab3d, - $ flag_LES) - - - - IMPLICIT NONE -c======================================================================= -c -c CAREFUL: THIS IS A VERSION TO BE USED WITH WRF !!! -c -c ... CHECK THE ****WRF lines -c -c======================================================================= -c -c subject: -c -------- -c -c Organisation of the physical parametrisations of the LMD -c martian atmospheric general circulation model. -c -c The GCM can be run without or with tracer transport -c depending on the value of Logical "tracer" in file "callphys.def" -c Tracers may be water vapor, ice OR chemical species OR dust particles -c -c SEE comments in initracer.F about numbering of tracer species... -c -c It includes: -c -c 1. Initialisation: -c 1.5 Calculation of mean mass and cp, R and thermal conduction coeff. -c 2. Calcul of the radiative tendencies : radiative transfer -c (longwave and shortwave) for CO2 and dust. -c 3. Gravity wave and subgrid scale topography drag : -c 4. Vertical diffusion (turbulent mixing): -c 5. convective adjustment -c 6. TRACERS : -c 6a. water and water ice -c 6b. call for photochemistry when tracers are chemical species -c 6c. other scheme for tracer (dust) transport (lifting, sedimentation) -c 6d. updates (CO2 pressure variations, surface budget) -c 7. condensation and sublimation of carbon dioxide. -c 8. Surface and sub-surface temperature calculations -c 9. Writing output files : -c - "startfi", "histfi" (if it's time) -c - Saving statistics (if "callstats = .true.") -c - Dumping eof (if "calleofdump = .true.") -c - Output any needed variables in "diagfi" -c 10. Diagnostic: mass conservation of tracers -c -c author: -c ------- -c Frederic Hourdin 15/10/93 -c Francois Forget 1994 -c Christophe Hourdin 02/1997 -c Subroutine completly rewritten by F.Forget (01/2000) -c Introduction of the photochemical module: S. Lebonnois (11/2002) -c Introduction of the thermosphere module: M. Angelats i Coll (2002) -c Water ice clouds: Franck Montmessin (update 06/2003) -c WRF version: Aymeric Spiga (01-03/2007) -c -c -c -c arguments: -c ---------- -c -c input: -c ------ -c ecri period (in dynamical timestep) to write output -c ngrid Size of the horizontal grid. -c All internal loops are performed on that grid. -c nlayer Number of vertical layers. -c nq Number of advected fields -c firstcall True at the first call -c lastcall True at the last call -c pday Number of days counted from the North. Spring -c equinoxe. -c ptime Universal time (0 part of the job of phyetat0 is done in inifis -c > remaining initializations are passed here from the WRF variables -c > beware, some operations were done by phyetat0 (ex: tracers) -c > if any problems, look in phyetat0 -c - tsurf(:)=wtsurf(:) - PRINT*,'check: tsurf ',tsurf(1),tsurf(ngridmx) - tsoil(:,:)=wtsoil(:,:) - PRINT*,'check: tsoil ',tsoil(1,1),tsoil(ngridmx,nsoilmx) - emis(:)=wemis(:) - PRINT*,'check: emis ',emis(1),emis(ngridmx) - q2(:,:)=wq2(:,:) - PRINT*,'check: q2 ',q2(1,1),q2(ngridmx,nlayermx+1) - qsurf(:,:)=wqsurf(:,:) - PRINT*,'check: qsurf ',qsurf(1,1),qsurf(ngridmx,nqmx) - co2ice(:)=wco2ice(:) - PRINT*,'check: co2 ',co2ice(1),co2ice(ngridmx) - day_ini=wday_ini - -c artificially filling dyn3d/control.h is also required -c > iphysiq is put in WRF to be set easily (cf ptimestep) -c > day_step is simply deduced: -c - day_step=daysec/ptimestep - PRINT*,'Call to LMD physics:',day_step,' per Martian day' -c - iphysiq=ptimestep -c - ecritphy=wecritphys - PRINT*,'Write LMD physics each:',ecritphy,' seconds' - !!PRINT*,ecri_phys - !!PRINT*,float(ecri_phys) ... - !!renvoient tous deux des nombres absurdes - !!pourtant callkeys.h est inclus ... - !! - !!donc ecritphys est passe en argument ... - PRINT*,'Write LMD physics each:',ecritphy,' seconds' -c -c ****WRF -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - - - if (pday.ne.day_ini) then - write(*,*) "***ERROR Pb de synchronisation entre phys et dyn" - write(*,*) "jour dynamique: ",pday - write(*,*) "jour physique: ",day_ini - stop - endif - - write (*,*) 'In physic day_ini =', day_ini - - - -c Initialize albedo and orbital calculation -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - CALL surfini(ngrid,co2ice,qsurf,albedo) - CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit) - -c initialisation soil -c ~~~~~~~~~~~~~~~~~~~ - IF (callsoil) THEN - CALL soil(ngrid,nsoilmx,firstcall,inertiedat, - s ptimestep,tsurf,tsoil,capcal,fluxgrd) - ELSE - PRINT*,'WARNING! Thermal conduction in the soil turned off' - DO ig=1,ngrid - capcal(ig)=1.e5 - fluxgrd(ig)=0. - ENDDO - ENDIF - icount=1 - - -c initialisation pour les traceurs -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - tracerdyn=tracer - IF (tracer) THEN - CALL initracer(qsurf,co2ice) - ENDIF ! end tracer - - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c ****WRF -c -c nosense in mesoscale modeling -c -cc Determining gridpoint near Viking Lander 1 (used for diagnostic only) -cc ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -c if(ngrid.ne.1) then -c latvl1= 22.27 -c lonvl1= -47.94 -c ig_vl1= 1+ int( (1.5-(latvl1-90.)*jjm/180.) -2 )*iim + -c & int(1.5+(lonvl1+180)*iim/360.) -c write(*,*) 'Viking Lander 1 GCM point: lat,lon', -c & lati(ig_vl1)*180/pi, long(ig_vl1)*180/pi -c end if -c -c ****WRF -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - - -c Initializing thermospheric parameters -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - if (callthermos) call param_read - -c Initializing R and Cp as constant - - if (.not.callthermos .and. .not.photochem) then - do l=1,nlayermx - do ig=1,ngridmx - rnew(ig,l)=r - cpnew(ig,l)=cpp - mmean(ig,l)=mugaz - enddo - enddo - endif - - ENDIF ! (end of "if firstcall") - -!!!!!!!!!!!!!!!!TEST TEST -! IF (nocalculphy) THEN -! -! write(*,*) 'tendencies are not recalculated !' -! pdu(:,:)=spdu(:,:) -! pdv(:,:)=spdv(:,:) -! pdt(:,:)=spdt(:,:) -! pdq(:,:,:)=spdq(:,:,:) -! pdpsrf(:)=spdpsrf(:) -! write(*,*) pdu(100,10), pdv(100,10), pdt(100,10) -! ELSE -!!!!!!!!!!!!!!!!TEST TEST - - -c ------------------------------------------ -c Initialisations at every physical timestep: -c ------------------------------------------ -c - IF (ngrid.NE.ngridmx) THEN - PRINT*,'STOP in PHYSIQ' - PRINT*,'Probleme de dimensions :' - PRINT*,'ngrid = ',ngrid - PRINT*,'ngridmx = ',ngridmx - STOP - ENDIF - - zday=pday+ptime - - - -c Computing Solar Longitude (Ls) : -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - if (season) then - PRINT*,'day',zday - CALL solarlong(zday,zls) - else - PRINT*,'day_ini',day_ini - call solarlong(float(day_ini),zls) - end if - - -c Initializing various variable -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - call zerophys(ngrid*nlayer, pdv) - call zerophys(ngrid*nlayer, pdu) - call zerophys(ngrid*nlayer, pdt) - call zerophys(ngrid*nlayer*nq, pdq) - call zerophys(ngrid, pdpsrf) - call zerophys(ngrid, zflubid) - call zerophys(ngrid, zdtsurf) - call zerophys(ngrid*nq, dqsurf) - igout=ngrid/2+1 - -c computing geopotentiel at interlayer levels -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -c ponderation des altitudes au niveau des couches en dp/p - - DO l=1,nlayer - DO ig=1,ngrid - zzlay(ig,l)=pphi(ig,l)/g - ENDDO - ENDDO - DO ig=1,ngrid - zzlev(ig,1)=0. - zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km... - ENDDO - DO l=2,nlayer - DO ig=1,ngrid - z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l)) - z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l)) - zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2) - ENDDO - ENDDO - - -! Potential temperature calculation not the same in physiq and dynamic - -c Computing potential temperature -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - DO l=1,nlayer - DO ig=1,ngrid - zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp - zh(ig,l)=pt(ig,l)/zpopsk(ig,l) - ENDDO - ENDDO - -c----------------------------------------------------------------------- -c 1.5 Calculation of mean mass, cp, and R -c --------------------------------------- - - if(photochem.or.callthermos) then - call concentrations(pplay,pt,pdt,pq,pdq,ptimestep) - endif -c----------------------------------------------------------------------- -c 2. Calcul of the radiative tendencies : -c --------------------------------------- - - - IF(callrad) THEN - IF( MOD(icount-1,iradia).EQ.0) THEN - - write (*,*) 'call radiative transfer' - -c Local Solar zenith angle -c ~~~~~~~~~~~~~~~~~~~~~~~~ - - CALL orbite(zls,dist_sol,declin) - - IF(diurnal) THEN - ztim1=SIN(declin) - ztim2=COS(declin)*COS(2.*pi*(zday-.5)) - ztim3=-COS(declin)*SIN(2.*pi*(zday-.5)) - - CALL solang(ngrid,sinlon,coslon,sinlat,coslat, - s ztim1,ztim2,ztim3, mu0,fract) - - ELSE - CALL mucorr(ngrid,declin, lati, mu0, fract,10000.,rad) - ENDIF - - - -c NLTE cooling from CO2 emission -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - IF(callnlte) CALL nltecool(ngrid,nlayer,pplay,pt,zdtnlte) - - -c Find number of layers for LTE radiation calculations - IF(MOD(iphysiq*(icount-1),day_step).EQ.0) THEN - CALL nlthermeq(ngrid,nlayer,pplev,pplay) - ENDIF - - -c Atmospheric dust opacity and aerosol distribution: -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -cc*****WRF - CALL meso_dustopacity(ngrid,nlayer,nq, - $ zday,pplay,pplev,zls,pq, - $ tauref,tau,aerosol) - -cc*****WRF -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - - -c Calling main radiative transfer scheme -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -c Transfer through dust and CO2, except NIR CO2 absorption - -c ---------- -c partie rajoutee par Franck, commentee pour l'instant -c ---------- -c -c if (ngridmx.eq.1) pclc(1)=1. !TEST for 1D simulation -c -c pclc_min=1. -c if (activice.and.naerkind.gt.1) then -c do ig=1,ngrid -c pclc_min=min(pclc_min,pclc(ig)) -c enddo -c endif -c -c IF(activice.and.naerkind.gt.1.and.pclc_min.lt.1) THEN -c Computing radiative tendencies accounting for a cloudy area (0< pclc(ngrid) <1) -c pclc is set in initracer (prescribed for the moment). -c two steps : 1/rad. tend. without clouds (aerosol(*,*,2)=0.) -c ~~~~~~~~~ 2/rad. tend. with clouds -c 3/final tendencies=average of 1/ and 2/ weighted by the -cloud area (pclc) -c -c -c 1/ -c call zerophys(nlayer*ngrid,aerosol(1,1,2)) !remettre -c CALL callradite(icount,ngrid,nlayer,aerosol,albedo, -c $ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, -c $ cldtlw,cldtsw,clsurf_lw,clsurf_sw,cltop_lw,cltop_sw) -c -c 2/ -c CALL h2oiceopacity(ngrid,nlayer,nq,pplay,pplev,pt,pq, -c $ tau,aerosol,zls,co2ice) -c -c CALL callradite(icount,ngrid,nlayer,aerosol,albedo, -c $ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, -c $ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw) -c -c 3/ -c do l=1,nlayer -c do ig=1,ngrid -c zdtlw(ig,l)=(1.-pclc(ig))*cldtlw(ig,l)+pclc(ig)*zdtlw(ig,l) -c zdtsw(ig,l)=(1.-pclc(ig))*cldtsw(ig,l)+pclc(ig)*zdtsw(ig,l) -c enddo -c enddo -c do ig=1,ngrid -c fluxsurf_lw(ig)=(1.-pclc(ig))*clsurf_lw(ig)+ -c $ pclc(ig)*fluxsurf_lw(ig) -c fluxtop_lw(ig)=(1.-pclc(ig))*cltop_lw(ig)+ -c $ pclc(ig)*fluxtop_lw(ig) -c -c fluxsurf_sw(ig,1)=(1.-pclc(ig))*clsurf_sw(ig,1)+ -c $ pclc(ig)*fluxsurf_sw(ig,1) -c fluxsurf_sw(ig,2)=(1.-pclc(ig))*clsurf_sw(ig,2)+ -c $ pclc(ig)*fluxsurf_sw(ig,2) -c -c fluxtop_sw(ig,1)=(1.-pclc(ig))*cltop_sw(ig,1)+ -c $ pclc(ig)*fluxtop_sw(ig,1) -c fluxtop_sw(ig,2)=(1.-pclc(ig))*cltop_sw(ig,2)+ -c $ pclc(ig)*fluxtop_sw(ig,2) -c enddo -c -c ELSE -c -c Atmospheric water ice opacity and particle distribution: -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -c -c if (activice.and.naerkind.gt.1) -c & CALL h2oiceopacity(ngrid,nlayer,nq,pplay,pplev,pt,pq, -c & tau,aerosol,zls,co2ice) -c -c -c ---------- -c fin partie rajoutee par Franck (plus ENDIF ci-dessous) -c ---------- - - -! DO ig=1,ngrid -! DO l=1,nlayer -! IF ( (pplev(ig,l+1) - pplev(ig,l) ) .gt. 0. ) -! . PRINT*,'pb1 ', -! . pplev(ig,l), -! . pplev(ig,l+1), -! . pplev(ig,l+1)-pplev(ig,l), -! . l, -! . ig -! ENDDO -! ENDDO - - CALL callradite(icount,ngrid,nlayer,aerosol,albedo, - $ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, - $ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw) - -! DO ig=1,ngrid -! zdtlw(ig)=zdtlw(1) -! zdtsw(ig)=zdtsw(1) -! fluxsurf_lw(ig)=fluxsurf_lw(1) -! fluxsurf_sw(ig,1)=fluxsurf_sw(1,1) -! fluxsurf_sw(ig,2)=fluxsurf_sw(1,2) -! fluxtop_lw(ig)=fluxtop_lw(1) -! fluxtop_sw(ig,1)=fluxtop_sw(1,1) -! fluxtop_sw(ig,2)=fluxtop_sw(1,2) -! ENDDO - -c ---------- -c ENDIF ! end of condition on the cloudy fraction -c ---------- - - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -ccccc -ccccc PARAM SLOPE : Insolation (direct + scattered) -ccccc - DO ig=1,ngrid - sl_the = theta_sl(ig) - IF (sl_the .ne. 0.) THEN - ztim1=fluxsurf_sw(ig,1)+fluxsurf_sw(ig,2) - DO l=1,2 - sl_lct = ptime*24. + 180.*long(ig)/pi/15. - sl_ra = pi*(1.0-sl_lct/12.) - sl_lat = 180.*lati(ig)/pi - sl_tau = tau(ig,1) - sl_alb = albedo(ig,l) - sl_psi = psi_sl(ig) - sl_fl0 = fluxsurf_sw(ig,l) - sl_di0 = 0. - if (mu0(ig) .gt. 0.) then - sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig))) - sl_di0 = sl_di0*1370./dist_sol/dist_sol - sl_di0 = sl_di0/ztim1 - sl_di0 = fluxsurf_sw(ig,l)*sl_di0 - endif - ! sait-on jamais (a cause des arrondis) - if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0 - !!!!!!!!!!!!!!!!!!!!!!!!!! - CALL param_slope( mu0(ig), declin, sl_ra, sl_lat, - & sl_tau, sl_alb, - & sl_the, sl_psi, sl_di0, sl_fl0, sl_flu) - !!!!!!!!!!!!!!!!!!!!!!!!!! - fluxsurf_sw(ig,l) = sl_flu - !! sl_ls = 180.*zls/pi - !! sl_lct = ptime*24. + 180.*long(ig)/pi/15. - !! sl_lat = 180.*lati(ig)/pi - !! sl_tau = tau(ig,1) - !! sl_alb = albedo(ig,l) - !! sl_the = theta_sl(ig) - !! sl_psi = psi_sl(ig) - !! sl_fl0 = fluxsurf_sw(ig,l) - !! CALL param_slope_full(sl_ls, sl_lct, sl_lat, - !! & sl_tau, sl_alb, - !! & sl_the, sl_psi, sl_fl0, sl_flu) - ENDDO - !!! compute correction on IR flux as well - sky= (1.+cos(pi*theta_sl(ig)/180.))/2. - fluxsurf_lw(ig)= fluxsurf_lw(ig)*sky - ENDIF - ENDDO -ccccc -ccccc PARAM SLOPE -ccccc -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - - -c CO2 near infrared absorption -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - call zerophys(ngrid*nlayer,zdtnirco2) - if (callnirco2) then - call nirco2abs (ngrid,nlayer,pplay,dist_sol, - . mu0,fract,declin, zdtnirco2) - endif - - -c Radiative flux from the sky absorbed by the surface (W.m-2) -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - DO ig=1,ngrid - fluxrad_save(ig)=emis(ig)*fluxsurf_lw(ig) - $ +fluxsurf_sw(ig,1)*(1.-albedo(ig,1)) - $ +fluxsurf_sw(ig,2)*(1.-albedo(ig,2)) - - !print*,'RAD ', fluxrad_save(ig) - !print*,'LW ', emis(ig)*fluxsurf_lw(ig) - !print*,'SW ', fluxsurf_sw(ig,2)*(1.-albedo(ig,2)) - - ENDDO - - - -c Net atmospheric radiative heating rate (K.s-1) -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - IF(callnlte) THEN - CALL blendrad(ngrid, nlayer, pplay, - & zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad) - ELSE - DO l=1,nlayer - DO ig=1,ngrid - dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l) - & +zdtnirco2(ig,l) - ENDDO - ENDDO - ENDIF - - !PRINT*,'zdtsw',zdtsw - !PRINT*,'zdtlw',zdtlw - !PRINT*,'zdtnirco2',zdtnirco2 - !PRINT*,'dtrad',dtrad - - - ENDIF ! mod(icount-1,iradia).eq.0 - -c Transformation of the radiative tendencies: -c ----------------------------------------------- - -c Net radiative surface flux (W.m-2) -c ~~~~~~~~~~~~~~~~~~~~~~~~~~ - DO ig=1,ngrid - zplanck(ig)=tsurf(ig)*tsurf(ig) - zplanck(ig)=emis(ig)* - $ stephan*zplanck(ig)*zplanck(ig) - fluxrad(ig)=fluxrad_save(ig)-zplanck(ig) -cccc -cccc param slope -cccc - sky= (1.+cos(pi*theta_sl(ig)/180.))/2. - fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig) -cccc -cccc -cccc - ENDDO - - - DO l=1,nlayer - DO ig=1,ngrid - pdt(ig,l)=pdt(ig,l)+dtrad(ig,l) - ENDDO - ENDDO - - ENDIF - - - -!c----------------------------------------------------------------------- -!c 3. Gravity wave and subgrid scale topography drag : -!c ------------------------------------------------- -! -! -! IF(calllott)THEN -! -! CALL calldrag_noro(ngrid,nlayer,ptimestep, -! & pplay,pplev,pt,pu,pv,zdtgw,zdugw,zdvgw) -! -! DO l=1,nlayer -! DO ig=1,ngrid -! pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l) -! pdu(ig,l)=pdu(ig,l)+zdugw(ig,l) -! pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l) -! ENDDO -! ENDDO -! ENDIF - -c----------------------------------------------------------------------- -c 4. Vertical diffusion (turbulent mixing): -c ----------------------------------------- -c - - IF(calldifv) THEN - - DO ig=1,ngrid - zflubid(ig)=fluxrad(ig)+fluxgrd(ig) - ENDDO - - CALL zerophys(ngrid*nlayer,zdum1) - CALL zerophys(ngrid*nlayer,zdum2) - do l=1,nlayer - do ig=1,ngrid - zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l) - enddo - enddo - - -c Calling vdif (Martian version WITH CO2 condensation) - CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk, - $ ptimestep,capcal,lwrite, - $ pplay,pplev,zzlay,zzlev,z0, - $ pu,pv,zh,pq,tsurf,emis,qsurf, - $ zdum1,zdum2,zdh,pdq,zflubid, - $ zdudif,zdvdif,zdhdif,zdtsdif,q2, - & zdqdif,zdqsdif) - - DO ig=1,ngrid - !! sensible heat flux in W/m2 - hfx(ig) = zflubid(ig)-capcal(ig)*zdtsdif(ig) - !! u star in similarity theory in m/s - ust(ig) = 0.4 - . * sqrt( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) ) - . / log( 1.E+0 + zzlay(ig,1)/z0 ) - ENDDO - - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!! LES LES - IF (flag_LES) THEN - - write (*,*) '************************************************' - write (*,*) '** LES mode: the difv part is only used to' - write (*,*) '** provide HFX and UST to the dynamics' - write (*,*) '** NB: - dudif, dvdif, dhdif, dqdif are set to 0' - write (*,*) '** - tsurf is updated' - write (*,*) '************************************************' - -! DO ig=1,ngrid -! !! sensible heat flux in W/m2 -! hfx(ig) = zflubid(ig)-capcal(ig)*zdtsdif(ig) -! !! u star in similarity theory in m/s -! ust(ig) = 0.4 -! . * sqrt( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) ) -! . / log( 1.E+0 + zzlay(ig,1)/z0 ) -! - DO l=1,nlayer - zdvdif(ig,l) = 0. - zdudif(ig,l) = 0. - zdhdif(ig,l) = 0. - DO iq=1, nq - zdqdif(ig,l,iq) = 0. - zdqsdif(ig,iq) = 0. !! sortir de la boucle - ENDDO - ENDDO -! -! ENDDO - !write (*,*) 'RAD ',fluxrad(igout) - !write (*,*) 'GRD ',fluxgrd(igout) - !write (*,*) 'dTs/dt ',capcal(igout)*zdtsurf(igout) - !write (*,*) 'HFX ', hfx(igout) - !write (*,*) 'UST ', ust(igout) - ENDIF -!!! LES LES -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - DO l=1,nlayer - DO ig=1,ngrid - - pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l) - pdu(ig,l)=pdu(ig,l)+zdudif(ig,l) - pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l) - - zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only - - ENDDO - ENDDO - - DO ig=1,ngrid - zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig) - ENDDO -!!!!gros caca : sans chaleur sensible -! DO ig=1,ngrid -! zdtsurf(ig)=zdtsurf(ig)+ -! s (fluxrad(ig)+fluxgrd(ig))/capcal(ig) -! ENDDO - - - if(tracer) then - DO iq=1, nq - DO l=1,nlayer - DO ig=1,ngrid - pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq) - ENDDO - ENDDO - ENDDO - DO iq=1, nq - DO ig=1,ngrid - dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq) - ENDDO - ENDDO - end if - - ELSE - DO ig=1,ngrid - zdtsurf(ig)=zdtsurf(ig)+ - s (fluxrad(ig)+fluxgrd(ig))/capcal(ig) - ENDDO -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - IF (flag_LES) THEN - write(*,*) 'LES mode !' - write(*,*) 'Please set calldifv to T in callphys.def' - STOP - ENDIF -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ENDIF - - -c----------------------------------------------------------------------- -c 5. Dry convective adjustment: -c ----------------------------- - - IF(calladj) THEN - - DO l=1,nlayer - DO ig=1,ngrid - zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l) - ENDDO - ENDDO - CALL zerophys(ngrid*nlayer,zduadj) - CALL zerophys(ngrid*nlayer,zdvadj) - CALL zerophys(ngrid*nlayer,zdhadj) - - CALL convadj(ngrid,nlayer,nq,ptimestep, - $ pplay,pplev,zpopsk, - $ pu,pv,zh,pq, - $ pdu,pdv,zdh,pdq, - $ zduadj,zdvadj,zdhadj, - $ zdqadj) - - - DO l=1,nlayer - DO ig=1,ngrid - pdu(ig,l)=pdu(ig,l)+zduadj(ig,l) - pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l) - pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l) - - zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only - ENDDO - ENDDO - - if(tracer) then - DO iq=1, nq - DO l=1,nlayer - DO ig=1,ngrid - pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq) - ENDDO - ENDDO - ENDDO - end if - ENDIF - -c----------------------------------------------------------------------- -c 6. Carbon dioxide condensation-sublimation: -c ------------------------------------------- - - IF(callcond) THEN - CALL newcondens(ngrid,nlayer,nq,ptimestep, - $ capcal,pplay,pplev,tsurf,pt, - $ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, - $ co2ice,albedo,emis, - $ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc, - $ fluxsurf_sw) - - DO l=1,nlayer - DO ig=1,ngrid - pdt(ig,l)=pdt(ig,l)+zdtc(ig,l) - pdv(ig,l)=pdv(ig,l)+zdvc(ig,l) - pdu(ig,l)=pdu(ig,l)+zduc(ig,l) - ENDDO - ENDDO - DO ig=1,ngrid - zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig) - ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep - ENDDO - - IF(tracer) THEN - DO iq=1, nq - DO l=1,nlayer - DO ig=1,ngrid - pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq) - ENDDO - ENDDO - ENDDO - END IF !(tracer) - - ENDIF !(callcond) - -c print*,'condens ok' -c----------------------------------------------------------------------- -c 7. Specific parameterizations for tracers -c: ----------------------------------------- - - if(tracer) then - -c 7a. Water and ice -c --------------- - -c --------------------------------------- -c Water ice condensation in the atmosphere -c ---------------------------------------- - IF (water) THEN - - call watercloud(ngrid,nlayer, ptimestep, - & pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt, - & pq,pdq,zdqcloud,qsurf,zdqscloud,zdtcloud, - & nq,naerkind,tau,icount,zls) - - if (activice) then -c Temperature variation due to latent heat release - DO l=1,nlayer - DO ig=1,ngrid - pdt(ig,l)=pdt(ig,l)+zdtcloud(ig,l) - ENDDO - ENDDO - endif - - IF (iceparty) then - iqmin=nq-1 - ELSE - iqmin=nq - ENDIF - - DO iq=iqmin,nq - DO l=1,nlayer - DO ig=1,ngrid - pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqcloud(ig,l,iq) - ENDDO - ENDDO - DO ig=1,ngrid - dqsurf(ig,iq)= dqsurf(ig,iq) + zdqscloud(ig,iq) - ENDDO - ENDDO - - END IF ! (water) - -c 7b. Chemical species -c ------------------ - -!c -------------- -!c photochemistry : -!c -------------- -! IF(photochem .or. thermochem) then -! call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0, -! . zzlay,zday,pq,pdq,zdqchim,zdqschim,zdqcloud,zdqscloud) -! -!c Photochemistry includes condensation of H2O2 -! -! do iq=nqchem_min,nq -! if (noms(iq).eq."h2o2") then -! DO l=1,nlayer -! DO ig=1,ngrid -! pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqchim(ig,l,iq) -! pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqcloud(ig,l,iq) -! ENDDO -! ENDDO -! else -! DO l=1,nlayer -! DO ig=1,ngrid -! pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqchim(ig,l,iq) -! ENDDO -! ENDDO -! endif -! ENDDO -! do iq=nqchem_min,nq -! if (noms(iq).eq."h2o2") then -! DO ig=1,ngrid -! dqsurf(ig,iq)= dqsurf(ig,iq) + zdqschim(ig,iq) -! dqsurf(ig,iq)= dqsurf(ig,iq) + zdqscloud(ig,iq) -! ENDDO -! else -! DO ig=1,ngrid -! dqsurf(ig,iq)= dqsurf(ig,iq) + zdqschim(ig,iq) -! ENDDO -! endif -! ENDDO -! -! END IF ! (photochem.or.thermochem) -!c print*,'photochem ok' - -c 7c. Aerosol particles -c ------------------- - -c ---------- -c Dust devil : -c ---------- - IF(callddevil) then - call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2, - & zdqdev,zdqsdev) - - if (dustbin.ge.1) then - do iq=1,nq - DO l=1,nlayer - DO ig=1,ngrid - pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq) - ENDDO - ENDDO - enddo - do iq=1,nq - DO ig=1,ngrid - dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq) - ENDDO - enddo - endif ! (test sur dustbin) - - END IF - -c ------------- -c Sedimentation : acts also on water ice -c ------------- - IF (sedimentation) THEN - call zerophys(ngrid*nlayer*nq, zdqsed) - call zerophys(ngrid*nq, zdqssed) - - if(doubleq) then - call callsedim2q(ngrid,nlayer, ptimestep, - & pplev,zzlev, pt, - & pq, pdq, zdqsed, zdqssed,nq) - else - call callsedim(ngrid,nlayer, ptimestep, - & pplev,zzlev, pt, - & pq, pdq, zdqsed, zdqssed,nq) - end if - - - DO iq=1, nq - DO l=1,nlayer - DO ig=1,ngrid - pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq) - ENDDO - ENDDO - ENDDO - DO iq=1, nq - DO ig=1,ngrid - dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq) - ENDDO - ENDDO - END IF ! (sedimentation) - -c print*,'sedim ok' - -c 7d. Updates -c --------- - - DO iq=1, nq - DO ig=1,ngrid - -c --------------------------------- -c Updating tracer budget on surface -c --------------------------------- - qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq) - - ENDDO ! (ig) - ENDDO ! (iq) - - END IF ! (tracer) - -c print*,'tracers ok' - - -!c----------------------------------------------------------------------- -!c 8.5 THERMOSPHERE CALCULATION -!c----------------------------------------------------------------------- -! -! if (callthermos) then -! call thermosphere(pplev,pplay,dist_sol, -! $ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay, -! & pt,pq,pu,pv,pdt,pdq, -! $ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff) -!c do iq=nqchem_min,nq -!c write(*,*) 'thermo iq,pq',iq,pq(690,1,iq) -!c enddo -! -! DO l=1,nlayer -! DO ig=1,ngrid -! dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l) -! pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l) -! & +zdteuv(ig,l) -! pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l) -! pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l) -! DO iq=1, nq -! pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq) -! ENDDO -! ENDDO -! ENDDO -! -! -! endif - -c----------------------------------------------------------------------- -c 8. Surface and sub-surface soil temperature -c----------------------------------------------------------------------- -c -c -c Surface temperature incrementation : -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - DO ig=1,ngrid - tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig) - ENDDO - - -c Prescription piege froid au pole sud (Except at high obliquity !!) -c temperature en surface = temperature equilibre de phases du CO2 - - IF (tracer.AND.water.AND.ngridmx.NE.1) THEN - !if (caps.and.(obliquit.lt.27.)) then - ! tsurf(ngrid)=1/(1/136.27-r/5.9e+5*alog(0.0095*ps(ngrid))) - !endif -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!!!! note WRF MESOSCALE AYMERIC -- mot cle "caps" -!!!!! watercaptag n'est plus utilise que dans vdifc -!!!!! ... pour que la sublimation ne soit pas stoppee -!!!!! ... dans la calotte permanente nord si qsurf=0 -!!!!! on desire garder cet effet regle par caps=T -!!!!! on a donc commente "if (caps.and.(obliquit.lt.27.))" ci-dessus -!!!!! --- remplacer ces lignes par qqch de plus approprie -!!!!! si on s attaque a la calotte polaire sud -!!!!! pas d'autre occurrence majeure du mot-cle "caps" -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -c ------------------------------------------------------------- -c Change of surface albedo (set to 0.4) in case of ground frost -c everywhere except on the north permanent cap and in regions -c covered by dry ice. -c ALWAYS PLACE these lines after newcondens !!! -c ------------------------------------------------------------- - do ig=1,ngrid - -c -------------- Version temporaire fit TES 2008 ------------ - if (co2ice(ig).eq.0.and.qsurf(ig,nqmx).gt.0.005) then - albedo(ig,1)=0.45 - albedo(ig,2)=0.45 - endif - -c if (co2ice(ig).eq.0.and.qsurf(ig,nqmx).gt.0.005) then -c if (.not.watercaptag(ig)) then -c albedo(ig,1)=0.4 -c albedo(ig,2)=0.4 -c endif -c endif -c -------------- version Francois --------------- -c if (co2ice(ig).eq.0.and. -c & ((qsurf(ig,nqmx).gt.0.005).or.(watercaptag(ig)))) then -c albedo(ig,1)=max(0.4,albedodat(ig)) -c albedo(ig,2)=albedo(ig,1) -c endif -c --------------------------------------------- - enddo ! of do ig=1,ngrid - ENDIF ! of IF (tracer.AND.water.AND.ngridmx.NE.1) - - -c print*,'tracer, water and 3D ok' - -c -c soil temperatures and subsurface heat flux: -c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - IF (callsoil) THEN - CALL soil(ngrid,nsoilmx,.false.,inertiedat, - s ptimestep,tsurf,tsoil,capcal,fluxgrd) - ENDIF - -c print*,'soil ok' - - - - - -c----------------------------------------------------------------------- -c 9. Writing output files -c ------------------------ - -c ------------------------------- -c Dynamical fields incrementation -c ------------------------------- -c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics) - - DO l=1,nlayer - DO ig=1,ngrid - zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep - zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep - zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep - ENDDO - ENDDO - DO iq=1, nq - DO l=1,nlayer - DO ig=1,ngrid - zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep - ENDDO - ENDDO - ENDDO - DO ig=1,ngrid - ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep !already in 7 - ENDDO - DO l=1,nlayer - DO ig=1,ngrid - zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig) - zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig) - ENDDO - ENDDO - ! Density calculation - DO l=1,nlayer - DO ig=1,ngrid - rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l)) - ENDDO - ENDDO - -c Sum of fluxes in solar spectral bands (for output only) - DO ig=1,ngrid - fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2) - fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2) - ENDDO -c ******* TEMPORAIRE ************************************************ - ztim1 = 999 - DO l=1,nlayer - DO ig=1,ngrid - if (pt(ig,l).lt.ztim1) then - ztim1 = pt(ig,l) - igmin = ig - lmin = l - end if - ENDDO - ENDDO - if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then - write(*,*) 'stability WARNING :' - write(*,*) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin), - & 'ig l =', igmin, lmin - end if -c ******************************************************************* - -c -------------------- -c Output on the screen -c -------------------- - - - IF (lwrite) THEN - PRINT*,'Global diagnostics for the physics' - PRINT*,'Variables and their increments x and dx/dt * dt' - WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps' - WRITE(*,'(2f10.5,2f15.5)') - s tsurf(igout),zdtsurf(igout)*ptimestep, - s pplev(igout,1),pdpsrf(igout)*ptimestep - WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT' - WRITE(*,'(i4,6f10.5)') (l, - s pu(igout,l),pdu(igout,l)*ptimestep, - s pv(igout,l),pdv(igout,l)*ptimestep, - s pt(igout,l),pdt(igout,l)*ptimestep, - s l=1,nlayer) - ENDIF - -cc****WRF -cc IF (ngrid.NE.1) THEN -c PRINT*,'check - values after update at ngrid/2+1' -c WRITE(*,'(a4,a6,2a10)') 'l','u','v','T' -c WRITE(*,'(i4,3f10.5)') (l, -c s zu(igout,l), -c s zv(igout,l), -c s zt(igout,l), -c s l=1,nlayer) -c -cc****WRF - - print*,'Ls =',zls*180./pi, - & ' tauref(700 Pa,lat=0) =',tauref(ngrid/2) -c & ' tau(Viking1) =',tau(ig_vl1,1) - - -c ------------------------------------------------------------------- -c Writing NetCDF file "RESTARTFI" at the end of the run -c ------------------------------------------------------------------- -c Remarque : On stocke restarfi -c juste avant que la dynamique ne le soit dans restart. -c entre maintenant et l'ecriture de restart, -c on aura itau = itau +1 et remise a jour de time. -c (lastcall = .true. lorsque itau+1 = itaufin) -c Donc on stocke avec time = time + dtvr - - IF(lastcall) THEN - ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec) - write(*,*)'pour physdem ztime_fin =',ztime_fin - call physdem1("restartfi.nc",long,lati,nsoilmx,nq, - . ptimestep,pday, - . ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf, - . area,albedodat,inertiedat,zmea,zstd,zsig, - . zgam,zthe) - ENDIF - - -ccc**************** WRF OUTPUT ************************** -ccc**************** WRF OUTPUT ************************** -ccc**************** WRF OUTPUT ************************** - do ig=1,ngrid - wtsurf(ig) = tsurf(ig) !! surface temperature - wco2ice(ig) = co2ice(ig) !! co2 ice - - !!! TEMP TEMP TEMP TEMP TEMP TEMP TEMP - !!! specific to WRF WRF WRF - !!! just to output water ice on surface - !!! [it might not be water ice on surface but OK] - !!! uncomment the Registry entry - qsurflast(ig) = qsurf(ig,nqmx) - - enddo -c -c THIS INCLUDE IS AUTOMATICALLY GENERATED FROM REGISTRY -c -#include "fill_save.inc" -c -ccc**************** WRF OUTPUT ************************** -ccc**************** WRF OUTPUT ************************** -ccc**************** WRF OUTPUT ************************** - - -cc----------------------------------- -cc you can still use meso_WRITEDIAGFI (e.g. for debugging purpose), -cc though this is not the default strategy now -cc----------------------------------- -cc please use cudt in namelist.input to set frequency of outputs -cc----------------------------------- -cc BEWARE: if at least one call to meso_WRITEDIAGFI is performed, -cc cudt cannot be 0 - otherwise you'll get a "Floating exception" -cc----------------------------------- -! call meso_WRITEDIAGFI(ngrid,"tauref", -! . "tauref","W.m-2",2, -! . tauref) -! call meso_WRITEDIAGFI(ngrid,"dtrad", -! . "dtrad","W.m-2",2, -! . dtrad) -c call meso_WRITEDIAGFI(ngrid,"tsurf", -c . "tsurf","K",2, -c . tsurf) -c -! call meso_WRITEDIAGFI(ngrid,"zt", -! . "zt","W.m-2",3, -! . zt) -! call meso_WRITEDIAGFI(ngrid,"zdtlw", -! . "zdtlw","W.m-2",2, -! . zdtlw) -! call meso_WRITEDIAGFI(ngrid,"zdtsw", -! . "zdtsw","W.m-2",2, -! . zdtsw) - - - icount=icount+1 - -!!!!!!!!!!!!!TEST TEST -! ENDIF -!!!!!!!!!!!!!TEST TEST - - write(*,*) 'now, back to the dynamical core...' - RETURN - END Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/meso_testphys1d.F =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/meso_testphys1d.F (revision 234) +++ (revision ) @@ -1,508 +1,0 @@ - - PROGRAM meso_testphys1d - IMPLICIT NONE - -c======================================================================= -c subject: -c -------- -c PROGRAM useful to run physical part of the martian GCM in a 1D column -c -c Can be compiled with a command like (e.g. for 25 layers) -c "makegcm -p mars -d 25 testphys1d" -c It requires the files "testphys1d.def" "callphys.def" -c and a file describing the sigma layers (e.g. "z2sig.def") -c -c author: Frederic Hourdin, R.Fournier,F.Forget -c ------- -c -c update: 12/06/2003 including chemistry (S. Lebonnois) -c and water ice (F. Montmessin) -c -c======================================================================= - -#include "dimensions.h" -#include "dimphys.h" -#include "dimradmars.h" -#include "comgeomfi.h" -#include "surfdat.h" -#include "comdiurn.h" -#include "callkeys.h" -#include "comcstfi.h" -#include "planete.h" -#include "comsaison.h" -#include "yomaer.h" -#include "aerdust.h" -#include "control.h" -#include "comvert.h" -#include "netcdf.inc" -#include "comg1d.h" -#include "watercap.h" -#include "fisice.h" -#include "logic.h" - -c -------------------------------------------------------------- -c Declarations -c -------------------------------------------------------------- -c - INTEGER unit ! unite de lecture de "testphys1d.def" - INTEGER unitstart ! unite d'ecriture de "startfi" - INTEGER nlayer,nlevel,nsoil,ndt - INTEGER ilayer,ilevel,isoil,idt,iq - LOGICAl firstcall,lastcall -c - INTEGER day0 ! date initial (sol ; =0 a Ls=0) - REAL day ! date durant le run - REAL time ! time (0 - area(1)=1.E+0 - -c le geopotentiel au sol est inutile en 1D car tout est controle -c par la pression de surface ---> - phisfi(1)=0.E+0 - -c "inifis" reproduit un certain nombre d'initialisations deja faites -c + lecture des clefs de callphys.def -c + calcul de la frequence d'appel au rayonnement -c + calcul des sinus et cosinus des longitude latitude - -!Mars possible matter with dtphys in input and include!!! - CALL meso_inifis(1,llm,day0,daysec,dtphys, - . lati,long,area,rad,g,r,cpp) -c Initialisation pour prendre en compte les vents en 1-D -c ------------------------------------------------------ - ptif=2.E+0*omeg*sinlat(1) - -c vent geostrophique - PRINT *,'composante vers l est du vent geostrophique (U) ?' - READ(unit,*) gru - PRINT *,'composante vers le nord du vent geostrophique (V) ?' - READ(unit,*) grv - -c Initialisation des vents au premier pas de temps - DO ilayer=1,nlayer - u(ilayer)=gru - v(ilayer)=grv - ENDDO - -c energie cinetique turbulente - DO ilevel=1,nlevel - q2(ilevel)=0.E+0 - ENDDO - -c glace de CO2 au sol -c ------------------- - co2ice=0.E+0 - PRINT *,'co2ice (kg.m-2)' - READ(unit,*) co2ice - -c -c emissivite -c ---------- - emis=emissiv - IF (co2ice.eq.1.E+0) THEN - emis=emisice(1) - IF(lati(1).LT.0) emis=emisice(2) - ENDIF - - - -c calcul des pressions et altitudes en utilisant les niveaux sigma -c ---------------------------------------------------------------- - -c Vertical Coordinates -c """""""""""""""""""" - PRINT *,'Hybrid coordinates ?' - READ(unit,*) hybrid - PRINT *,'Hybrid =', hybrid - - CALL disvert - - DO ilevel=1,nlevel - plev(ilevel)=ap(ilevel)+psurf*bp(ilevel) - ENDDO - - DO ilayer=1,nlayer - play(ilayer)=aps(ilayer)+psurf*bps(ilayer) - ENDDO - - DO ilayer=1,nlayer - zlay(ilayer)=-200.E+0 *r*log(play(ilayer)/plev(1)) - & /g - ENDDO - - -c profil de temperature au premier appel -c -------------------------------------- - pks=psurf**rcp - -c altitude en km dans profile: on divise zlay par 1000 - tmp1(0)=0.E+0 - DO ilayer=1,nlayer - tmp1(ilayer)=zlay(ilayer)/1000.E+0 - ENDDO - call profile(unit,nlayer+1,tmp1,tmp2) - - tsurf=tmp2(0) - DO ilayer=1,nlayer - temp(ilayer)=tmp2(ilayer) - ENDDO - - - -c temperature du sous-sol -c ~~~~~~~~~~~~~~~~~~~~~~~ - DO isoil=1,nsoil - tsoil(isoil)=tsurf - ENDDO - -c Initialisation des traceurs -c --------------------------- - - DO iq=1,nqmx - DO ilayer=1,nlayer - q(ilayer,iq) = 0. - ENDDO - ENDDO - - if (photochem.or.callthermos) then - write(*,*) 'Especes chimiques initialisees' - ! thermo=0: initialisation pour toutes les couches - thermo=0 - call inichim_newstart(q,psurf,sig,nqmx,lati,long,area, - $ thermo) - endif - watercaptag(ngridmx)=.false. - - DO iq=1,nqmx-1 - qsurf(iq) = 0. - ENDDO - - -c Initialisation pour les sorties GRADS dans "g1d.dat" et "g1d.ctl" -c ---------------------------------------------------------------- -c (effectuee avec "writeg1d" a partir de "physiq.F" -c ou tout autre subroutine physique, y compris celle ci). - - g1d_nlayer=nlayer - g1d_nomfich='g1d.dat' - g1d_unitfich=40 - g1d_nomctl='g1d.ctl' - g1d_unitctl=41 - g1d_premier=.true. - g2d_premier=.true. - -c Ecriture de "startfi" -c -------------------- -c (Ce fichier sera aussitot relu au premier -c appel de "physiq", mais il est necessaire pour passer -c les variables purement physiques a "physiq"... - - call physdem1("startfi.nc",long,lati,nsoilmx,nqmx, - . dtphys,float(day0), - . time,tsurf,tsoil,co2ice,emis,q2,qsurf, - . area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe) -c======================================================================= -c BOUCLE TEMPORELLE DU MODELE 1D -c======================================================================= -c - firstcall=.true. - lastcall=.false. - - DO idt=1,ndt -c IF (idt.eq.ndt) lastcall=.true. - IF (idt.eq.ndt-day_step-1) then !test - lastcall=.true. - call solarlong(day*1.0,zls) - write(103,*) 'Ls=',zls*180./pi - write(103,*) 'Lat=', lati(1)*180./pi - write(103,*) 'Tau=', tauvis/700*psurf - write(103,*) 'RunEnd - Atmos. Temp. File' - write(103,*) 'RunEnd - Atmos. Temp. File' - write(104,*) 'Ls=',zls*180./pi - write(104,*) 'Lat=', lati(1) - write(104,*) 'Tau=', tauvis/700*psurf - write(104,*) 'RunEnd - Atmos. Temp. File' - ENDIF - -c calcul du geopotentiel -c ~~~~~~~~~~~~~~~~~~~~~ - DO ilayer=1,nlayer - s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp - h(ilayer)=cpp*temp(ilayer)/(pks*s(ilayer)) - ENDDO - phi(1)=pks*h(1)*(1.E+0-s(1)) - DO ilayer=2,nlayer - phi(ilayer)=phi(ilayer-1)+ - & pks*(h(ilayer-1)+h(ilayer))*.5E+0 - & *(s(ilayer-1)-s(ilayer)) - - ENDDO - -c appel de la physique -c -------------------- - - CALL meso_physiq (1,llm,nqmx, - , firstcall,lastcall, - , day,time,dtphys, - , plev,play,phi, - , u, v,temp, q, - , w, -C - sorties - s du, dv, dtemp, dq,dpsurf,tracerdyn) - -c evolution du vent : modele 1D -c ----------------------------- - -c la physique calcule les derivees temporelles de u et v. -c on y rajoute betement un effet Coriolis. -c -c DO ilayer=1,nlayer -c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv) -c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru) -c ENDDO - -c Pour certain test : pas de coriolis a l'equateur -c if(lati(1).eq.0.) then - DO ilayer=1,nlayer - du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4 - dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4 - ENDDO -c end if -c -c -c Calcul du temps au pas de temps suivant -c --------------------------------------- - firstcall=.false. - time=time+dtphys/daysec - IF (time.gt.1.E+0) then - time=time-1.E+0 - day=day+1 - ENDIF - -c calcul des vitesses et temperature au pas de temps suivant -c ---------------------------------------------------------- - - DO ilayer=1,nlayer - u(ilayer)=u(ilayer)+dtphys*du(ilayer) - v(ilayer)=v(ilayer)+dtphys*dv(ilayer) - temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer) - ENDDO - -c calcul des pressions au pas de temps suivant -c ---------------------------------------------------------- - - psurf=psurf+dtphys*dpsurf ! evolution de la pression de surface - DO ilevel=1,nlevel - plev(ilevel)=ap(ilevel)+psurf*bp(ilevel) - ENDDO - DO ilayer=1,nlayer - play(ilayer)=aps(ilayer)+psurf*bps(ilayer) - ENDDO - - ENDDO ! fin de la boucle temporelle - -c ======================================================== -c GESTION DES SORTIE -c ======================================================== - -c fermeture pour conclure les sorties format grads dans "g1d.dat" -c et "g1d.ctl" (effectuee avec "writeg1d" a partir de "physiq.F" -c ou tout autre subroutine physique, y compris celle ci). - -c CALL endg1d(1,nlayer,zphi/(g*1000.),ndt) - CALL endg1d(1,nlayer,zlay/1000.,ndt) - -c ======================================================== - END - -c*********************************************************************** -c*********************************************************************** -c Subroutines Bidons utilise seulement en 3D, mais -c necessaire a la compilation de testphys1d en 1-D - - subroutine gr_fi_dyn - RETURN - END - - subroutine iniwrite - RETURN - END - -c*********************************************************************** -c*********************************************************************** - -#include "../dyn3d/disvert.F" Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/physiq.F =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/physiq.F (revision 235) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/physiq.F (revision 235) @@ -0,0 +1,1569 @@ + SUBROUTINE physiq(ngrid,nlayer,nq, + $ firstcall,lastcall, + $ pday,ptime,ptimestep, + $ pplev,pplay,pphi, + $ pu,pv,pt,pq, + $ pw, + $ pdu,pdv,pdt,pdq,pdpsrf,tracerdyn, + $ wtsurf,wtsoil,wemis,wq2,wqsurf,wco2ice, + $ wday_ini, + $ output_tab2d, output_tab3d, + $ flag_LES) + + + IMPLICIT NONE +c======================================================================= +c +c CAREFUL: THIS IS A VERSION TO BE USED WITH WRF !!! +c +c ... CHECK THE ****WRF lines +c +c======================================================================= +c +c subject: +c -------- +c +c Organisation of the physical parametrisations of the LMD +c martian atmospheric general circulation model. +c +c The GCM can be run without or with tracer transport +c depending on the value of Logical "tracer" in file "callphys.def" +c Tracers may be water vapor, ice OR chemical species OR dust particles +c +c SEE comments in initracer.F about numbering of tracer species... +c +c It includes: +c +c 1. Initialisation: +c 1.5 Calculation of mean mass and cp, R and thermal conduction coeff. +c 2. Calcul of the radiative tendencies : radiative transfer +c (longwave and shortwave) for CO2 and dust. +c 3. Gravity wave and subgrid scale topography drag : +c 4. Vertical diffusion (turbulent mixing): +c 5. convective adjustment +c 6. TRACERS : +c 6a. water and water ice +c 6b. call for photochemistry when tracers are chemical species +c 6c. other scheme for tracer (dust) transport (lifting, sedimentation) +c 6d. updates (CO2 pressure variations, surface budget) +c 7. condensation and sublimation of carbon dioxide. +c 8. Surface and sub-surface temperature calculations +c 9. Writing output files : +c - "startfi", "histfi" (if it's time) +c - Saving statistics (if "callstats = .true.") +c - Dumping eof (if "calleofdump = .true.") +c - Output any needed variables in "diagfi" +c 10. Diagnostic: mass conservation of tracers +c +c author: +c ------- +c Frederic Hourdin 15/10/93 +c Francois Forget 1994 +c Christophe Hourdin 02/1997 +c Subroutine completly rewritten by F.Forget (01/2000) +c Introduction of the photochemical module: S. Lebonnois (11/2002) +c Introduction of the thermosphere module: M. Angelats i Coll (2002) +c Water ice clouds: Franck Montmessin (update 06/2003) +c WRF version: Aymeric Spiga (01-03/2007) +c +c +c +c arguments: +c ---------- +c +c input: +c ------ +c ecri period (in dynamical timestep) to write output +c ngrid Size of the horizontal grid. +c All internal loops are performed on that grid. +c nlayer Number of vertical layers. +c nq Number of advected fields +c firstcall True at the first call +c lastcall True at the last call +c pday Number of days counted from the North. Spring +c equinoxe. +c ptime Universal time (0 part of the job of phyetat0 is done in inifis +c > remaining initializations are passed here from the WRF variables +c > beware, some operations were done by phyetat0 (ex: tracers) +c > if any problems, look in phyetat0 +c + tsurf(:)=wtsurf(:) + PRINT*,'check: tsurf ',tsurf(1),tsurf(ngridmx) + tsoil(:,:)=wtsoil(:,:) + PRINT*,'check: tsoil ',tsoil(1,1),tsoil(ngridmx,nsoilmx) + emis(:)=wemis(:) + PRINT*,'check: emis ',emis(1),emis(ngridmx) + q2(:,:)=wq2(:,:) + PRINT*,'check: q2 ',q2(1,1),q2(ngridmx,nlayermx+1) + qsurf(:,:)=wqsurf(:,:) + PRINT*,'check: qsurf ',qsurf(1,1),qsurf(ngridmx,nqmx) + co2ice(:)=wco2ice(:) + PRINT*,'check: co2 ',co2ice(1),co2ice(ngridmx) + day_ini=wday_ini + +c artificially filling dyn3d/control.h is also required +c > iphysiq is put in WRF to be set easily (cf ptimestep) +c > day_step is simply deduced: +c + day_step=daysec/ptimestep + PRINT*,'Call to LMD physics:',day_step,' per Martian day' +c + iphysiq=ptimestep +c + ecritphy=8.e18 !! a dummy low frequency + PRINT*,'Write LMD physics each:',ecritphy,' seconds' + !!PRINT*,ecri_phys + !!PRINT*,float(ecri_phys) ... + !!renvoient tous deux des nombres absurdes + !!pourtant callkeys.h est inclus ... + !! + !!donc ecritphys est passe en argument ... + PRINT*,'Write LMD physics each:',ecritphy,' seconds' +c +c ****WRF +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + + if (pday.ne.day_ini) then + write(*,*) "***ERROR Pb de synchronisation entre phys et dyn" + write(*,*) "jour dynamique: ",pday + write(*,*) "jour physique: ",day_ini + stop + endif + + write (*,*) 'In physic day_ini =', day_ini + + + +c Initialize albedo and orbital calculation +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + CALL surfini(ngrid,co2ice,qsurf,albedo) + CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit) + +c initialisation soil +c ~~~~~~~~~~~~~~~~~~~ + IF (callsoil) THEN + CALL soil(ngrid,nsoilmx,firstcall,inertiedat, + s ptimestep,tsurf,tsoil,capcal,fluxgrd) + ELSE + PRINT*,'WARNING! Thermal conduction in the soil turned off' + DO ig=1,ngrid + capcal(ig)=1.e5 + fluxgrd(ig)=0. + ENDDO + ENDIF + icount=1 + + +c initialisation pour les traceurs +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + tracerdyn=tracer + IF (tracer) THEN + CALL initracer(qsurf,co2ice) + ENDIF ! end tracer + + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c ****WRF +c +c nosense in mesoscale modeling +c +cc Determining gridpoint near Viking Lander 1 (used for diagnostic only) +cc ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +c if(ngrid.ne.1) then +c latvl1= 22.27 +c lonvl1= -47.94 +c ig_vl1= 1+ int( (1.5-(latvl1-90.)*jjm/180.) -2 )*iim + +c & int(1.5+(lonvl1+180)*iim/360.) +c write(*,*) 'Viking Lander 1 GCM point: lat,lon', +c & lati(ig_vl1)*180/pi, long(ig_vl1)*180/pi +c end if +c +c ****WRF +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + +c Initializing thermospheric parameters +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + if (callthermos) call param_read + +c Initializing R and Cp as constant + + if (.not.callthermos .and. .not.photochem) then + do l=1,nlayermx + do ig=1,ngridmx + rnew(ig,l)=r + cpnew(ig,l)=cpp + mmean(ig,l)=mugaz + enddo + enddo + endif + + ENDIF ! (end of "if firstcall") + +!!!!!!!!!!!!!!!!TEST TEST +! IF (nocalculphy) THEN +! +! write(*,*) 'tendencies are not recalculated !' +! pdu(:,:)=spdu(:,:) +! pdv(:,:)=spdv(:,:) +! pdt(:,:)=spdt(:,:) +! pdq(:,:,:)=spdq(:,:,:) +! pdpsrf(:)=spdpsrf(:) +! write(*,*) pdu(100,10), pdv(100,10), pdt(100,10) +! ELSE +!!!!!!!!!!!!!!!!TEST TEST + + +c ------------------------------------------ +c Initialisations at every physical timestep: +c ------------------------------------------ +c + IF (ngrid.NE.ngridmx) THEN + PRINT*,'STOP in PHYSIQ' + PRINT*,'Probleme de dimensions :' + PRINT*,'ngrid = ',ngrid + PRINT*,'ngridmx = ',ngridmx + STOP + ENDIF + + zday=pday+ptime + + + +c Computing Solar Longitude (Ls) : +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + if (season) then + PRINT*,'day',zday + CALL solarlong(zday,zls) + else + PRINT*,'day_ini',day_ini + call solarlong(float(day_ini),zls) + end if + + +c Initializing various variable +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + call zerophys(ngrid*nlayer, pdv) + call zerophys(ngrid*nlayer, pdu) + call zerophys(ngrid*nlayer, pdt) + call zerophys(ngrid*nlayer*nq, pdq) + call zerophys(ngrid, pdpsrf) + call zerophys(ngrid, zflubid) + call zerophys(ngrid, zdtsurf) + call zerophys(ngrid*nq, dqsurf) + igout=ngrid/2+1 + +c computing geopotentiel at interlayer levels +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +c ponderation des altitudes au niveau des couches en dp/p + + DO l=1,nlayer + DO ig=1,ngrid + zzlay(ig,l)=pphi(ig,l)/g + ENDDO + ENDDO + DO ig=1,ngrid + zzlev(ig,1)=0. + zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km... + ENDDO + DO l=2,nlayer + DO ig=1,ngrid + z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l)) + z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l)) + zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2) + ENDDO + ENDDO + + +! Potential temperature calculation not the same in physiq and dynamic + +c Computing potential temperature +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + DO l=1,nlayer + DO ig=1,ngrid + zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp + zh(ig,l)=pt(ig,l)/zpopsk(ig,l) + ENDDO + ENDDO + +c----------------------------------------------------------------------- +c 1.5 Calculation of mean mass, cp, and R +c --------------------------------------- + + if(photochem.or.callthermos) then + call concentrations(pplay,pt,pdt,pq,pdq,ptimestep) + endif +c----------------------------------------------------------------------- +c 2. Calcul of the radiative tendencies : +c --------------------------------------- + + + IF(callrad) THEN + IF( MOD(icount-1,iradia).EQ.0) THEN + + write (*,*) 'call radiative transfer' + +c Local Solar zenith angle +c ~~~~~~~~~~~~~~~~~~~~~~~~ + + CALL orbite(zls,dist_sol,declin) + + IF(diurnal) THEN + ztim1=SIN(declin) + ztim2=COS(declin)*COS(2.*pi*(zday-.5)) + ztim3=-COS(declin)*SIN(2.*pi*(zday-.5)) + + CALL solang(ngrid,sinlon,coslon,sinlat,coslat, + s ztim1,ztim2,ztim3, mu0,fract) + + ELSE + CALL mucorr(ngrid,declin, lati, mu0, fract,10000.,rad) + ENDIF + + + +c NLTE cooling from CO2 emission +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + IF(callnlte) CALL nltecool(ngrid,nlayer,pplay,pt,zdtnlte) + + +c Find number of layers for LTE radiation calculations + IF(MOD(iphysiq*(icount-1),day_step).EQ.0) THEN + CALL nlthermeq(ngrid,nlayer,pplev,pplay) + ENDIF + + +c Atmospheric dust opacity and aerosol distribution: +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +cc*****WRF + CALL meso_dustopacity(ngrid,nlayer,nq, + $ zday,pplay,pplev,zls,pq, + $ tauref,tau,aerosol) + +cc*****WRF +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + +c Calling main radiative transfer scheme +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +c Transfer through dust and CO2, except NIR CO2 absorption + +c ---------- +c partie rajoutee par Franck, commentee pour l'instant +c ---------- +c +c if (ngridmx.eq.1) pclc(1)=1. !TEST for 1D simulation +c +c pclc_min=1. +c if (activice.and.naerkind.gt.1) then +c do ig=1,ngrid +c pclc_min=min(pclc_min,pclc(ig)) +c enddo +c endif +c +c IF(activice.and.naerkind.gt.1.and.pclc_min.lt.1) THEN +c Computing radiative tendencies accounting for a cloudy area (0< pclc(ngrid) <1) +c pclc is set in initracer (prescribed for the moment). +c two steps : 1/rad. tend. without clouds (aerosol(*,*,2)=0.) +c ~~~~~~~~~ 2/rad. tend. with clouds +c 3/final tendencies=average of 1/ and 2/ weighted by the +cloud area (pclc) +c +c +c 1/ +c call zerophys(nlayer*ngrid,aerosol(1,1,2)) !remettre +c CALL callradite(icount,ngrid,nlayer,aerosol,albedo, +c $ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, +c $ cldtlw,cldtsw,clsurf_lw,clsurf_sw,cltop_lw,cltop_sw) +c +c 2/ +c CALL h2oiceopacity(ngrid,nlayer,nq,pplay,pplev,pt,pq, +c $ tau,aerosol,zls,co2ice) +c +c CALL callradite(icount,ngrid,nlayer,aerosol,albedo, +c $ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, +c $ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw) +c +c 3/ +c do l=1,nlayer +c do ig=1,ngrid +c zdtlw(ig,l)=(1.-pclc(ig))*cldtlw(ig,l)+pclc(ig)*zdtlw(ig,l) +c zdtsw(ig,l)=(1.-pclc(ig))*cldtsw(ig,l)+pclc(ig)*zdtsw(ig,l) +c enddo +c enddo +c do ig=1,ngrid +c fluxsurf_lw(ig)=(1.-pclc(ig))*clsurf_lw(ig)+ +c $ pclc(ig)*fluxsurf_lw(ig) +c fluxtop_lw(ig)=(1.-pclc(ig))*cltop_lw(ig)+ +c $ pclc(ig)*fluxtop_lw(ig) +c +c fluxsurf_sw(ig,1)=(1.-pclc(ig))*clsurf_sw(ig,1)+ +c $ pclc(ig)*fluxsurf_sw(ig,1) +c fluxsurf_sw(ig,2)=(1.-pclc(ig))*clsurf_sw(ig,2)+ +c $ pclc(ig)*fluxsurf_sw(ig,2) +c +c fluxtop_sw(ig,1)=(1.-pclc(ig))*cltop_sw(ig,1)+ +c $ pclc(ig)*fluxtop_sw(ig,1) +c fluxtop_sw(ig,2)=(1.-pclc(ig))*cltop_sw(ig,2)+ +c $ pclc(ig)*fluxtop_sw(ig,2) +c enddo +c +c ELSE +c +c Atmospheric water ice opacity and particle distribution: +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +c +c if (activice.and.naerkind.gt.1) +c & CALL h2oiceopacity(ngrid,nlayer,nq,pplay,pplev,pt,pq, +c & tau,aerosol,zls,co2ice) +c +c +c ---------- +c fin partie rajoutee par Franck (plus ENDIF ci-dessous) +c ---------- + + +! DO ig=1,ngrid +! DO l=1,nlayer +! IF ( (pplev(ig,l+1) - pplev(ig,l) ) .gt. 0. ) +! . PRINT*,'pb1 ', +! . pplev(ig,l), +! . pplev(ig,l+1), +! . pplev(ig,l+1)-pplev(ig,l), +! . l, +! . ig +! ENDDO +! ENDDO + + CALL callradite(icount,ngrid,nlayer,aerosol,albedo, + $ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, + $ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw) + +! DO ig=1,ngrid +! zdtlw(ig)=zdtlw(1) +! zdtsw(ig)=zdtsw(1) +! fluxsurf_lw(ig)=fluxsurf_lw(1) +! fluxsurf_sw(ig,1)=fluxsurf_sw(1,1) +! fluxsurf_sw(ig,2)=fluxsurf_sw(1,2) +! fluxtop_lw(ig)=fluxtop_lw(1) +! fluxtop_sw(ig,1)=fluxtop_sw(1,1) +! fluxtop_sw(ig,2)=fluxtop_sw(1,2) +! ENDDO + +c ---------- +c ENDIF ! end of condition on the cloudy fraction +c ---------- + + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +ccccc +ccccc PARAM SLOPE : Insolation (direct + scattered) +ccccc + DO ig=1,ngrid + sl_the = theta_sl(ig) + IF (sl_the .ne. 0.) THEN + ztim1=fluxsurf_sw(ig,1)+fluxsurf_sw(ig,2) + DO l=1,2 + sl_lct = ptime*24. + 180.*long(ig)/pi/15. + sl_ra = pi*(1.0-sl_lct/12.) + sl_lat = 180.*lati(ig)/pi + sl_tau = tau(ig,1) + sl_alb = albedo(ig,l) + sl_psi = psi_sl(ig) + sl_fl0 = fluxsurf_sw(ig,l) + sl_di0 = 0. + if (mu0(ig) .gt. 0.) then + sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig))) + sl_di0 = sl_di0*1370./dist_sol/dist_sol + sl_di0 = sl_di0/ztim1 + sl_di0 = fluxsurf_sw(ig,l)*sl_di0 + endif + ! sait-on jamais (a cause des arrondis) + if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0 + !!!!!!!!!!!!!!!!!!!!!!!!!! + CALL param_slope( mu0(ig), declin, sl_ra, sl_lat, + & sl_tau, sl_alb, + & sl_the, sl_psi, sl_di0, sl_fl0, sl_flu) + !!!!!!!!!!!!!!!!!!!!!!!!!! + fluxsurf_sw(ig,l) = sl_flu + !! sl_ls = 180.*zls/pi + !! sl_lct = ptime*24. + 180.*long(ig)/pi/15. + !! sl_lat = 180.*lati(ig)/pi + !! sl_tau = tau(ig,1) + !! sl_alb = albedo(ig,l) + !! sl_the = theta_sl(ig) + !! sl_psi = psi_sl(ig) + !! sl_fl0 = fluxsurf_sw(ig,l) + !! CALL param_slope_full(sl_ls, sl_lct, sl_lat, + !! & sl_tau, sl_alb, + !! & sl_the, sl_psi, sl_fl0, sl_flu) + ENDDO + !!! compute correction on IR flux as well + sky= (1.+cos(pi*theta_sl(ig)/180.))/2. + fluxsurf_lw(ig)= fluxsurf_lw(ig)*sky + ENDIF + ENDDO +ccccc +ccccc PARAM SLOPE +ccccc +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + +c CO2 near infrared absorption +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + call zerophys(ngrid*nlayer,zdtnirco2) + if (callnirco2) then + call nirco2abs (ngrid,nlayer,pplay,dist_sol, + . mu0,fract,declin, zdtnirco2) + endif + + +c Radiative flux from the sky absorbed by the surface (W.m-2) +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + DO ig=1,ngrid + fluxrad_save(ig)=emis(ig)*fluxsurf_lw(ig) + $ +fluxsurf_sw(ig,1)*(1.-albedo(ig,1)) + $ +fluxsurf_sw(ig,2)*(1.-albedo(ig,2)) + + !print*,'RAD ', fluxrad_save(ig) + !print*,'LW ', emis(ig)*fluxsurf_lw(ig) + !print*,'SW ', fluxsurf_sw(ig,2)*(1.-albedo(ig,2)) + + ENDDO + + + +c Net atmospheric radiative heating rate (K.s-1) +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + IF(callnlte) THEN + CALL blendrad(ngrid, nlayer, pplay, + & zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad) + ELSE + DO l=1,nlayer + DO ig=1,ngrid + dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l) + & +zdtnirco2(ig,l) + ENDDO + ENDDO + ENDIF + + !PRINT*,'zdtsw',zdtsw + !PRINT*,'zdtlw',zdtlw + !PRINT*,'zdtnirco2',zdtnirco2 + !PRINT*,'dtrad',dtrad + + + ENDIF ! mod(icount-1,iradia).eq.0 + +c Transformation of the radiative tendencies: +c ----------------------------------------------- + +c Net radiative surface flux (W.m-2) +c ~~~~~~~~~~~~~~~~~~~~~~~~~~ + DO ig=1,ngrid + zplanck(ig)=tsurf(ig)*tsurf(ig) + zplanck(ig)=emis(ig)* + $ stephan*zplanck(ig)*zplanck(ig) + fluxrad(ig)=fluxrad_save(ig)-zplanck(ig) +cccc +cccc param slope +cccc + sky= (1.+cos(pi*theta_sl(ig)/180.))/2. + fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig) +cccc +cccc +cccc + ENDDO + + + DO l=1,nlayer + DO ig=1,ngrid + pdt(ig,l)=pdt(ig,l)+dtrad(ig,l) + ENDDO + ENDDO + + ENDIF + + + +!c----------------------------------------------------------------------- +!c 3. Gravity wave and subgrid scale topography drag : +!c ------------------------------------------------- +! +! +! IF(calllott)THEN +! +! CALL calldrag_noro(ngrid,nlayer,ptimestep, +! & pplay,pplev,pt,pu,pv,zdtgw,zdugw,zdvgw) +! +! DO l=1,nlayer +! DO ig=1,ngrid +! pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l) +! pdu(ig,l)=pdu(ig,l)+zdugw(ig,l) +! pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l) +! ENDDO +! ENDDO +! ENDIF + +c----------------------------------------------------------------------- +c 4. Vertical diffusion (turbulent mixing): +c ----------------------------------------- +c + + IF(calldifv) THEN + + DO ig=1,ngrid + zflubid(ig)=fluxrad(ig)+fluxgrd(ig) + ENDDO + + CALL zerophys(ngrid*nlayer,zdum1) + CALL zerophys(ngrid*nlayer,zdum2) + do l=1,nlayer + do ig=1,ngrid + zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l) + enddo + enddo + + +c Calling vdif (Martian version WITH CO2 condensation) + CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk, + $ ptimestep,capcal,lwrite, + $ pplay,pplev,zzlay,zzlev,z0, + $ pu,pv,zh,pq,tsurf,emis,qsurf, + $ zdum1,zdum2,zdh,pdq,zflubid, + $ zdudif,zdvdif,zdhdif,zdtsdif,q2, + & zdqdif,zdqsdif) + + DO ig=1,ngrid + !! sensible heat flux in W/m2 + hfx(ig) = zflubid(ig)-capcal(ig)*zdtsdif(ig) + !! u star in similarity theory in m/s + ust(ig) = 0.4 + . * sqrt( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) ) + . / log( 1.E+0 + zzlay(ig,1)/z0 ) + ENDDO + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! LES LES + IF (flag_LES) THEN + + write (*,*) '************************************************' + write (*,*) '** LES mode: the difv part is only used to' + write (*,*) '** provide HFX and UST to the dynamics' + write (*,*) '** NB: - dudif, dvdif, dhdif, dqdif are set to 0' + write (*,*) '** - tsurf is updated' + write (*,*) '************************************************' + +! DO ig=1,ngrid +! !! sensible heat flux in W/m2 +! hfx(ig) = zflubid(ig)-capcal(ig)*zdtsdif(ig) +! !! u star in similarity theory in m/s +! ust(ig) = 0.4 +! . * sqrt( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) ) +! . / log( 1.E+0 + zzlay(ig,1)/z0 ) +! + DO l=1,nlayer + zdvdif(ig,l) = 0. + zdudif(ig,l) = 0. + zdhdif(ig,l) = 0. + DO iq=1, nq + zdqdif(ig,l,iq) = 0. + zdqsdif(ig,iq) = 0. !! sortir de la boucle + ENDDO + ENDDO +! +! ENDDO + !write (*,*) 'RAD ',fluxrad(igout) + !write (*,*) 'GRD ',fluxgrd(igout) + !write (*,*) 'dTs/dt ',capcal(igout)*zdtsurf(igout) + !write (*,*) 'HFX ', hfx(igout) + !write (*,*) 'UST ', ust(igout) + ENDIF +!!! LES LES +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + DO l=1,nlayer + DO ig=1,ngrid + + pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l) + pdu(ig,l)=pdu(ig,l)+zdudif(ig,l) + pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l) + + zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only + + ENDDO + ENDDO + + DO ig=1,ngrid + zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig) + ENDDO +!!!!gros caca : sans chaleur sensible +! DO ig=1,ngrid +! zdtsurf(ig)=zdtsurf(ig)+ +! s (fluxrad(ig)+fluxgrd(ig))/capcal(ig) +! ENDDO + + + if(tracer) then + DO iq=1, nq + DO l=1,nlayer + DO ig=1,ngrid + pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq) + ENDDO + ENDDO + ENDDO + DO iq=1, nq + DO ig=1,ngrid + dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq) + ENDDO + ENDDO + end if + + ELSE + DO ig=1,ngrid + zdtsurf(ig)=zdtsurf(ig)+ + s (fluxrad(ig)+fluxgrd(ig))/capcal(ig) + ENDDO +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + IF (flag_LES) THEN + write(*,*) 'LES mode !' + write(*,*) 'Please set calldifv to T in callphys.def' + STOP + ENDIF +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ENDIF + + +c----------------------------------------------------------------------- +c 5. Dry convective adjustment: +c ----------------------------- + + IF(calladj) THEN + + DO l=1,nlayer + DO ig=1,ngrid + zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l) + ENDDO + ENDDO + CALL zerophys(ngrid*nlayer,zduadj) + CALL zerophys(ngrid*nlayer,zdvadj) + CALL zerophys(ngrid*nlayer,zdhadj) + + CALL convadj(ngrid,nlayer,nq,ptimestep, + $ pplay,pplev,zpopsk, + $ pu,pv,zh,pq, + $ pdu,pdv,zdh,pdq, + $ zduadj,zdvadj,zdhadj, + $ zdqadj) + + + DO l=1,nlayer + DO ig=1,ngrid + pdu(ig,l)=pdu(ig,l)+zduadj(ig,l) + pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l) + pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l) + + zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only + ENDDO + ENDDO + + if(tracer) then + DO iq=1, nq + DO l=1,nlayer + DO ig=1,ngrid + pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq) + ENDDO + ENDDO + ENDDO + end if + ENDIF + +c----------------------------------------------------------------------- +c 6. Carbon dioxide condensation-sublimation: +c ------------------------------------------- + + IF(callcond) THEN + CALL newcondens(ngrid,nlayer,nq,ptimestep, + $ capcal,pplay,pplev,tsurf,pt, + $ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, + $ co2ice,albedo,emis, + $ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc, + $ fluxsurf_sw) + + DO l=1,nlayer + DO ig=1,ngrid + pdt(ig,l)=pdt(ig,l)+zdtc(ig,l) + pdv(ig,l)=pdv(ig,l)+zdvc(ig,l) + pdu(ig,l)=pdu(ig,l)+zduc(ig,l) + ENDDO + ENDDO + DO ig=1,ngrid + zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig) + ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep + ENDDO + + IF(tracer) THEN + DO iq=1, nq + DO l=1,nlayer + DO ig=1,ngrid + pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq) + ENDDO + ENDDO + ENDDO + END IF !(tracer) + + ENDIF !(callcond) + +c print*,'condens ok' +c----------------------------------------------------------------------- +c 7. Specific parameterizations for tracers +c: ----------------------------------------- + + if(tracer) then + +c 7a. Water and ice +c --------------- + +c --------------------------------------- +c Water ice condensation in the atmosphere +c ---------------------------------------- + IF (water) THEN + + call watercloud(ngrid,nlayer, ptimestep, + & pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt, + & pq,pdq,zdqcloud,qsurf,zdqscloud,zdtcloud, + & nq,naerkind,tau,icount,zls) + + if (activice) then +c Temperature variation due to latent heat release + DO l=1,nlayer + DO ig=1,ngrid + pdt(ig,l)=pdt(ig,l)+zdtcloud(ig,l) + ENDDO + ENDDO + endif + + IF (iceparty) then + iqmin=nq-1 + ELSE + iqmin=nq + ENDIF + + DO iq=iqmin,nq + DO l=1,nlayer + DO ig=1,ngrid + pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqcloud(ig,l,iq) + ENDDO + ENDDO + DO ig=1,ngrid + dqsurf(ig,iq)= dqsurf(ig,iq) + zdqscloud(ig,iq) + ENDDO + ENDDO + + END IF ! (water) + +c 7b. Chemical species +c ------------------ + +!c -------------- +!c photochemistry : +!c -------------- +! IF(photochem .or. thermochem) then +! call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0, +! . zzlay,zday,pq,pdq,zdqchim,zdqschim,zdqcloud,zdqscloud) +! +!c Photochemistry includes condensation of H2O2 +! +! do iq=nqchem_min,nq +! if (noms(iq).eq."h2o2") then +! DO l=1,nlayer +! DO ig=1,ngrid +! pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqchim(ig,l,iq) +! pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqcloud(ig,l,iq) +! ENDDO +! ENDDO +! else +! DO l=1,nlayer +! DO ig=1,ngrid +! pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqchim(ig,l,iq) +! ENDDO +! ENDDO +! endif +! ENDDO +! do iq=nqchem_min,nq +! if (noms(iq).eq."h2o2") then +! DO ig=1,ngrid +! dqsurf(ig,iq)= dqsurf(ig,iq) + zdqschim(ig,iq) +! dqsurf(ig,iq)= dqsurf(ig,iq) + zdqscloud(ig,iq) +! ENDDO +! else +! DO ig=1,ngrid +! dqsurf(ig,iq)= dqsurf(ig,iq) + zdqschim(ig,iq) +! ENDDO +! endif +! ENDDO +! +! END IF ! (photochem.or.thermochem) +!c print*,'photochem ok' + +c 7c. Aerosol particles +c ------------------- + +c ---------- +c Dust devil : +c ---------- + IF(callddevil) then + call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2, + & zdqdev,zdqsdev) + + if (dustbin.ge.1) then + do iq=1,nq + DO l=1,nlayer + DO ig=1,ngrid + pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq) + ENDDO + ENDDO + enddo + do iq=1,nq + DO ig=1,ngrid + dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq) + ENDDO + enddo + endif ! (test sur dustbin) + + END IF + +c ------------- +c Sedimentation : acts also on water ice +c ------------- + IF (sedimentation) THEN + call zerophys(ngrid*nlayer*nq, zdqsed) + call zerophys(ngrid*nq, zdqssed) + + if(doubleq) then + call callsedim2q(ngrid,nlayer, ptimestep, + & pplev,zzlev, pt, + & pq, pdq, zdqsed, zdqssed,nq) + else + call callsedim(ngrid,nlayer, ptimestep, + & pplev,zzlev, pt, + & pq, pdq, zdqsed, zdqssed,nq) + end if + + + DO iq=1, nq + DO l=1,nlayer + DO ig=1,ngrid + pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq) + ENDDO + ENDDO + ENDDO + DO iq=1, nq + DO ig=1,ngrid + dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq) + ENDDO + ENDDO + END IF ! (sedimentation) + +c print*,'sedim ok' + +c 7d. Updates +c --------- + + DO iq=1, nq + DO ig=1,ngrid + +c --------------------------------- +c Updating tracer budget on surface +c --------------------------------- + qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq) + + ENDDO ! (ig) + ENDDO ! (iq) + + END IF ! (tracer) + +c print*,'tracers ok' + + +!c----------------------------------------------------------------------- +!c 8.5 THERMOSPHERE CALCULATION +!c----------------------------------------------------------------------- +! +! if (callthermos) then +! call thermosphere(pplev,pplay,dist_sol, +! $ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay, +! & pt,pq,pu,pv,pdt,pdq, +! $ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff) +!c do iq=nqchem_min,nq +!c write(*,*) 'thermo iq,pq',iq,pq(690,1,iq) +!c enddo +! +! DO l=1,nlayer +! DO ig=1,ngrid +! dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l) +! pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l) +! & +zdteuv(ig,l) +! pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l) +! pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l) +! DO iq=1, nq +! pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq) +! ENDDO +! ENDDO +! ENDDO +! +! +! endif + +c----------------------------------------------------------------------- +c 8. Surface and sub-surface soil temperature +c----------------------------------------------------------------------- +c +c +c Surface temperature incrementation : +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + DO ig=1,ngrid + tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig) + ENDDO + + +c Prescription piege froid au pole sud (Except at high obliquity !!) +c temperature en surface = temperature equilibre de phases du CO2 + + IF (tracer.AND.water.AND.ngridmx.NE.1) THEN + !if (caps.and.(obliquit.lt.27.)) then + ! tsurf(ngrid)=1/(1/136.27-r/5.9e+5*alog(0.0095*ps(ngrid))) + !endif +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!! note WRF MESOSCALE AYMERIC -- mot cle "caps" +!!!!! watercaptag n'est plus utilise que dans vdifc +!!!!! ... pour que la sublimation ne soit pas stoppee +!!!!! ... dans la calotte permanente nord si qsurf=0 +!!!!! on desire garder cet effet regle par caps=T +!!!!! on a donc commente "if (caps.and.(obliquit.lt.27.))" ci-dessus +!!!!! --- remplacer ces lignes par qqch de plus approprie +!!!!! si on s attaque a la calotte polaire sud +!!!!! pas d'autre occurrence majeure du mot-cle "caps" +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +c ------------------------------------------------------------- +c Change of surface albedo (set to 0.4) in case of ground frost +c everywhere except on the north permanent cap and in regions +c covered by dry ice. +c ALWAYS PLACE these lines after newcondens !!! +c ------------------------------------------------------------- + do ig=1,ngrid + +c -------------- Version temporaire fit TES 2008 ------------ + if (co2ice(ig).eq.0.and.qsurf(ig,nqmx).gt.0.005) then + albedo(ig,1)=0.45 + albedo(ig,2)=0.45 + endif + +c if (co2ice(ig).eq.0.and.qsurf(ig,nqmx).gt.0.005) then +c if (.not.watercaptag(ig)) then +c albedo(ig,1)=0.4 +c albedo(ig,2)=0.4 +c endif +c endif +c -------------- version Francois --------------- +c if (co2ice(ig).eq.0.and. +c & ((qsurf(ig,nqmx).gt.0.005).or.(watercaptag(ig)))) then +c albedo(ig,1)=max(0.4,albedodat(ig)) +c albedo(ig,2)=albedo(ig,1) +c endif +c --------------------------------------------- + enddo ! of do ig=1,ngrid + ENDIF ! of IF (tracer.AND.water.AND.ngridmx.NE.1) + + +c print*,'tracer, water and 3D ok' + +c +c soil temperatures and subsurface heat flux: +c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + IF (callsoil) THEN + CALL soil(ngrid,nsoilmx,.false.,inertiedat, + s ptimestep,tsurf,tsoil,capcal,fluxgrd) + ENDIF + +c print*,'soil ok' + + + + + +c----------------------------------------------------------------------- +c 9. Writing output files +c ------------------------ + +c ------------------------------- +c Dynamical fields incrementation +c ------------------------------- +c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics) + + DO l=1,nlayer + DO ig=1,ngrid + zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep + zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep + zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep + ENDDO + ENDDO + DO iq=1, nq + DO l=1,nlayer + DO ig=1,ngrid + zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep + ENDDO + ENDDO + ENDDO + DO ig=1,ngrid + ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep !already in 7 + ENDDO + DO l=1,nlayer + DO ig=1,ngrid + zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig) + zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig) + ENDDO + ENDDO + ! Density calculation + DO l=1,nlayer + DO ig=1,ngrid + rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l)) + ENDDO + ENDDO + +c Sum of fluxes in solar spectral bands (for output only) + DO ig=1,ngrid + fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2) + fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2) + ENDDO +c ******* TEMPORAIRE ************************************************ + ztim1 = 999 + DO l=1,nlayer + DO ig=1,ngrid + if (pt(ig,l).lt.ztim1) then + ztim1 = pt(ig,l) + igmin = ig + lmin = l + end if + ENDDO + ENDDO + if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then + write(*,*) 'stability WARNING :' + write(*,*) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin), + & 'ig l =', igmin, lmin + end if +c ******************************************************************* + +c -------------------- +c Output on the screen +c -------------------- + + + IF (lwrite) THEN + PRINT*,'Global diagnostics for the physics' + PRINT*,'Variables and their increments x and dx/dt * dt' + WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps' + WRITE(*,'(2f10.5,2f15.5)') + s tsurf(igout),zdtsurf(igout)*ptimestep, + s pplev(igout,1),pdpsrf(igout)*ptimestep + WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT' + WRITE(*,'(i4,6f10.5)') (l, + s pu(igout,l),pdu(igout,l)*ptimestep, + s pv(igout,l),pdv(igout,l)*ptimestep, + s pt(igout,l),pdt(igout,l)*ptimestep, + s l=1,nlayer) + ENDIF + +cc****WRF +cc IF (ngrid.NE.1) THEN +c PRINT*,'check - values after update at ngrid/2+1' +c WRITE(*,'(a4,a6,2a10)') 'l','u','v','T' +c WRITE(*,'(i4,3f10.5)') (l, +c s zu(igout,l), +c s zv(igout,l), +c s zt(igout,l), +c s l=1,nlayer) +c +cc****WRF + + print*,'Ls =',zls*180./pi, + & ' tauref(700 Pa,lat=0) =',tauref(ngrid/2) +c & ' tau(Viking1) =',tau(ig_vl1,1) + + +c ------------------------------------------------------------------- +c Writing NetCDF file "RESTARTFI" at the end of the run +c ------------------------------------------------------------------- +c Remarque : On stocke restarfi +c juste avant que la dynamique ne le soit dans restart. +c entre maintenant et l'ecriture de restart, +c on aura itau = itau +1 et remise a jour de time. +c (lastcall = .true. lorsque itau+1 = itaufin) +c Donc on stocke avec time = time + dtvr + + IF(lastcall) THEN + ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec) + write(*,*)'pour physdem ztime_fin =',ztime_fin + call physdem1("restartfi.nc",long,lati,nsoilmx,nq, + . ptimestep,pday, + . ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf, + . area,albedodat,inertiedat,zmea,zstd,zsig, + . zgam,zthe) + ENDIF + + +ccc**************** WRF OUTPUT ************************** +ccc**************** WRF OUTPUT ************************** +ccc**************** WRF OUTPUT ************************** + do ig=1,ngrid + wtsurf(ig) = tsurf(ig) !! surface temperature + wco2ice(ig) = co2ice(ig) !! co2 ice + + !!! TEMP TEMP TEMP TEMP TEMP TEMP TEMP + !!! specific to WRF WRF WRF + !!! just to output water ice on surface + !!! [it might not be water ice on surface but OK] + !!! uncomment the Registry entry + qsurflast(ig) = qsurf(ig,nqmx) + + enddo +c +c THIS INCLUDE IS AUTOMATICALLY GENERATED FROM REGISTRY +c +#include "fill_save.inc" +c +ccc**************** WRF OUTPUT ************************** +ccc**************** WRF OUTPUT ************************** +ccc**************** WRF OUTPUT ************************** + + +cc----------------------------------- +cc you can still use meso_WRITEDIAGFI (e.g. for debugging purpose), +cc though this is not the default strategy now +cc----------------------------------- +cc please use cudt in namelist.input to set frequency of outputs +cc----------------------------------- +cc BEWARE: if at least one call to meso_WRITEDIAGFI is performed, +cc cudt cannot be 0 - otherwise you'll get a "Floating exception" +cc----------------------------------- +! call meso_WRITEDIAGFI(ngrid,"tauref", +! . "tauref","W.m-2",2, +! . tauref) +! call meso_WRITEDIAGFI(ngrid,"dtrad", +! . "dtrad","W.m-2",2, +! . dtrad) +c call meso_WRITEDIAGFI(ngrid,"tsurf", +c . "tsurf","K",2, +c . tsurf) +c +! call meso_WRITEDIAGFI(ngrid,"zt", +! . "zt","W.m-2",3, +! . zt) +! call meso_WRITEDIAGFI(ngrid,"zdtlw", +! . "zdtlw","W.m-2",2, +! . zdtlw) +! call meso_WRITEDIAGFI(ngrid,"zdtsw", +! . "zdtsw","W.m-2",2, +! . zdtsw) + + + icount=icount+1 + +!!!!!!!!!!!!!TEST TEST +! ENDIF +!!!!!!!!!!!!!TEST TEST + + write(*,*) 'now, back to the dynamical core...' + RETURN + END Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/testphys1d.F =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/testphys1d.F (revision 235) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd/libf/phymars/testphys1d.F (revision 235) @@ -0,0 +1,508 @@ + + PROGRAM testphys1d + IMPLICIT NONE + +c======================================================================= +c subject: +c -------- +c PROGRAM useful to run physical part of the martian GCM in a 1D column +c +c Can be compiled with a command like (e.g. for 25 layers) +c "makegcm -p mars -d 25 testphys1d" +c It requires the files "testphys1d.def" "callphys.def" +c and a file describing the sigma layers (e.g. "z2sig.def") +c +c author: Frederic Hourdin, R.Fournier,F.Forget +c ------- +c +c update: 12/06/2003 including chemistry (S. Lebonnois) +c and water ice (F. Montmessin) +c +c======================================================================= + +#include "dimensions.h" +#include "dimphys.h" +#include "dimradmars.h" +#include "comgeomfi.h" +#include "surfdat.h" +#include "comdiurn.h" +#include "callkeys.h" +#include "comcstfi.h" +#include "planete.h" +#include "comsaison.h" +#include "yomaer.h" +#include "aerdust.h" +#include "control.h" +#include "comvert.h" +#include "netcdf.inc" +#include "comg1d.h" +#include "watercap.h" +#include "fisice.h" +#include "logic.h" + +c -------------------------------------------------------------- +c Declarations +c -------------------------------------------------------------- +c + INTEGER unit ! unite de lecture de "testphys1d.def" + INTEGER unitstart ! unite d'ecriture de "startfi" + INTEGER nlayer,nlevel,nsoil,ndt + INTEGER ilayer,ilevel,isoil,idt,iq + LOGICAl firstcall,lastcall +c + INTEGER day0 ! date initial (sol ; =0 a Ls=0) + REAL day ! date durant le run + REAL time ! time (0 + area(1)=1.E+0 + +c le geopotentiel au sol est inutile en 1D car tout est controle +c par la pression de surface ---> + phisfi(1)=0.E+0 + +c "inifis" reproduit un certain nombre d'initialisations deja faites +c + lecture des clefs de callphys.def +c + calcul de la frequence d'appel au rayonnement +c + calcul des sinus et cosinus des longitude latitude + +!Mars possible matter with dtphys in input and include!!! + CALL inifis(1,llm,day0,daysec,dtphys, + . lati,long,area,rad,g,r,cpp) +c Initialisation pour prendre en compte les vents en 1-D +c ------------------------------------------------------ + ptif=2.E+0*omeg*sinlat(1) + +c vent geostrophique + PRINT *,'composante vers l est du vent geostrophique (U) ?' + READ(unit,*) gru + PRINT *,'composante vers le nord du vent geostrophique (V) ?' + READ(unit,*) grv + +c Initialisation des vents au premier pas de temps + DO ilayer=1,nlayer + u(ilayer)=gru + v(ilayer)=grv + ENDDO + +c energie cinetique turbulente + DO ilevel=1,nlevel + q2(ilevel)=0.E+0 + ENDDO + +c glace de CO2 au sol +c ------------------- + co2ice=0.E+0 + PRINT *,'co2ice (kg.m-2)' + READ(unit,*) co2ice + +c +c emissivite +c ---------- + emis=emissiv + IF (co2ice.eq.1.E+0) THEN + emis=emisice(1) + IF(lati(1).LT.0) emis=emisice(2) + ENDIF + + + +c calcul des pressions et altitudes en utilisant les niveaux sigma +c ---------------------------------------------------------------- + +c Vertical Coordinates +c """""""""""""""""""" + PRINT *,'Hybrid coordinates ?' + READ(unit,*) hybrid + PRINT *,'Hybrid =', hybrid + + CALL disvert + + DO ilevel=1,nlevel + plev(ilevel)=ap(ilevel)+psurf*bp(ilevel) + ENDDO + + DO ilayer=1,nlayer + play(ilayer)=aps(ilayer)+psurf*bps(ilayer) + ENDDO + + DO ilayer=1,nlayer + zlay(ilayer)=-200.E+0 *r*log(play(ilayer)/plev(1)) + & /g + ENDDO + + +c profil de temperature au premier appel +c -------------------------------------- + pks=psurf**rcp + +c altitude en km dans profile: on divise zlay par 1000 + tmp1(0)=0.E+0 + DO ilayer=1,nlayer + tmp1(ilayer)=zlay(ilayer)/1000.E+0 + ENDDO + call profile(unit,nlayer+1,tmp1,tmp2) + + tsurf=tmp2(0) + DO ilayer=1,nlayer + temp(ilayer)=tmp2(ilayer) + ENDDO + + + +c temperature du sous-sol +c ~~~~~~~~~~~~~~~~~~~~~~~ + DO isoil=1,nsoil + tsoil(isoil)=tsurf + ENDDO + +c Initialisation des traceurs +c --------------------------- + + DO iq=1,nqmx + DO ilayer=1,nlayer + q(ilayer,iq) = 0. + ENDDO + ENDDO + + if (photochem.or.callthermos) then + write(*,*) 'Especes chimiques initialisees' + ! thermo=0: initialisation pour toutes les couches + thermo=0 + call inichim_newstart(q,psurf,sig,nqmx,lati,long,area, + $ thermo) + endif + watercaptag(ngridmx)=.false. + + DO iq=1,nqmx-1 + qsurf(iq) = 0. + ENDDO + + +c Initialisation pour les sorties GRADS dans "g1d.dat" et "g1d.ctl" +c ---------------------------------------------------------------- +c (effectuee avec "writeg1d" a partir de "physiq.F" +c ou tout autre subroutine physique, y compris celle ci). + + g1d_nlayer=nlayer + g1d_nomfich='g1d.dat' + g1d_unitfich=40 + g1d_nomctl='g1d.ctl' + g1d_unitctl=41 + g1d_premier=.true. + g2d_premier=.true. + +c Ecriture de "startfi" +c -------------------- +c (Ce fichier sera aussitot relu au premier +c appel de "physiq", mais il est necessaire pour passer +c les variables purement physiques a "physiq"... + + call physdem1("startfi.nc",long,lati,nsoilmx,nqmx, + . dtphys,float(day0), + . time,tsurf,tsoil,co2ice,emis,q2,qsurf, + . area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe) +c======================================================================= +c BOUCLE TEMPORELLE DU MODELE 1D +c======================================================================= +c + firstcall=.true. + lastcall=.false. + + DO idt=1,ndt +c IF (idt.eq.ndt) lastcall=.true. + IF (idt.eq.ndt-day_step-1) then !test + lastcall=.true. + call solarlong(day*1.0,zls) + write(103,*) 'Ls=',zls*180./pi + write(103,*) 'Lat=', lati(1)*180./pi + write(103,*) 'Tau=', tauvis/700*psurf + write(103,*) 'RunEnd - Atmos. Temp. File' + write(103,*) 'RunEnd - Atmos. Temp. File' + write(104,*) 'Ls=',zls*180./pi + write(104,*) 'Lat=', lati(1) + write(104,*) 'Tau=', tauvis/700*psurf + write(104,*) 'RunEnd - Atmos. Temp. File' + ENDIF + +c calcul du geopotentiel +c ~~~~~~~~~~~~~~~~~~~~~ + DO ilayer=1,nlayer + s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp + h(ilayer)=cpp*temp(ilayer)/(pks*s(ilayer)) + ENDDO + phi(1)=pks*h(1)*(1.E+0-s(1)) + DO ilayer=2,nlayer + phi(ilayer)=phi(ilayer-1)+ + & pks*(h(ilayer-1)+h(ilayer))*.5E+0 + & *(s(ilayer-1)-s(ilayer)) + + ENDDO + +c appel de la physique +c -------------------- + + CALL physiq (1,llm,nqmx, + , firstcall,lastcall, + , day,time,dtphys, + , plev,play,phi, + , u, v,temp, q, + , w, +C - sorties + s du, dv, dtemp, dq,dpsurf,tracerdyn) + +c evolution du vent : modele 1D +c ----------------------------- + +c la physique calcule les derivees temporelles de u et v. +c on y rajoute betement un effet Coriolis. +c +c DO ilayer=1,nlayer +c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv) +c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru) +c ENDDO + +c Pour certain test : pas de coriolis a l'equateur +c if(lati(1).eq.0.) then + DO ilayer=1,nlayer + du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4 + dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4 + ENDDO +c end if +c +c +c Calcul du temps au pas de temps suivant +c --------------------------------------- + firstcall=.false. + time=time+dtphys/daysec + IF (time.gt.1.E+0) then + time=time-1.E+0 + day=day+1 + ENDIF + +c calcul des vitesses et temperature au pas de temps suivant +c ---------------------------------------------------------- + + DO ilayer=1,nlayer + u(ilayer)=u(ilayer)+dtphys*du(ilayer) + v(ilayer)=v(ilayer)+dtphys*dv(ilayer) + temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer) + ENDDO + +c calcul des pressions au pas de temps suivant +c ---------------------------------------------------------- + + psurf=psurf+dtphys*dpsurf ! evolution de la pression de surface + DO ilevel=1,nlevel + plev(ilevel)=ap(ilevel)+psurf*bp(ilevel) + ENDDO + DO ilayer=1,nlayer + play(ilayer)=aps(ilayer)+psurf*bps(ilayer) + ENDDO + + ENDDO ! fin de la boucle temporelle + +c ======================================================== +c GESTION DES SORTIE +c ======================================================== + +c fermeture pour conclure les sorties format grads dans "g1d.dat" +c et "g1d.ctl" (effectuee avec "writeg1d" a partir de "physiq.F" +c ou tout autre subroutine physique, y compris celle ci). + +c CALL endg1d(1,nlayer,zphi/(g*1000.),ndt) + CALL endg1d(1,nlayer,zlay/1000.,ndt) + +c ======================================================== + END + +c*********************************************************************** +c*********************************************************************** +c Subroutines Bidons utilise seulement en 3D, mais +c necessaire a la compilation de testphys1d en 1-D + + subroutine gr_fi_dyn + RETURN + END + + subroutine iniwrite + RETURN + END + +c*********************************************************************** +c*********************************************************************** + +#include "../dyn3d/disvert.F"