Context Navigation

meso_physiq.F @ 119

Last change on this file since 119 was 86, checked in by aslmd, 14 years ago

*
mars + LMD_MM_MARS * Precompilation flag MESOSCALE for better transparency

* in shared phymars between GCM and mesoscale model

M 85 mars/libf/phymars/meso_physiq.F
M 85 mars/libf/phymars/meso_inifis.F
Added a pre-compilation flag MESOSCALE so that the LMDZ.MARS GCM
will compile without stating errors because of mesoscale routines.

M 85 mars/libf/phymars/newcondens.F
M 85 mars/libf/phymars/testphys1d.F
M 85 mars/libf/phymars/dustlift.F
D 85 mars/libf/phymars/meso_testphys1d.F
D 85 mars/libf/phymars/meso_dustlift.F
D 85 mars/libf/phymars/meso_newcondens.F
Now, this MESOSCALE precompilation flag can be used to lower
the number of meso_* routines when adaptations for mesoscale
applications are not very extended.
--> Three meso_* routines were deleted and changes are
now impacted under the MESOSCALE flag in the original GCM routines
--> Completely transparent for GCM compilation since it is devoid
of the -DMESOSCALE option
--> Very good for syncing because changes in dustlift, newcondens
will be directly available in the mesoscale model

M 84 mesoscale/LMD_MM_MARS/makemeso
Changed meso_testphys1d in testphys1d

M 84 mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_pgf
M 84 mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_mpifort
M 84 mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_ifort
M 84 mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_g95
M 84 mesoscale/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new/makegcm_mpi
Added the option -DMESOSCALE in these scripts

*
LMD_MM_MARS * Various minor changes related to water cycle and plotting routines

* Also included the GW test case

A 0 mesoscale/LMDZ.MARS.new/myGCM/DEFS_JB/callphys.def.orig
M 84 mesoscale/NOTES.txt
D 84 mesoscale/LMD_MM_MARS/SRC/ARWpost/idl
M 84 mesoscale/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM
M 84 mesoscale/LMD_MM_MARS/SIMU/gnome_launch.meso
M 85 mesoscale/PLOT/MINIMAL/map_latlon.pro
D 85 mesoscale/PLOT/SPEC/LES/getget.pro
M 85 mesoscale/PLOT/SPEC/MAP/map_uvt.pro
A + - mesoscale/PLOT/SPEC/getget.pro
A 0 mesoscale/PLOT/RESERVE/obsolete
A 0 mesoscale/TESTS/TESTGW.tar.gz
M 84 000-USERS

File size: 67.8 KB

Rev	Line
[42]	1	SUBROUTINE meso_physiq(ngrid,nlayer,nq,
	2	$ firstcall,lastcall,
	3	$ wday_ini,
	4	$ pday,ptime,ptimestep,
	5	$ pplev,pplay,pphi,
[55]	6	$ pu,pv,pt,pq,pw,
	7	$ wtnom,
[42]	8	$ pdu,pdv,pdt,pdq,pdpsrf,tracerdyn,
	9	$ wtsurf,wtsoil,wemis,wq2,wqsurf,wco2ice,
	10	$ wisoil,wdsoil,
	11	$ wecritphys,
[86]	12	#ifdef MESOSCALE
[42]	13	$ output_tab2d, output_tab3d,
[86]	14	#endif
[42]	15	$ flag_LES)
	16
	17
	18
	19	IMPLICIT NONE
	20	c=======================================================================
	21	c
	22	c CAREFUL: THIS IS A VERSION TO BE USED WITH WRF !!!
	23	c
	24	c ... CHECK THE ****WRF lines
	25	c
	26	c=======================================================================
	27	c
	28	c subject:
	29	c --------
	30	c
	31	c Organisation of the physical parametrisations of the LMD
	32	c martian atmospheric general circulation model.
	33	c
	34	c The GCM can be run without or with tracer transport
	35	c depending on the value of Logical "tracer" in file "callphys.def"
	36	c Tracers may be water vapor, ice OR chemical species OR dust particles
	37	c
	38	c SEE comments in initracer.F about numbering of tracer species...
	39	c
	40	c It includes:
	41	c
	42	c 1. Initialization:
	43	c 1.1 First call initializations
	44	c 1.2 Initialization for every call to physiq
	45	c 1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
	46	c 2. Compute radiative transfer tendencies
	47	c (longwave and shortwave) for CO2 and aerosols.
	48	c 3. Gravity wave and subgrid scale topography drag :
	49	c 4. Vertical diffusion (turbulent mixing):
	50	c 5. Convective adjustment
	51	c 6. Condensation and sublimation of carbon dioxide.
	52	c 7. TRACERS :
	53	c 7a. water and water ice
	54	c 7b. call for photochemistry when tracers are chemical species
	55	c 7c. other scheme for tracer (dust) transport (lifting, sedimentation)
	56	c 7d. updates (CO2 pressure variations, surface budget)
	57	c 8. Contribution to tendencies due to thermosphere
	58	c 9. Surface and sub-surface temperature calculations
	59	c 10. Write outputs :
	60	c - "startfi", "histfi" (if it's time)
	61	c - Saving statistics (if "callstats = .true.")
	62	c - Dumping eof (if "calleofdump = .true.")
	63	c - Output any needed variables in "diagfi"
	64	c 11. Diagnostic: mass conservation of tracers
	65	c
	66	c author:
	67	c -------
	68	c Frederic Hourdin 15/10/93
	69	c Francois Forget 1994
	70	c Christophe Hourdin 02/1997
	71	c Subroutine completly rewritten by F.Forget (01/2000)
	72	c Introduction of the photochemical module: S. Lebonnois (11/2002)
	73	c Introduction of the thermosphere module: M. Angelats i Coll (2002)
	74	c Water ice clouds: Franck Montmessin (update 06/2003)
	75	c Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
	76	c Nb: See callradite.F for more information.
	77	c new WRF version: Aymeric Spiga (01/2009)
	78	c
	79	c
	80	c arguments:
	81	c ----------
	82	c
	83	c input:
	84	c ------
	85	c ecri period (in dynamical timestep) to write output
	86	c ngrid Size of the horizontal grid.
	87	c All internal loops are performed on that grid.
	88	c nlayer Number of vertical layers.
	89	c nq Number of advected fields
	90	c firstcall True at the first call
	91	c lastcall True at the last call
	92	c pday Number of days counted from the North. Spring
	93	c equinoxe.
	94	c ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
	95	c ptimestep timestep (s)
	96	c pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa)
	97	c pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
	98	c pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2)
	99	c pu(ngrid,nlayer) u component of the wind (ms-1)
	100	c pv(ngrid,nlayer) v component of the wind (ms-1)
	101	c pt(ngrid,nlayer) Temperature (K)
	102	c pq(ngrid,nlayer,nq) Advected fields
	103	c pudyn(ngrid,nlayer) \
	104	c pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the
	105	c ptdyn(ngrid,nlayer) / corresponding variables
	106	c pqdyn(ngrid,nlayer,nq) /
	107	c pw(ngrid,?) vertical velocity
	108	c
	109	c
	110	c ****WRF
	111	c day_ini,tsurf,tsoil,emis,q2,qsurf,co2ice are inputs
	112	c and locally saved variables
	113	c (no need to call phyetat0)
	114	c
	115	c
	116	c output:
	117	c -------
	118	c
	119	c pdu(ngrid,nlayermx) \
	120	c pdv(ngrid,nlayermx) \ Temporal derivative of the corresponding
	121	c pdt(ngrid,nlayermx) / variables due to physical processes.
	122	c pdq(ngrid,nlayermx,nqmx) /
	123	c pdpsrf(ngrid) /
	124	c tracerdyn call tracer in dynamical part of GCM ?
	125
	126	c
	127	c=======================================================================
	128	c
	129	c 0. Declarations :
	130	c ------------------
	131
	132	#include "dimensions.h"
	133	#include "dimphys.h"
	134	#include "comgeomfi.h"
	135	#include "surfdat.h"
	136	#include "comsoil.h" !!! new soil common
	137	#include "comdiurn.h"
	138	#include "callkeys.h"
	139	#include "comcstfi.h"
	140	#include "planete.h"
	141	#include "comsaison.h"
	142	#include "control.h"
	143	#include "dimradmars.h"
	144	#include "comg1d.h"
	145	#include "tracer.h"
	146	#include "nlteparams.h"
	147
	148	#include "chimiedata.h"
	149	#include "watercap.h"
	150	#include "param.h"
	151	#include "param_v3.h"
	152	#include "conc.h"
	153
	154	#include "netcdf.inc"
	155
[86]	156	!!!!**** SPECIFIC TO MESOSCALE
	157	#ifdef MESOSCALE
[42]	158	#include "meso_slope.h"
	159	#include "wrf_output_2d.h"
	160	#include "wrf_output_3d.h"
[86]	161	#endif
[42]	162
[55]	163	#include "advtrac.h" !!! this is necessary for tracers (in dyn3d)
[42]	164
	165	c Arguments :
	166	c -----------
	167
	168	c inputs:
	169	c -------
	170	INTEGER ngrid,nlayer,nq
	171	REAL ptimestep
	172	REAL pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
	173	REAL pphi(ngridmx,nlayer)
	174	REAL pu(ngridmx,nlayer),pv(ngridmx,nlayer)
	175	REAL pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
	176	REAL pw(ngridmx,nlayer) !Mars pvervel transmit par dyn3d
	177	REAL zh(ngridmx,nlayermx) ! potential temperature (K)
	178	LOGICAL firstcall,lastcall
[86]	179	!!! ****WRF WRF specific to mesoscale
[42]	180	INTEGER wday_ini
	181	REAL wtsurf(ngridmx) ! input only ay firstcall - output
	182	REAL wtsoil(ngridmx,nsoilmx)
	183	REAL wisoil(ngridmx,nsoilmx) !! new soil scheme
	184	REAL wdsoil(ngridmx,nsoilmx) !! new soil scheme
	185	REAL wco2ice(ngridmx)
	186	REAL wemis(ngridmx)
	187	REAL wqsurf(ngridmx,nqmx)
	188	REAL wq2(ngridmx,nlayermx+1)
	189	REAL wecritphys
[86]	190	#ifdef MESOSCALE
[42]	191	REAL output_tab2d(ngridmx,n2d)
	192	REAL output_tab3d(ngridmx,nlayer,n3d)
[86]	193	#endif
[42]	194	REAL sl_ls, sl_lct, sl_lat, sl_tau, sl_alb, sl_the, sl_psi
	195	REAL sl_fl0, sl_flu
	196	REAL sl_ra, sl_di0
	197	REAL sky
	198	REAL hfx(ngridmx) !! pour LES avec isfflx!=0
	199	REAL ust(ngridmx) !! pour LES avec isfflx!=0
	200	LOGICAL flag_LES !! pour LES avec isfflx!=0
[77]	201	REAL qsurfice(ngridmx) !! pour diagnostics
[42]	202	real alpha,lay1 ! coefficients for building layers
	203	integer iloop
[55]	204	INTEGER tracerset !!! this corresponds to config%mars
[86]	205	!!! ****WRF WRF specific to mesoscale
[42]	206	REAL pday
	207	REAL ptime
	208	logical tracerdyn
[55]	209	CHARACTER (len=20) :: wtnom(nqmx) ! tracer name
[42]	210
	211	c outputs:
	212	c --------
	213	c physical tendencies
	214	REAL pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
	215	REAL pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
	216	REAL pdpsrf(ngridmx) ! surface pressure tendency
	217
	218
	219	c Local saved variables:
	220	c ----------------------
	221	c aerosol (dust or ice) extinction optical depth at reference wavelength
	222	c "longrefvis" set in dimradmars.h , for one of the "naerkind" kind of
	223	c aerosol optical properties :
	224	REAL aerosol(ngridmx,nlayermx,naerkind)
	225
	226	INTEGER day_ini ! Initial date of the run (sol since Ls=0)
	227	INTEGER icount ! counter of calls to physiq during the run.
	228	REAL tsurf(ngridmx) ! Surface temperature (K)
	229	REAL tsoil(ngridmx,nsoilmx) ! sub-surface temperatures (K)
	230	REAL co2ice(ngridmx) ! co2 ice surface layer (kg.m-2)
	231	REAL albedo(ngridmx,2) ! Surface albedo in each solar band
	232	REAL emis(ngridmx) ! Thermal IR surface emissivity
	233	REAL dtrad(ngridmx,nlayermx) ! Net atm. radiative heating rate (K.s-1)
	234	REAL fluxrad_sky(ngridmx) ! rad. flux from sky absorbed by surface (W.m-2)
	235	REAL fluxrad(ngridmx) ! Net radiative surface flux (W.m-2)
	236	REAL capcal(ngridmx) ! surface heat capacity (J m-2 K-1)
	237	REAL fluxgrd(ngridmx) ! surface conduction flux (W.m-2)
	238	REAL qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2)
	239	REAL q2(ngridmx,nlayermx+1) ! Turbulent Kinetic Energy
	240	INTEGER ig_vl1 ! Grid Point near VL1 (for diagnostic)
	241
	242	c Variables used by the water ice microphysical scheme:
	243	REAL rice(ngridmx,nlayermx) ! Water ice geometric mean radius (m)
	244	REAL nuice(ngridmx,nlayermx) ! Estimated effective variance
	245	! of the size distribution
	246	c Albedo of deposited surface ice
[82]	247	!!REAL, PARAMETER :: alb_surfice = 0.4 ! 0.45
	248	REAL, PARAMETER :: alb_surfice = 0.45 !!TESTS_JB
[42]	249
	250	SAVE day_ini, icount
	251	SAVE aerosol, tsurf,tsoil
	252	SAVE co2ice,albedo,emis, q2
	253	SAVE capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
	254	SAVE ig_vl1
	255
	256	REAL stephan
	257	DATA stephan/5.67e-08/ ! Stephan Boltzman constant
	258	SAVE stephan
	259
	260	c Local variables :
	261	c -----------------
	262
	263	REAL CBRT
	264	EXTERNAL CBRT
	265
	266	CHARACTER*80 fichier
	267	INTEGER l,ig,ierr,igout,iq,i, tapphys
	268
	269	REAL fluxsurf_lw(ngridmx) !incident LW (IR) surface flux (W.m-2)
	270	REAL fluxsurf_sw(ngridmx,2) !incident SW (solar) surface flux (W.m-2)
	271	REAL fluxtop_lw(ngridmx) !Outgoing LW (IR) flux to space (W.m-2)
	272	REAL fluxtop_sw(ngridmx,2) !Outgoing SW (solar) flux to space (W.m-2)
	273	REAL tauref(ngridmx) ! Reference column optical depth at 700 Pa
	274	! (used if active=F)
	275	REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point
	276	REAL zls ! solar longitude (rad)
	277	REAL zday ! date (time since Ls=0, in martian days)
	278	REAL zzlay(ngridmx,nlayermx) ! altitude at the middle of the layers
	279	REAL zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
	280	REAL latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic)
	281
	282	c Tendancies due to various processes:
	283	REAL dqsurf(ngridmx,nqmx)
	284	REAL zdtlw(ngridmx,nlayermx) ! (K/s)
	285	REAL zdtsw(ngridmx,nlayermx) ! (K/s)
	286	REAL cldtlw(ngridmx,nlayermx) ! (K/s) LW heating rate for clear area
	287	REAL cldtsw(ngridmx,nlayermx) ! (K/s) SW heating rate for clear area
	288	REAL zdtnirco2(ngridmx,nlayermx) ! (K/s)
	289	REAL zdtnlte(ngridmx,nlayermx) ! (K/s)
	290	REAL zdtsurf(ngridmx) ! (K/s)
	291	REAL zdtcloud(ngridmx,nlayermx)
	292	REAL zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx) ! (m.s-2)
	293	REAL zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx) ! (K/s)
	294	REAL zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx) ! (m.s-2)
	295	REAL zdhadj(ngridmx,nlayermx) ! (K/s)
	296	REAL zdtgw(ngridmx,nlayermx) ! (K/s)
	297	REAL zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx) ! (m.s-2)
	298	REAL zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
	299	REAL zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
	300
	301	REAL zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
	302	REAL zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
	303	REAL zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
	304	REAL zdqadj(ngridmx,nlayermx,nqmx)
	305	REAL zdqc(ngridmx,nlayermx,nqmx)
	306	REAL zdqcloud(ngridmx,nlayermx,nqmx)
	307	REAL zdqscloud(ngridmx,nqmx)
	308	REAL zdqchim(ngridmx,nlayermx,nqmx)
	309	REAL zdqschim(ngridmx,nqmx)
	310
	311	REAL zdteuv(ngridmx,nlayermx) ! (K/s)
	312	REAL zdtconduc(ngridmx,nlayermx) ! (K/s)
	313	REAL zdumolvis(ngridmx,nlayermx)
	314	REAL zdvmolvis(ngridmx,nlayermx)
	315	real zdqmoldiff(ngridmx,nlayermx,nqmx)
	316
	317	c Local variable for local intermediate calcul:
	318	REAL zflubid(ngridmx)
	319	REAL zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
	320	REAL zdum1(ngridmx,nlayermx)
	321	REAL zdum2(ngridmx,nlayermx)
	322	REAL ztim1,ztim2,ztim3, z1,z2
	323	REAL ztime_fin
	324	REAL zdh(ngridmx,nlayermx)
	325	INTEGER length
	326	PARAMETER (length=100)
	327
	328	c local variables only used for diagnostic (output in file "diagfi" or "stats")
	329	c -----------------------------------------------------------------------------
	330	REAL ps(ngridmx), zt(ngridmx,nlayermx)
	331	REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
	332	REAL zq(ngridmx,nlayermx,nqmx)
	333	REAL fluxtop_sw_tot(ngridmx), fluxsurf_sw_tot(ngridmx)
	334	character*2 str2
	335	character*5 str5
	336	real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
	337	REAL ccn(ngridmx,nlayermx) ! Cloud condensation nuclei
	338	! (particules kg-1)
[83]	339	SAVE ccn !! in case iradia != 1
[42]	340	real rdust(ngridmx,nlayermx) ! dust geometric mean radius (m)
	341	real qtot1,qtot2 ! total aerosol mass
	342	integer igmin, lmin
	343	logical tdiag
	344
	345	real co2col(ngridmx) ! CO2 column
	346	REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
	347	REAL zstress(ngrid), cd
	348	real hco2(nqmx),tmean, zlocal(nlayermx)
	349	real rho(ngridmx,nlayermx) ! density
	350	real vmr(ngridmx,nlayermx) ! volume mixing ratio
	351	REAL mtot(ngridmx) ! Total mass of water vapor (kg/m2)
	352	REAL icetot(ngridmx) ! Total mass of water ice (kg/m2)
	353	REAL rave(ngridmx) ! Mean water ice effective radius (m)
	354	REAL opTES(ngridmx,nlayermx)! abs optical depth at 825 cm-1
	355	REAL tauTES(ngridmx) ! column optical depth at 825 cm-1
	356	REAL Qabsice ! Water ice absorption coefficient
	357
	358
	359	REAL time_phys
	360
	361	c=======================================================================
[86]	362	#ifdef MESOSCALE
[42]	363
	364	c 1. Initialisation:
	365	c -----------------
	366
	367	c 1.1 Initialisation only at first call
	368	c ---------------------------------------
	369	IF (firstcall) THEN
	370
	371	c variables set to 0
	372	c ~~~~~~~~~~~~~~~~~~
	373	call zerophys(ngridnlayernaerkind,aerosol)
	374	call zerophys(ngrid*nlayer,dtrad)
	375	call zerophys(ngrid,fluxrad)
	376
	377	c read startfi
	378	c ~~~~~~~~~~~~
	379	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	380	c ****WRF
	381	c
	382	c No need to use startfi.nc
	383	c > part of the job of phyetat0 is done in inifis
	384	c > remaining initializations are passed here from the WRF variables
	385	c > beware, some operations were done by phyetat0 (ex: tracers)
	386	c > if any problems, look in phyetat0
	387	c
	388	tsurf(:)=wtsurf(:)
	389	PRINT*,'check: tsurf ',tsurf(1),tsurf(ngridmx)
	390	tsoil(:,:)=wtsoil(:,:)
	391	PRINT*,'check: tsoil ',tsoil(1,1),tsoil(ngridmx,nsoilmx)
	392	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	393	!!!new physics
	394	inertiedat(:,:)=wisoil(:,:)
	395	PRINT*,'check: inert ',inertiedat(1,1),inertiedat(ngridmx,nsoilmx)
	396	mlayer(0:nsoilmx-1)=wdsoil(1,:)
	397	PRINT*,'check: layer ', mlayer
	398	!!!!!!!!!!!!!!!!! DONE in soil_setting.F
	399	! 1.5 Build layer(); following the same law as mlayer()
	400	! Assuming layer distribution follows mid-layer law:
	401	! layer(k)=lay1alpha*(k-1)
	402	lay1=sqrt(mlayer(0)*mlayer(1))
	403	alpha=mlayer(1)/mlayer(0)
	404	do iloop=1,nsoilmx
	405	layer(iloop)=lay1(alpha*(iloop-1))
	406	enddo
	407	!!!!!!!!!!!!!!!!! DONE in soil_setting.F
[55]	408	tnom(:)=wtnom(:) !! est rempli dans advtrac.h
	409	PRINT*,'check: tracernames ', tnom
[42]	410	!!!new physics
	411	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	412	emis(:)=wemis(:)
	413	PRINT*,'check: emis ',emis(1),emis(ngridmx)
	414	q2(:,:)=wq2(:,:)
	415	PRINT*,'check: q2 ',q2(1,1),q2(ngridmx,nlayermx+1)
	416	qsurf(:,:)=wqsurf(:,:)
	417	PRINT*,'check: qsurf ',qsurf(1,1),qsurf(ngridmx,nqmx)
	418	co2ice(:)=wco2ice(:)
	419	PRINT*,'check: co2 ',co2ice(1),co2ice(ngridmx)
	420	day_ini=wday_ini
	421
	422	c artificially filling dyn3d/control.h is also required
	423	c > iphysiq is put in WRF to be set easily (cf ptimestep)
	424	c > day_step is simply deduced:
	425	c
	426	day_step=daysec/ptimestep
	427	PRINT*,'Call to LMD physics:',day_step,' per Martian day'
	428	c
	429	iphysiq=ptimestep
	430	c
	431	ecritphy=wecritphys
	432	PRINT*,'Write LMD physics each:',ecritphy,' seconds'
	433	!!PRINT*,ecri_phys
	434	!!PRINT*,float(ecri_phys) ...
	435	!!renvoient tous deux des nombres absurdes
	436	!!pourtant callkeys.h est inclus ...
	437	!!
	438	!!donc ecritphys est passe en argument ...
	439	PRINT*,'Write LMD physics each:',ecritphy,' seconds'
	440	c
[61]	441	!DO iq=1, nq
	442	! PRINT*, tnom(iq), pq(:,:,iq)
	443	!ENDDO
[80]	444
[55]	445	c
[42]	446	c ****WRF
	447	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	448
	449
	450
	451
	452	!! Read netcdf initial physical parameters.
	453	! CALL phyetat0 ("startfi.nc",0,0,
	454	! & nsoilmx,nq,
	455	! & day_ini,time_phys,
	456	! & tsurf,tsoil,emis,q2,qsurf,co2ice)
	457
	458	if (pday.ne.day_ini) then
	459	write(,) "PHYSIQ: ERROR: bad synchronization between ",
	460	& "physics and dynamics"
	461	write(,) "dynamics day: ",pday
	462	write(,) "physics day: ",day_ini
	463	stop
	464	endif
	465
	466	write (,) 'In physiq day_ini =', day_ini
	467
	468	c Initialize albedo and orbital calculation
	469	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	470	CALL surfini(ngrid,co2ice,qsurf,albedo)
	471	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
	472
	473	c initialize soil
	474	c ~~~~~~~~~~~~~~~
	475	IF (callsoil) THEN
	476	CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
	477	s ptimestep,tsurf,tsoil,capcal,fluxgrd)
	478	ELSE
	479	PRINT*,
	480	& 'PHYSIQ WARNING! Thermal conduction in the soil turned off'
	481	DO ig=1,ngrid
	482	capcal(ig)=1.e5
	483	fluxgrd(ig)=0.
	484	ENDDO
	485	ENDIF
	486	icount=1
	487
	488
	489	c initialize tracers
	490	c ~~~~~~~~~~~~~~~~~~
	491	tracerdyn=tracer
	492	IF (tracer) THEN
	493	CALL initracer(qsurf,co2ice)
	494	ENDIF ! end tracer
	495
[80]	496	!!!!!! WRF WRF WRF MARS MARS
	497	!!!!!! TEST TEST TEST TEST AS+JBM 28/02/11
	498	!!!!!! TEST TEST TEST TEST AS+JBM 28/02/11
	499	!!!!!! TEST TEST TEST TEST AS+JBM 28/02/11
	500	!!!!
	501	!!!! principe: une option 'caps=T' specifique au mesoscale
	502	!!!! ... en vue d'un meso_initracer ????
	503	!!!!
	504	!!!! depots permanents => albedo TES du PDS
	505	!!!! depots saisonniers => alb_surfice (~0.4, cf plus bas)
	506	!!!! [!!!! y compris pour les depots saisonniers sur les depots permanents]
	507	!!!!
	508	!!!! --> todo: il faut garder les depots saisonniers qui viennent
	509	!!!! du GCM lorsqu'ils sont consequents
	510	!!!!
	511	IF ( caps .and. (igcm_h2o_ice .ne. 0) ) THEN
	512	PRINT *, 'OVERWRITING watercaptag DEFINITION in INITRACER'
	513	PRINT *, 'lat>70 et alb>0.26 => watercaptag=T'
	514	!! Perennial H20 north cap defined by watercaptag=true (allows surface to be
	515	!! hollowed by sublimation in vdifc).
	516	do ig=1,ngridmx
	517	qsurf(ig,igcm_h2o_ice)=0. !! on jette les inputs GCM
	518	if ( (lati(ig)*180./pi.gt.70.) .and.
	519	. (albedodat(ig).ge.0.26) ) then
	520	watercaptag(ig)=.true.
	521	dryness(ig) = 1.
	522	else
	523	watercaptag(ig)=.false.
	524	dryness(ig) = 1.
	525	endif ! (lati, albedodat)
	526	end do ! (ngridmx)
	527	ELSE ! (caps)
	528	print *,'Blork !!!'
	529	print *,'caps=T avec water=F ????'
	530	ENDIF ! (caps)
	531	!!!!!! TEST TEST TEST TEST AS+JBM 28/02/11
	532	!!!!!! TEST TEST TEST TEST AS+JBM 28/02/11
	533	!!!!!! TEST TEST TEST TEST AS+JBM 28/02/11
[42]	534
[80]	535
[42]	536	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	537	c ****WRF
	538	c
	539	c nosense in mesoscale modeling
	540	c
	541	cc Determining gridpoint near Viking Lander 1 (used for diagnostic only)
	542	cc ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	543	c
	544	c if(ngrid.ne.1) then
	545	c latvl1= 22.27
	546	c lonvl1= -47.94
	547	c ig_vl1= 1+ int( (1.5-(latvl1-90.)jjm/180.) -2 )iim +
	548	c & int(1.5+(lonvl1+180)*iim/360.)
	549	c write(,) 'Viking Lander 1 GCM point: lat,lon',
	550	c & lati(ig_vl1)180/pi, long(ig_vl1)180/pi
	551	c end if
	552	c ****WRF
	553	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	554
	555	!!!
	556	!!! WRF WRF WRF commented for smaller executables
	557	!!!
	558	!c Initialize thermospheric parameters
	559	!c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	560	!
	561	! if (callthermos) call param_read
	562
	563
	564	c Initialize R and Cp as constant
	565
	566	if (.not.callthermos .and. .not.photochem) then
	567	do l=1,nlayermx
	568	do ig=1,ngridmx
	569	rnew(ig,l)=r
	570	cpnew(ig,l)=cpp
	571	mmean(ig,l)=mugaz
	572	enddo
	573	enddo
	574	endif
	575
	576	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
	577	write(,)"physiq: water_param Surface ice alb:",alb_surfice
	578	ENDIF
	579
	580	ENDIF ! (end of "if firstcall")
	581
	582	c ---------------------------------------------------
	583	c 1.2 Initializations done at every physical timestep:
	584	c ---------------------------------------------------
	585	c
	586	IF (ngrid.NE.ngridmx) THEN
	587	PRINT*,'STOP in PHYSIQ'
	588	PRINT*,'Probleme de dimensions :'
	589	PRINT*,'ngrid = ',ngrid
	590	PRINT*,'ngridmx = ',ngridmx
	591	STOP
	592	ENDIF
	593
	594	c Initialize various variables
	595	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	596	call zerophys(ngrid*nlayer, pdv)
	597	call zerophys(ngrid*nlayer, pdu)
	598	call zerophys(ngrid*nlayer, pdt)
	599	call zerophys(ngridnlayernq, pdq)
	600	call zerophys(ngrid, pdpsrf)
	601	call zerophys(ngrid, zflubid)
	602	call zerophys(ngrid, zdtsurf)
	603	call zerophys(ngrid*nq, dqsurf)
	604	igout=ngrid/2+1
	605
	606
	607	zday=pday+ptime ! compute time, in sols (and fraction thereof)
	608
	609	c Compute Solar Longitude (Ls) :
	610	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	611	if (season) then
	612	call solarlong(zday,zls)
	613	else
	614	call solarlong(float(day_ini),zls)
	615	end if
	616
	617	c Compute geopotential at interlayers
	618	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	619	c ponderation des altitudes au niveau des couches en dp/p
	620
	621	DO l=1,nlayer
	622	DO ig=1,ngrid
	623	zzlay(ig,l)=pphi(ig,l)/g
	624	ENDDO
	625	ENDDO
	626	DO ig=1,ngrid
	627	zzlev(ig,1)=0.
	628	zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km...
	629	ENDDO
	630	DO l=2,nlayer
	631	DO ig=1,ngrid
	632	z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
	633	z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
	634	zzlev(ig,l)=(z1zzlay(ig,l-1)+z2zzlay(ig,l))/(z1+z2)
	635	ENDDO
	636	ENDDO
	637
	638
	639	! Potential temperature calculation not the same in physiq and dynamic
	640
	641	c Compute potential temperature
	642	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	643	DO l=1,nlayer
	644	DO ig=1,ngrid
	645	zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
	646	zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
	647	ENDDO
	648	ENDDO
	649
	650	!!!
	651	!!! WRF WRF WRF commented for smaller executables
	652	!!!
	653	!c-----------------------------------------------------------------------
	654	!c 1.2.5 Compute mean mass, cp, and R
	655	!c --------------------------------
	656	!
	657	! if(photochem.or.callthermos) then
	658	! call concentrations(pplay,pt,pdt,pq,pdq,ptimestep)
	659	! endif
	660
	661	c-----------------------------------------------------------------------
	662	c 2. Compute radiative tendencies :
	663	c------------------------------------
	664
	665
	666	IF (callrad) THEN
	667	IF( MOD(icount-1,iradia).EQ.0) THEN
	668
	669	write (,) 'call radiative transfer'
	670
	671	c Local Solar zenith angle
	672	c ~~~~~~~~~~~~~~~~~~~~~~~~
	673	CALL orbite(zls,dist_sol,declin)
	674
	675	IF(diurnal) THEN
	676	ztim1=SIN(declin)
	677	ztim2=COS(declin)COS(2.pi*(zday-.5))
	678	ztim3=-COS(declin)SIN(2.pi*(zday-.5))
	679
	680	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
	681	s ztim1,ztim2,ztim3, mu0,fract)
	682
	683	ELSE
	684	CALL mucorr(ngrid,declin, lati, mu0, fract,10000.,rad)
	685	ENDIF
	686
	687	c NLTE cooling from CO2 emission
	688	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	689
	690	IF(callnlte) CALL nltecool(ngrid,nlayer,pplay,pt,zdtnlte)
	691
	692	c Find number of layers for LTE radiation calculations
	693	IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
	694	& CALL nlthermeq(ngrid,nlayer,pplev,pplay)
	695
	696	c Note: Dustopacity.F has been transferred to callradite.F
	697
	698	c Call main radiative transfer scheme
	699	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	700	c Transfer through CO2 (except NIR CO2 absorption)
	701	c and aerosols (dust and water ice)
	702
	703	c Radiative transfer
	704	c ------------------
	705	cc
	706	cc **WRF: desormais dust_opacity est dans callradite -- modifications
	707	cc nveaux arguments: tauref,tau,aerosol,rice,nuice
	708	cc
	709	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
	710	$ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
	711	$ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
	712	& tauref,tau,aerosol,ccn,rdust,rice,nuice)
	713
	714	c write(,) icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
	715	c $ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
	716	c $ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
	717	c & tauref,tau,aerosol,rice,nuice
	718	c write(,) fluxsurf_lw
	719
	720	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	721	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	722	ccccc
	723	ccccc PARAM SLOPE : Insolation (direct + scattered)
	724	ccccc
	725	DO ig=1,ngrid
	726	sl_the = theta_sl(ig)
	727	IF (sl_the .ne. 0.) THEN
	728	ztim1=fluxsurf_sw(ig,1)+fluxsurf_sw(ig,2)
	729	DO l=1,2
	730	sl_lct = ptime24. + 180.long(ig)/pi/15.
	731	sl_ra = pi*(1.0-sl_lct/12.)
	732	sl_lat = 180.*lati(ig)/pi
	733	sl_tau = tau(ig,1)
	734	sl_alb = albedo(ig,l)
	735	sl_psi = psi_sl(ig)
	736	sl_fl0 = fluxsurf_sw(ig,l)
	737	sl_di0 = 0.
	738	if (mu0(ig) .gt. 0.) then
	739	sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
	740	sl_di0 = sl_di0*1370./dist_sol/dist_sol
	741	sl_di0 = sl_di0/ztim1
	742	sl_di0 = fluxsurf_sw(ig,l)*sl_di0
	743	endif
	744	! sait-on jamais (a cause des arrondis)
	745	if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
	746	!!!!!!!!!!!!!!!!!!!!!!!!!!
	747	CALL meso_param_slope( mu0(ig), declin, sl_ra, sl_lat,
	748	& sl_tau, sl_alb,
	749	& sl_the, sl_psi, sl_di0, sl_fl0, sl_flu)
	750	!!!!!!!!!!!!!!!!!!!!!!!!!!
	751	fluxsurf_sw(ig,l) = sl_flu
	752	!! sl_ls = 180.*zls/pi
	753	!! sl_lct = ptime24. + 180.long(ig)/pi/15.
	754	!! sl_lat = 180.*lati(ig)/pi
	755	!! sl_tau = tau(ig,1)
	756	!! sl_alb = albedo(ig,l)
	757	!! sl_the = theta_sl(ig)
	758	!! sl_psi = psi_sl(ig)
	759	!! sl_fl0 = fluxsurf_sw(ig,l)
	760	!! CALL param_slope_full(sl_ls, sl_lct, sl_lat,
	761	!! & sl_tau, sl_alb,
	762	!! & sl_the, sl_psi, sl_fl0, sl_flu)
	763	ENDDO
	764	!!! compute correction on IR flux as well
	765	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
	766	fluxsurf_lw(ig)= fluxsurf_lw(ig)*sky
	767	ENDIF
	768	ENDDO
	769	ccccc
	770	ccccc PARAM SLOPE
	771	ccccc
	772	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	773	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	774
	775
	776	c CO2 near infrared absorption
	777	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	778	call zerophys(ngrid*nlayer,zdtnirco2)
	779	if (callnirco2) then
	780	call nirco2abs (ngrid,nlayer,pplay,dist_sol,
	781	. mu0,fract,declin, zdtnirco2)
	782	endif
	783
	784	c Radiative flux from the sky absorbed by the surface (W.m-2)
	785	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	786	DO ig=1,ngrid
	787	fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
	788	$ +fluxsurf_sw(ig,1)*(1.-albedo(ig,1))
	789	$ +fluxsurf_sw(ig,2)*(1.-albedo(ig,2))
	790
	791	!print*,'RAD ', fluxrad_sky(ig)
	792	!print,'LW ', emis(ig)fluxsurf_lw(ig)
	793	!print,'SW ', fluxsurf_sw(ig,2)(1.-albedo(ig,2))
	794
	795	ENDDO
	796
	797
	798	c Net atmospheric radiative heating rate (K.s-1)
	799	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	800	IF(callnlte) THEN
	801	CALL blendrad(ngrid, nlayer, pplay,
	802	& zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
	803	ELSE
	804	DO l=1,nlayer
	805	DO ig=1,ngrid
	806	dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
	807	& +zdtnirco2(ig,l)
	808	ENDDO
	809	ENDDO
	810	ENDIF
	811
	812
	813
	814	ENDIF ! of if(mod(icount-1,iradia).eq.0)
	815
	816	c Transformation of the radiative tendencies:
	817	c -------------------------------------------
	818
	819	c Net radiative surface flux (W.m-2)
	820	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
	821	c
	822	DO ig=1,ngrid
	823	zplanck(ig)=tsurf(ig)*tsurf(ig)
	824	zplanck(ig)=emis(ig)*
	825	$ stephanzplanck(ig)zplanck(ig)
	826	fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
	827	cccc
	828	cccc param slope
	829	cccc
	830	sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
	831	fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig)
	832	cccc
	833	cccc
	834	cccc
	835	ENDDO
	836
	837
	838	DO l=1,nlayer
	839	DO ig=1,ngrid
	840	pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
	841	ENDDO
	842	ENDDO
	843
	844	ENDIF ! of IF (callrad)
	845
	846	!!!
	847	!!! WRF WRF WRF commented for smaller executables
	848	!!!
	849	!c-----------------------------------------------------------------------
	850	!c 3. Gravity wave and subgrid scale topography drag :
	851	!c -------------------------------------------------
	852	!
	853	!
	854	! IF(calllott)THEN
	855	!
	856	! CALL calldrag_noro(ngrid,nlayer,ptimestep,
	857	! & pplay,pplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
	858	!
	859	! DO l=1,nlayer
	860	! DO ig=1,ngrid
	861	! pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
	862	! pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
	863	! pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
	864	! ENDDO
	865	! ENDDO
	866	! ENDIF
	867
	868	c-----------------------------------------------------------------------
	869	c 4. Vertical diffusion (turbulent mixing):
	870	c -----------------------------------------
	871	c
	872	IF (calldifv) THEN
	873
	874
	875	DO ig=1,ngrid
	876	zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
	877	!write (,), fluxrad(ig), fluxgrd(ig), zflubid(ig)
	878	ENDDO
	879
	880	CALL zerophys(ngrid*nlayer,zdum1)
	881	CALL zerophys(ngrid*nlayer,zdum2)
	882	do l=1,nlayer
	883	do ig=1,ngrid
	884	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
	885	enddo
	886	enddo
	887
	888	c Calling vdif (Martian version WITH CO2 condensation)
	889	CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
	890	$ ptimestep,capcal,lwrite,
	891	$ pplay,pplev,zzlay,zzlev,z0,
	892	$ pu,pv,zh,pq,tsurf,emis,qsurf,
	893	$ zdum1,zdum2,zdh,pdq,zflubid,
	894	$ zdudif,zdvdif,zdhdif,zdtsdif,q2,
	895	& zdqdif,zdqsdif)
	896
	897	DO ig=1,ngrid
	898	!! sensible heat flux in W/m2
	899	hfx(ig) = zflubid(ig)-capcal(ig)*zdtsdif(ig)
	900	!! u star in similarity theory in m/s
	901	ust(ig) = 0.4
	902	. * sqrt( pu(ig,1)pu(ig,1) + pv(ig,1)pv(ig,1) )
	903	. / log( 1.E+0 + zzlay(ig,1)/z0 )
	904	ENDDO
	905
	906	! write (,) 'PHYS HFX cp zdts', hfx(100), zflubid(100),
	907	! . capcal(100),
	908	! . zdtsdif(100)
	909	! write (,) 'PHYS UST', ust(100)
	910
	911	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	912	!!! LES LES
	913	IF (flag_LES) THEN
	914
	915	write (,) '************************************************'
	916	write (,) '** LES mode: the difv part is only used to'
	917	write (,) '** provide HFX and UST to the dynamics'
	918	write (,) '** NB: - dudif, dvdif, dhdif, dqdif are set to 0'
	919	write (,) '** - tsurf is updated'
	920	write (,) '************************************************'
	921
	922	!!!
	923	!!! WRF WRF LES LES : en fait le subgrid scale n'etait pas mis a zero !!
	924	!!!
	925	DO ig=1,ngrid
	926	! !! sensible heat flux in W/m2
	927	! hfx(ig) = zflubid(ig)-capcal(ig)*zdtsdif(ig)
	928	! !! u star in similarity theory in m/s
	929	! ust(ig) = 0.4
	930	! . * sqrt( pu(ig,1)pu(ig,1) + pv(ig,1)pv(ig,1) )
	931	! . / log( 1.E+0 + zzlay(ig,1)/z0 )
	932	!
	933	DO l=1,nlayer
	934	zdvdif(ig,l) = 0.
	935	zdudif(ig,l) = 0.
	936	zdhdif(ig,l) = 0.
	937	DO iq=1, nq
	938	zdqdif(ig,l,iq) = 0.
	939	zdqsdif(ig,iq) = 0. !! sortir de la boucle
	940	ENDDO
	941	ENDDO
	942	!
	943	ENDDO
	944	!write (,) 'RAD ',fluxrad(igout)
	945	!write (,) 'GRD ',fluxgrd(igout)
	946	!write (,) 'dTs/dt ',capcal(igout)*zdtsurf(igout)
	947	!write (,) 'HFX ', hfx(igout)
	948	!write (,) 'UST ', ust(igout)
	949	ENDIF
	950	!!! LES LES
	951	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	952
	953	DO l=1,nlayer
	954	DO ig=1,ngrid
	955	pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
	956	pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
	957	pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
	958
	959	zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
	960
	961	ENDDO
	962	ENDDO
	963
	964	DO ig=1,ngrid
	965	zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
	966	ENDDO
	967
	968	if (tracer) then
	969	DO iq=1, nq
	970	DO l=1,nlayer
	971	DO ig=1,ngrid
	972	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq)
	973	ENDDO
	974	ENDDO
	975	ENDDO
	976	DO iq=1, nq
	977	DO ig=1,ngrid
	978	dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
	979	ENDDO
	980	ENDDO
	981	end if ! of if (tracer)
	982
	983	ELSE
	984	DO ig=1,ngrid
	985	zdtsurf(ig)=zdtsurf(ig)+
	986	s (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
	987	ENDDO
	988	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	989	IF (flag_LES) THEN
	990	write(,) 'LES mode !'
	991	write(,) 'Please set calldifv to T in callphys.def'
	992	STOP
	993	ENDIF
	994	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	995	ENDIF ! of IF (calldifv)
	996
	997
	998	c-----------------------------------------------------------------------
	999	c 5. Dry convective adjustment:
	1000	c -----------------------------
	1001
	1002	IF(calladj) THEN
	1003
	1004	DO l=1,nlayer
	1005	DO ig=1,ngrid
	1006	zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
	1007	ENDDO
	1008	ENDDO
	1009	CALL zerophys(ngrid*nlayer,zduadj)
	1010	CALL zerophys(ngrid*nlayer,zdvadj)
	1011	CALL zerophys(ngrid*nlayer,zdhadj)
	1012
	1013	CALL convadj(ngrid,nlayer,nq,ptimestep,
	1014	$ pplay,pplev,zpopsk,
	1015	$ pu,pv,zh,pq,
	1016	$ pdu,pdv,zdh,pdq,
	1017	$ zduadj,zdvadj,zdhadj,
	1018	$ zdqadj)
	1019
	1020	DO l=1,nlayer
	1021	DO ig=1,ngrid
	1022	pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
	1023	pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
	1024	pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
	1025
	1026	zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
	1027	ENDDO
	1028	ENDDO
	1029
	1030	if(tracer) then
	1031	DO iq=1, nq
	1032	DO l=1,nlayer
	1033	DO ig=1,ngrid
	1034	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq)
	1035	ENDDO
	1036	ENDDO
	1037	ENDDO
	1038	end if
	1039	ENDIF ! of IF(calladj)
	1040
	1041	c-----------------------------------------------------------------------
	1042	c 6. Carbon dioxide condensation-sublimation:
	1043	c -------------------------------------------
	1044
	1045	IF (callcond) THEN
[86]	1046	CALL newcondens(ngrid,nlayer,nq,ptimestep,
[42]	1047	$ capcal,pplay,pplev,tsurf,pt,
	1048	$ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
	1049	$ co2ice,albedo,emis,
	1050	$ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
	1051	$ fluxsurf_sw,zls)
	1052
	1053	DO l=1,nlayer
	1054	DO ig=1,ngrid
	1055	pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
	1056	pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
	1057	pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
	1058	ENDDO
	1059	ENDDO
	1060	DO ig=1,ngrid
	1061	zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
	1062	!!!**WRF: newphys: ici la pression n'est plus mise a jour ds le GCM
	1063	!!!**WRF: mais il faut retablir ca dans le cas du mesoscale ?
	1064	!!!**WRF: ...non probablement OK
	1065	! ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
	1066	ENDDO
	1067
	1068	IF (tracer) THEN
	1069	DO iq=1, nq
	1070	DO l=1,nlayer
	1071	DO ig=1,ngrid
	1072	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq)
	1073	ENDDO
	1074	ENDDO
	1075	ENDDO
	1076	ENDIF ! of IF (tracer)
	1077
	1078	ENDIF ! of IF (callcond)
	1079
	1080	c-----------------------------------------------------------------------
	1081	c 7. Specific parameterizations for tracers
	1082	c: -----------------------------------------
	1083
	1084	if (tracer) then
	1085
	1086	c 7a. Water and ice
	1087	c ---------------
	1088
	1089	c ---------------------------------------
	1090	c Water ice condensation in the atmosphere
	1091	c ----------------------------------------
	1092	IF (water) THEN
	1093
	1094	c **WRF: new arguments here rnuclei,rice,nuice
	1095	c plus no more iqmin +igcm_h2o_vap replaces iq, what are the consequences?
	1096	c checks needed when tracers simulations
	1097
	1098	call watercloud(ngrid,nlayer,ptimestep,
	1099	& pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
	1100	& pq,pdq,zdqcloud,zdqscloud,zdtcloud,
	1101	& nq,naerkind,tau,
	1102	& ccn,rdust,rice,nuice)
	1103	if (activice) then
	1104	c Temperature variation due to latent heat release
	1105	DO l=1,nlayer
	1106	DO ig=1,ngrid
	1107	pdt(ig,l)=pdt(ig,l)+zdtcloud(ig,l)
	1108	ENDDO
	1109	ENDDO
	1110	endif
	1111
	1112	! increment water vapour and ice atmospheric tracers tendencies
	1113	IF (water) THEN
	1114	DO l=1,nlayer
	1115	DO ig=1,ngrid
	1116	pdq(ig,l,igcm_h2o_vap)=pdq(ig,l,igcm_h2o_vap)+
	1117	& zdqcloud(ig,l,igcm_h2o_vap)
	1118	pdq(ig,l,igcm_h2o_ice)=pdq(ig,l,igcm_h2o_ice)+
	1119	& zdqcloud(ig,l,igcm_h2o_ice)
	1120	ENDDO
	1121	ENDDO
	1122	ENDIF ! of IF (water) THEN
	1123	! Increment water ice surface tracer tendency
	1124	DO ig=1,ngrid
	1125	dqsurf(ig,igcm_h2o_ice)=dqsurf(ig,igcm_h2o_ice)+
	1126	& zdqscloud(ig,igcm_h2o_ice)
	1127	ENDDO
	1128
	1129	END IF ! of IF (water)
	1130
	1131	c 7b. Chemical species
	1132	c ------------------
	1133
	1134	!!!
	1135	!!! WRF WRF WRF commented for smaller executables
	1136	!!!
	1137	!c --------------
	1138	!c photochemistry :
	1139	!c --------------
	1140	! IF (photochem .or. thermochem) then
	1141	! call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
	1142	! & zzlay,zday,pq,pdq,rice,
	1143	! & zdqchim,zdqschim,zdqcloud,zdqscloud)
	1144	!!NB: Photochemistry includes condensation of H2O2
	1145	!
	1146	! ! increment values of tracers:
	1147	! DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
	1148	! ! tracers is zero anyways
	1149	! DO l=1,nlayer
	1150	! DO ig=1,ngrid
	1151	! pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
	1152	! ENDDO
	1153	! ENDDO
	1154	! ENDDO ! of DO iq=1,nq
	1155	! ! add condensation tendency for H2O2
	1156	! if (igcm_h2o2.ne.0) then
	1157	! DO l=1,nlayer
	1158	! DO ig=1,ngrid
	1159	! pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
	1160	! & +zdqcloud(ig,l,igcm_h2o2)
	1161	! ENDDO
	1162	! ENDDO
	1163	! endif
	1164	!
	1165	! ! increment surface values of tracers:
	1166	! DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
	1167	! ! tracers is zero anyways
	1168	! DO ig=1,ngrid
	1169	! dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
	1170	! ENDDO
	1171	! ENDDO ! of DO iq=1,nq
	1172	! ! add condensation tendency for H2O2
	1173	! if (igcm_h2o2.ne.0) then
	1174	! DO ig=1,ngrid
	1175	! dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
	1176	! & +zdqscloud(ig,igcm_h2o2)
	1177	! ENDDO
	1178	! endif
	1179	!
	1180	! END IF ! of IF (photochem.or.thermochem)
	1181
	1182	c 7c. Aerosol particles
	1183	c -------------------
	1184
	1185	c ----------
	1186	c Dust devil :
	1187	c ----------
	1188	IF(callddevil) then
	1189	call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2,
	1190	& zdqdev,zdqsdev)
	1191
	1192	if (dustbin.ge.1) then
	1193	do iq=1,nq
	1194	DO l=1,nlayer
	1195	DO ig=1,ngrid
	1196	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
	1197	ENDDO
	1198	ENDDO
	1199	enddo
	1200	do iq=1,nq
	1201	DO ig=1,ngrid
	1202	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
	1203	ENDDO
	1204	enddo
	1205	endif ! of if (dustbin.ge.1)
	1206
	1207	END IF ! of IF (callddevil)
	1208
	1209	c -------------
	1210	c Sedimentation : acts also on water ice
	1211	c -------------
	1212	IF (sedimentation) THEN
	1213	!call zerophys(ngridnlayernq, zdqsed)
	1214	zdqsed(1:ngrid,1:nlayer,1:nq)=0
	1215	!call zerophys(ngrid*nq, zdqssed)
	1216	zdqssed(1:ngrid,1:nq)=0
	1217
	1218	c
	1219	c **WRF: new arguments rnuclei, rice, need checks
	1220	c
	1221	call callsedim(ngrid,nlayer, ptimestep,
	1222	& pplev,zzlev, pt, rdust, rice,
	1223	& pq, pdq, zdqsed, zdqssed,nq)
	1224
	1225	DO iq=1, nq
	1226	DO l=1,nlayer
	1227	DO ig=1,ngrid
	1228	pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
	1229	ENDDO
	1230	ENDDO
	1231	ENDDO
	1232	DO iq=1, nq
	1233	DO ig=1,ngrid
	1234	dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
	1235	ENDDO
	1236	ENDDO
	1237	END IF ! of IF (sedimentation)
	1238
	1239	c 7d. Updates
	1240	c ---------
	1241
	1242	DO iq=1, nq
	1243	DO ig=1,ngrid
	1244
	1245	c ---------------------------------
	1246	c Updating tracer budget on surface
	1247	c ---------------------------------
	1248	qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
	1249
	1250	ENDDO ! (ig)
	1251	ENDDO ! (iq)
	1252
	1253	endif ! of if (tracer)
	1254
	1255	!!!
	1256	!!! WRF WRF WRF commented for smaller executables
	1257	!!!
	1258	!c-----------------------------------------------------------------------
	1259	!c 8. THERMOSPHERE CALCULATION
	1260	!c-----------------------------------------------------------------------
	1261	!
	1262	! if (callthermos) then
	1263	! call thermosphere(pplev,pplay,dist_sol,
	1264	! $ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
	1265	! & pt,pq,pu,pv,pdt,pdq,
	1266	! $ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
	1267	!
	1268	! DO l=1,nlayer
	1269	! DO ig=1,ngrid
	1270	! dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
	1271	! pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
	1272	! & +zdteuv(ig,l)
	1273	! pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
	1274	! pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
	1275	! DO iq=1, nq
	1276	! pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
	1277	! ENDDO
	1278	! ENDDO
	1279	! ENDDO
	1280	!
	1281	! endif ! of if (callthermos)
	1282
	1283	c-----------------------------------------------------------------------
	1284	c 9. Surface and sub-surface soil temperature
	1285	c-----------------------------------------------------------------------
	1286	c
	1287	c
	1288	c 9.1 Increment Surface temperature:
	1289	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1290
	1291	DO ig=1,ngrid
	1292	tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig)
	1293	ENDDO
	1294
	1295	ccc
	1296	ccc **WRF very specific to GCM
	1297	ccc
	1298	c Prescribe a cold trap at south pole (except at high obliquity !!)
	1299	c Temperature at the surface is set there to be the temperature
	1300	c corresponding to equilibrium temperature between phases of CO2
	1301
	1302	IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
[73]	1303	! if (caps.and.(obliquit.lt.27.)) then
	1304	! ! NB: Updated surface pressure, at grid point 'ngrid', is
	1305	! ! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
	1306	! tsurf(ngrid)=1./(1./136.27-r/5.9e+5alog(0.0095
	1307	! & (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
	1308	! endif
	1309	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1310	!!!!! note WRF MESOSCALE AYMERIC -- mot cle "caps"
	1311	!!!!! watercaptag n'est plus utilise que dans vdifc
	1312	!!!!! ... pour que la sublimation ne soit pas stoppee
	1313	!!!!! ... dans la calotte permanente nord si qsurf=0
	1314	!!!!! on desire garder cet effet regle par caps=T
	1315	!!!!! on a donc commente "if (caps.and.(obliquit.lt.27.))" ci-dessus
	1316	!!!!! --- remplacer ces lignes par qqch de plus approprie
	1317	!!!!! si on s attaque a la calotte polaire sud
	1318	!!!!! pas d'autre occurrence majeure du mot-cle "caps"
	1319	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1320
[42]	1321	c -------------------------------------------------------------
	1322	c Change of surface albedo (set to 0.4) in case of ground frost
	1323	c everywhere except on the north permanent cap and in regions
	1324	c covered by dry ice.
	1325	c ALWAYS PLACE these lines after newcondens !!!
	1326	c -------------------------------------------------------------
	1327	c
	1328	c **WRF : OK avec le mesoscale, pas d'indices bizarres au pole
	1329	c
	1330	do ig=1,ngrid
	1331	if ((co2ice(ig).eq.0).and.
	1332	& (qsurf(ig,igcm_h2o_ice).gt.0.005)) then
	1333	albedo(ig,1) = alb_surfice
	1334	albedo(ig,2) = alb_surfice
	1335	endif
	1336	enddo ! of do ig=1,ngrid
	1337	ENDIF ! of IF (tracer.AND.water.AND.(ngridmx.NE.1))
	1338
	1339	c
	1340	c 9.2 Compute soil temperatures and subsurface heat flux:
	1341	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1342	IF (callsoil) THEN
	1343	CALL soil(ngrid,nsoilmx,.false.,inertiedat,
	1344	& ptimestep,tsurf,tsoil,capcal,fluxgrd)
	1345	ENDIF
	1346
	1347	c-----------------------------------------------------------------------
	1348	c 10. Write output files
	1349	c ----------------------
	1350
	1351	c -------------------------------
	1352	c Dynamical fields incrementation
	1353	c -------------------------------
	1354	c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
	1355	! temperature, zonal and meridional wind
	1356	DO l=1,nlayer
	1357	DO ig=1,ngrid
	1358	zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep
	1359	zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep
	1360	zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep
	1361	ENDDO
	1362	ENDDO
	1363
	1364	! tracers
	1365	DO iq=1, nq
	1366	DO l=1,nlayer
	1367	DO ig=1,ngrid
	1368	zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
	1369	ENDDO
	1370	ENDDO
	1371	ENDDO
	1372
	1373	! surface pressure
	1374	DO ig=1,ngrid
	1375	ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
	1376	ENDDO
	1377
	1378	! pressure
	1379	DO l=1,nlayer
	1380	DO ig=1,ngrid
	1381	zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
	1382	zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
	1383	ENDDO
	1384	ENDDO
	1385
	1386	! Density
	1387	DO l=1,nlayer
	1388	DO ig=1,ngrid
	1389	rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
	1390	ENDDO
	1391	ENDDO
	1392
	1393	c Compute surface stress : (NB: z0 is a common in planete.h)
	1394	c DO ig=1,ngrid
	1395	c cd = (0.4/log(zzlay(ig,1)/z0))**2
	1396	c zstress(ig) = rho(ig,1)cd(zu(ig,1)2 + zv(ig,1)2)
	1397	c ENDDO
	1398
	1399	c Sum of fluxes in solar spectral bands (for output only)
	1400	DO ig=1,ngrid
	1401	fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
	1402	fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
	1403	ENDDO
	1404	c ***** TEST ****************************************************
	1405	ztim1 = 999
	1406	DO l=1,nlayer
	1407	DO ig=1,ngrid
	1408	if (pt(ig,l).lt.ztim1) then
	1409	ztim1 = pt(ig,l)
	1410	igmin = ig
	1411	lmin = l
	1412	end if
	1413	ENDDO
	1414	ENDDO
	1415	if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
	1416	write(,) 'PHYSIQ: stability WARNING :'
	1417	write(,) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
	1418	& 'ig l =', igmin, lmin
	1419	end if
	1420	c *******************************************************************
	1421
	1422	c ---------------------
	1423	c Outputs to the screen
	1424	c ---------------------
	1425
	1426	IF (lwrite) THEN
	1427	PRINT*,'Global diagnostics for the physics'
	1428	PRINT,'Variables and their increments x and dx/dt dt'
	1429	WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
	1430	WRITE(*,'(2f10.5,2f15.5)')
	1431	s tsurf(igout),zdtsurf(igout)*ptimestep,
	1432	s pplev(igout,1),pdpsrf(igout)*ptimestep
	1433	WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
	1434	WRITE(*,'(i4,6f10.5)') (l,
	1435	s pu(igout,l),pdu(igout,l)*ptimestep,
	1436	s pv(igout,l),pdv(igout,l)*ptimestep,
	1437	s pt(igout,l),pdt(igout,l)*ptimestep,
	1438	s l=1,nlayer)
	1439	ENDIF ! of IF (lwrite)
	1440
	1441	IF (ngrid.NE.1) THEN
	1442	print,'Ls =',zls180./pi,
	1443	& ' tauref(700 Pa,lat=0) =',tauref(ngrid/2)!,
	1444	! & ' tau(Viking1) =',tau(ig_vl1,1)
	1445
	1446
	1447	c -------------------------------------------------------------------
	1448	c Writing NetCDF file "RESTARTFI" at the end of the run
	1449	c -------------------------------------------------------------------
	1450	c Note: 'restartfi' is stored just before dynamics are stored
	1451	c in 'restart'. Between now and the writting of 'restart',
	1452	c there will have been the itau=itau+1 instruction and
	1453	c a reset of 'time' (lastacll = .true. when itau+1= itaufin)
	1454	c thus we store for time=time+dtvr
	1455
	1456
	1457	!!!
	1458	!!! WRF WRF WRF WRF
	1459	!!!
	1460	! IF(lastcall) THEN
	1461	! ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec)
	1462	! write(,)'PHYSIQ: for physdem ztime_fin =',ztime_fin
	1463	! call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
	1464	! . ptimestep,pday,
	1465	! . ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf,
	1466	! . area,albedodat,inertiedat,zmea,zstd,zsig,
	1467	! . zgam,zthe)
	1468	! ENDIF
	1469
	1470
	1471
	1472	c -------------------------------------------------------------------
	1473	c Calculation of diagnostic variables written in both stats and
	1474	c diagfi files
	1475	c -------------------------------------------------------------------
	1476
	1477	if (tracer) then
	1478	if (water) then
	1479
	1480	!!
	1481	!!***WRF: ok, des nouveaux trucs cools de la nouvelle physique
	1482	!!
	1483	call zerophys(ngrid,mtot)
	1484	call zerophys(ngrid,icetot)
	1485	call zerophys(ngrid,rave)
	1486	call zerophys(ngrid,tauTES)
	1487	do ig=1,ngrid
	1488	do l=1,nlayermx
	1489	mtot(ig) = mtot(ig) +
	1490	& zq(ig,l,igcm_h2o_vap) *
	1491	& (pplev(ig,l) - pplev(ig,l+1)) / g
	1492	icetot(ig) = icetot(ig) +
	1493	& zq(ig,l,igcm_h2o_ice) *
	1494	& (pplev(ig,l) - pplev(ig,l+1)) / g
	1495	rave(ig) = rave(ig) +
	1496	& zq(ig,l,igcm_h2o_ice) *
	1497	& (pplev(ig,l) - pplev(ig,l+1)) / g *
	1498	& rice(ig,l) * (1.+nuice_ref)
	1499	c Computing abs optical depth at 825 cm-1 in each
	1500	c layer to simulate NEW TES retrieval
	1501	Qabsice = min(
	1502	& max(0.4e6rice(ig,l)(1.+nuice_ref)-0.05 ,0.),1.2
	1503	& )
	1504	opTES(ig,l)= 0.75 * Qabsice *
	1505	& zq(ig,l,igcm_h2o_ice) *
	1506	& (pplev(ig,l) - pplev(ig,l+1)) / g
	1507	& / (rho_ice * rice(ig,l) * (1.+nuice_ref))
	1508	tauTES(ig)=tauTES(ig)+ opTES(ig,l)
	1509	enddo
	1510	rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
	1511	if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
	1512	enddo
	1513
	1514	endif ! of if (water)
	1515	endif ! of if (tracer)
	1516
	1517	c -----------------------------------------------------------------
	1518	c WSTATS: Saving statistics
	1519	c -----------------------------------------------------------------
	1520	c ("stats" stores and accumulates 8 key variables in file "stats.nc"
	1521	c which can later be used to make the statistic files of the run:
	1522	c "stats") only possible in 3D runs !
	1523
	1524
	1525	IF (callstats) THEN
	1526
	1527	write(,) 'callstats'
	1528
	1529	! call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
	1530	! call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
	1531	! call wstats(ngrid,"co2ice","CO2 ice cover",
	1532	! & "kg.m-2",2,co2ice)
	1533	! call wstats(ngrid,"fluxsurf_lw",
	1534	! & "Thermal IR radiative flux to surface","W.m-2",2,
	1535	! & fluxsurf_lw)
	1536	! call wstats(ngrid,"fluxsurf_sw",
	1537	! & "Solar radiative flux to surface","W.m-2",2,
	1538	! & fluxsurf_sw_tot)
	1539	! call wstats(ngrid,"fluxtop_lw",
	1540	! & "Thermal IR radiative flux to space","W.m-2",2,
	1541	! & fluxtop_lw)
	1542	! call wstats(ngrid,"fluxtop_sw",
	1543	! & "Solar radiative flux to space","W.m-2",2,
	1544	! & fluxtop_sw_tot)
	1545	! call wstats(ngrid,"taudustvis",
	1546	! & "Dust optical depth"," ",2,tau(1,1))
	1547	! call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
	1548	! call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
	1549	! call wstats(ngrid,"v","Meridional (North-South) wind",
	1550	! & "m.s-1",3,zv)
	1551	!c call wstats(ngrid,"w","Vertical (down-up) wind",
	1552	!c & "m.s-1",3,pw)
	1553	! call wstats(ngrid,"rho","Atmospheric density","none",3,rho)
	1554	!c call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
	1555	!c call wstats(ngrid,"q2",
	1556	!c & "Boundary layer eddy kinetic energy",
	1557	!c & "m2.s-2",3,q2)
	1558	!c call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
	1559	!c & emis)
	1560	!c call wstats(ngrid,"ssurf","Surface stress","N.m-2",
	1561	!c & 2,zstress)
	1562	!
	1563	! if (tracer) then
	1564	! if (water) then
	1565	! vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
	1566	! & *mugaz/mmol(igcm_h2o_vap)
	1567	! call wstats(ngrid,"vmr_h2ovapor",
	1568	! & "H2O vapor volume mixing ratio","mol/mol",
	1569	! & 3,vmr)
	1570	! vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
	1571	! & *mugaz/mmol(igcm_h2o_ice)
	1572	! call wstats(ngrid,"vmr_h2oice",
	1573	! & "H2O ice volume mixing ratio","mol/mol",
	1574	! & 3,vmr)
	1575	!
	1576	! call wstats(ngrid,"mtot",
	1577	! & "total mass of water vapor","kg/m2",
	1578	! & 2,mtot)
	1579	! call wstats(ngrid,"icetot",
	1580	! & "total mass of water ice","kg/m2",
	1581	! & 2,icetot)
	1582	!c If activice is true, tauTES is computed in aeropacity.F;
	1583	! if (.not.activice) then
	1584	! call wstats(ngrid,"tauTES",
	1585	! & "tau abs 825 cm-1","",
	1586	! & 2,tauTES)
	1587	! endif ! of if (activice)
	1588	!
	1589	! endif ! of if (water)
	1590	!
	1591	! if (thermochem.or.photochem) then
	1592	! do iq=1,nq
	1593	! if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or.
	1594	! . (noms(iq).eq."co").or.(noms(iq).eq."n2").or.
	1595	! . (noms(iq).eq."h2").or.
	1596	! . (noms(iq).eq."o3")) then
	1597	! do l=1,nlayer
	1598	! do ig=1,ngrid
	1599	! vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
	1600	! end do
	1601	! end do
	1602	! call wstats(ngrid,"vmr_"//trim(noms(iq)),
	1603	! . "Volume mixing ratio","mol/mol",3,vmr)
	1604	! endif
	1605	! enddo
	1606	! endif ! of if (thermochem.or.photochem)
	1607	!
	1608	! endif ! of if (tracer)
	1609	!
	1610	! IF(lastcall) THEN
	1611	! write (,) "Writing stats..."
	1612	! call mkstats(ierr)
	1613	! ENDIF
	1614
	1615	ENDIF !if callstats
	1616
	1617	c (Store EOF for Mars Climate database software)
	1618	IF (calleofdump) THEN
	1619	CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
	1620	ENDIF
	1621
	1622	ccc************** WRF OUTPUT ************************
	1623	ccc************** WRF OUTPUT ************************
	1624	ccc************** WRF OUTPUT ************************
[81]	1625	!do ig=1,ngrid
	1626	! wtsurf(ig) = tsurf(ig) !! surface temperature
	1627	! wco2ice(ig) = co2ice(ig) !! co2 ice
	1628	!
	1629	! !!! specific to WRF WRF WRF
	1630	! !!! just to output water ice on surface
	1631	! !!! uncomment the Registry entry
	1632	! IF (igcm_h2o_ice .ne. 0) qsurfice(ig) = qsurf(ig,igcm_h2o_ice)
	1633	!
	1634	! !!! "VMR_ICE" "VOL. MIXING RATIO ICE" "ppm"
	1635	! IF (igcm_h2o_ice .ne. 0) THEN
	1636	! vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)*mugaz/mmol(igcm_h2o_ice)
	1637	! ENDIF
	1638	!
	1639	!enddo
	1640	wtsurf(1:ngrid) = tsurf(1:ngrid) !! surface temperature
	1641	wco2ice(1:ngrid) = co2ice(1:ngrid) !! co2 ice
[82]	1642	mtot(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
[81]	1643	IF (igcm_h2o_ice .ne. 0) THEN
	1644	qsurfice(1:ngrid) = qsurf(1:ngrid,igcm_h2o_ice)
	1645	vmr=1.e6 * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
	1646	. * mugaz / mmol(igcm_h2o_ice)
	1647	ENDIF
[73]	1648
[42]	1649	c
	1650	c THIS INCLUDE IS AUTOMATICALLY GENERATED FROM REGISTRY
	1651	c
	1652	#include "fill_save.inc"
	1653	c
	1654	ccc************** WRF OUTPUT ************************
	1655	ccc************** WRF OUTPUT ************************
	1656	ccc************** WRF OUTPUT ************************
	1657
	1658
	1659	cc-----------------------------------
	1660	cc you can still use meso_WRITEDIAGFI (e.g. for debugging purpose),
	1661	cc though this is not the default strategy now
	1662	cc-----------------------------------
	1663	cc please use cudt in namelist.input to set frequency of outputs
	1664	cc-----------------------------------
	1665	cc BEWARE: if at least one call to meso_WRITEDIAGFI is performed,
	1666	cc cudt cannot be 0 - otherwise you'll get a "Floating exception"
	1667	cc-----------------------------------
	1668	! call meso_WRITEDIAGFI(ngrid,"tauref",
	1669	! . "tauref","W.m-2",2,
	1670	! . tauref)
	1671	! call meso_WRITEDIAGFI(ngrid,"dtrad",
	1672	! . "dtrad","W.m-2",2,
	1673	! . dtrad)
	1674	c call meso_WRITEDIAGFI(ngrid,"tsurf",
	1675	c . "tsurf","K",2,
	1676	c . tsurf)
	1677	c
	1678	! call meso_WRITEDIAGFI(ngrid,"zt",
	1679	! . "zt","W.m-2",3,
	1680	! . zt)
	1681	! call meso_WRITEDIAGFI(ngrid,"zdtlw",
	1682	! . "zdtlw","W.m-2",2,
	1683	! . zdtlw)
	1684	! call meso_WRITEDIAGFI(ngrid,"zdtsw",
	1685	! . "zdtsw","W.m-2",2,
	1686	! . zdtsw)
	1687
	1688
	1689	!!
	1690	!! ***WRF: everything below is kept for reference
	1691	!!
	1692	!
	1693	!c ==========================================================
	1694	!c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
	1695	!c any variable for diagnostic (output with period
	1696	!c "ecritphy", set in "run.def")
	1697	!c ==========================================================
	1698	!c WRITEDIAGFI can ALSO be called from any other subroutines
	1699	!c for any variables !!
	1700	! call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
	1701	! & emis)
	1702	! call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
	1703	! & tsurf)
	1704	! call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
	1705	! call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
	1706	! & co2ice)
	1707	!c call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
	1708	!c & "K",2,zt(1,7))
	1709	! call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
	1710	! & fluxsurf_lw)
	1711	! call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
	1712	! & fluxsurf_sw_tot)
	1713	! call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
	1714	! & fluxtop_lw)
	1715	! call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
	1716	! & fluxtop_sw_tot)
	1717	! call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
	1718	!c call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
	1719	!c call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
	1720	!c call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
	1721	! call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
	1722	!c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
	1723	!c call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
	1724	!c call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay)
	1725	!c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
	1726	!c & zstress)
	1727	!
	1728	!c ----------------------------------------------------------
	1729	!c Outputs of the CO2 cycle
	1730	!c ----------------------------------------------------------
	1731	!
	1732	! if (tracer.and.(igcm_co2.ne.0)) then
	1733	!! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
	1734	!! & "kg/kg",2,zq(1,1,igcm_co2))
	1735	!! call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
	1736	!! & "kg/kg",3,zq(1,1,igcm_co2))
	1737	!
	1738	! ! Compute co2 column
	1739	! call zerophys(ngrid,co2col)
	1740	! do l=1,nlayermx
	1741	! do ig=1,ngrid
	1742	! co2col(ig)=co2col(ig)+
	1743	! & zq(ig,l,igcm_co2)*(pplev(ig,l)-pplev(ig,l+1))/g
	1744	! enddo
	1745	! enddo
	1746	! call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
	1747	! & co2col)
	1748	! endif ! of if (tracer.and.(igcm_co2.ne.0))
	1749	!
	1750	!c ----------------------------------------------------------
	1751	!c Outputs of the water cycle
	1752	!c ----------------------------------------------------------
	1753	! if (tracer) then
	1754	! if (water) then
	1755	!
	1756	! CALL WRITEDIAGFI(ngridmx,'mtot',
	1757	! & 'total mass of water vapor',
	1758	! & 'kg/m2',2,mtot)
	1759	! CALL WRITEDIAGFI(ngridmx,'icetot',
	1760	! & 'total mass of water ice',
	1761	! & 'kg/m2',2,icetot)
	1762	!c If activice is true, tauTES is computed in aeropacity.F;
	1763	! if (.not.activice) then
	1764	! CALL WRITEDIAGFI(ngridmx,'tauTES',
	1765	! & 'tau abs 825 cm-1',
	1766	! & '',2,tauTES)
	1767	! endif
	1768	!
	1769	! call WRITEDIAGFI(ngridmx,'h2o_ice_s',
	1770	! & 'surface h2o_ice',
	1771	! & 'kg.m-2',2,qsurf(1,igcm_h2o_ice))
	1772	!
	1773	! if (activice) then
	1774	!c call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
	1775	!c & 'w.m-2',3,zdtsw)
	1776	!c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
	1777	!c & 'w.m-2',3,zdtlw)
	1778	! endif !(activice)
	1779	! endif !(water)
	1780	!
	1781	!
	1782	! if (water.and..not.photochem) then
	1783	! iq=nq
	1784	!c write(str2(1:2),'(i2.2)') iq
	1785	!c call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud',
	1786	!c & 'kg.m-2',2,zdqscloud(1,iq))
	1787	!c call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim',
	1788	!c & 'kg/kg',3,zdqchim(1,1,iq))
	1789	!c call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif',
	1790	!c & 'kg/kg',3,zdqdif(1,1,iq))
	1791	!c call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj',
	1792	!c & 'kg/kg',3,zdqadj(1,1,iq))
	1793	!c call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c',
	1794	!c & 'kg/kg',3,zdqc(1,1,iq))
	1795	! endif !(water.and..not.photochem)
	1796	! endif
	1797	!
	1798	!c ----------------------------------------------------------
	1799	!c Outputs of the dust cycle
	1800	!c ----------------------------------------------------------
	1801	!
	1802	! call WRITEDIAGFI(ngridmx,'taudustvis',
	1803	! & 'Dust optical depth',' ',2,tau(1,1))
	1804	!
	1805	! if (tracer.and.(dustbin.ne.0)) then
	1806	! call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
	1807	! if (doubleq) then
	1808	! call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
	1809	! & 'kg.m-2',2,qsurf(1,1))
	1810	! call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
	1811	! & 'N.m-2',2,qsurf(1,2))
	1812	! call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
	1813	! & 'kg.m-2.s-1',2,zdqsdev(1,1))
	1814	! call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
	1815	! & 'kg.m-2.s-1',2,zdqssed(1,1))
	1816	! do l=1,nlayer
	1817	! do ig=1, ngrid
	1818	! reff(ig,l)= ref_r0 *
	1819	! & (r3n_qpq(ig,l,1)/max(pq(ig,l,2),0.01))*(1./3.)
	1820	! reff(ig,l)=min(max(reff(ig,l),1.e-10),500.e-6)
	1821	! end do
	1822	! end do
	1823	! call WRITEDIAGFI(ngridmx,'reff','reff','m',3,reff)
	1824	! else
	1825	! do iq=1,dustbin
	1826	! write(str2(1:2),'(i2.2)') iq
	1827	! call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
	1828	! & 'kg/kg',3,zq(1,1,iq))
	1829	! call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf',
	1830	! & 'kg.m-2',2,qsurf(1,iq))
	1831	! end do
	1832	! endif ! (doubleq)
	1833	! end if ! (tracer.and.(dustbin.ne.0))
	1834	!
	1835	!c ----------------------------------------------------------
	1836	!c Output in netcdf file "diagsoil.nc" for subterranean
	1837	!c variables (output every "ecritphy", as for writediagfi)
	1838	!c ----------------------------------------------------------
	1839	!
	1840	! ! Write soil temperature
	1841	!! call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
	1842	!! & 3,tsoil)
	1843	! ! Write surface temperature
	1844	!! call writediagsoil(ngrid,"tsurf","Surface temperature","K",
	1845	!! & 2,tsurf)
	1846	!
	1847	!c ==========================================================
	1848	!c END OF WRITEDIAGFI
	1849	!c ==========================================================
	1850
	1851	ELSE ! if(ngrid.eq.1)
	1852
	1853	print,'Ls =',zls180./pi,
	1854	& ' tauref(700 Pa) =',tauref
	1855	c ----------------------------------------------------------------------
	1856	c Output in grads file "g1d" (ONLY when using testphys1d)
	1857	c (output at every X physical timestep)
	1858	c ----------------------------------------------------------------------
	1859	c
	1860	c CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
	1861	c CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
	1862	c CALL writeg1d(ngrid,1,ps,'ps','Pa')
	1863
	1864	c CALL writeg1d(ngrid,nlayer,zt,'T','K')
	1865	c CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
	1866	c CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
	1867	c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
	1868
	1869	!! or output in diagfi.nc (for testphys1d)
	1870	! call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps)
	1871	! call WRITEDIAGFI(ngridmx,'temp','Temperature',
	1872	! & 'K',1,zt)
	1873	!
	1874	! if(tracer) then
	1875	!c CALL writeg1d(ngrid,1,tau,'tau','SI')
	1876	! do iq=1,nq
	1877	!c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg')
	1878	! call WRITEDIAGFI(ngridmx,trim(noms(iq)),
	1879	! & trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
	1880	! end do
	1881	! end if
	1882	!
	1883	! zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g
	1884	!
	1885	! do l=2,nlayer-1
	1886	! tmean=zt(1,l)
	1887	! if(zt(1,l).ne.zt(1,l-1))
	1888	! & tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
	1889	! zlocal(l)= zlocal(l-1)
	1890	! & -log(pplay(1,l)/pplay(1,l-1))rnew(1,l)tmean/g
	1891	! enddo
	1892	! zlocal(nlayer)= zlocal(nlayer-1)-
	1893	! & log(pplay(1,nlayer)/pplay(1,nlayer-1))*
	1894	! & rnew(1,nlayer)*tmean/g
	1895
	1896	END IF ! if(ngrid.ne.1)
	1897
	1898	icount=icount+1
	1899	write(,) 'now, back to the dynamical core...'
[86]	1900	#endif
[42]	1901	RETURN
	1902	END

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format