Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

callradite_mod.F

Last change on this file was 3468, checked in by emillour, 4 weeks ago
Mars PCM: Remove obsolete/depreciated lwrite flag (which would trigger some very specific extra text outputs), in code and in reference callphys.def files. EM
File size: 24.0 KB

Rev	Line
[1711]	1	MODULE callradite_mod
	2
	3	IMPLICIT NONE
	4
	5	CONTAINS
	6
[38]	7	SUBROUTINE callradite(icount,ngrid,nlayer,nq,zday,ls,pq,albedo,
	8	$ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
[2685]	9	$ dtlw,dtsw,fluxsurf_lw,fluxsurf_dn_sw,fluxsurf_up_sw,
	10	$ fluxtop_lw,fluxtop_dn_sw,fluxtop_up_sw,
	11	& tau_pref_scenario,tau_pref_gcm,
[2643]	12	& tau,aerosol,dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef,
[2447]	13	$ taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,co2ice,
	14	$ rstormdust,rtopdust,totstormfract,clearatm,dsords,dsotop,
[2628]	15	$ nohmons,clearsky,totcloudfrac)
[38]	16
[1918]	17	use aeropacity_mod, only: aeropacity
[1969]	18	use updatereffrad_mod, only: updatereffrad
[1047]	19	use dimradmars_mod, only: ndomainsz, nflev, nsun, nir
[1246]	20	use dimradmars_mod, only: naerkind, name_iaer,
	21	& iaer_dust_conrath,iaer_dust_doubleq,
[2447]	22	& iaer_dust_submicron, iaer_h2o_ice, iaer_co2_ice,
[2678]	23	& iaer_stormdust_doubleq,iaer_topdust_doubleq,flux_1AU
[1047]	24	use yomlw_h, only: gcp, nlaylte
[1524]	25	use comcstfi_h, only: g,cpp
	26	use time_phylmdz_mod, only: daysec
[1983]	27	use lwmain_mod, only: lwmain
	28	use swmain_mod, only: swmain
[2409]	29	use dust_param_mod, only: doubleq, active, submicron
[1047]	30	IMPLICIT NONE
[38]	31	c=======================================================================
	32	c subject:
	33	c --------
	34	c Subroutine designed to call the main canonic
	35	c radiative transfer subroutine "lwmain" et "swmain"
	36	c to compute radiative heating and cooling rate and
	37	c radiative fluxes to the surface.
	38	c
	39	c These calculations are only valid on the part of the atmosphere
[2246]	40	c where Local Thermal Equilibrium (LTE) is verified. In practice
	41	c the calculations are only performed for the first "nlaylte"
[38]	42	c parameters (nlaylte is calculated by subroutine "nlthermeq"
[1047]	43	c and stored in module "yomlw_h").
[38]	44	c
	45	c The purpose of this subroutine is to:
	46	c 1) Make some initial calculation at first call
	47	c 2) Split the calculation in several sub-grid
	48	c ("sub-domain") to save memory and
	49	c be able run on a workstation at high resolution
[1047]	50	c The sub-grid size is defined in dimradmars_mod
[38]	51	c 3) Compute the 3D scattering parameters depending on the
	52	c size distribution of the different tracers (added by JBM)
	53	c 4) call "lwmain" and "swmain"
	54	c
	55	c
	56	c authors:
	57	c ------
	58	c Francois Forget / Christophe Hourdin / J.-B. Madeleine (2009)
	59	c
	60	c
	61	c 3D scattering scheme user's guide (J.-B. Madeleine)
	62	c ---------------------------------
	63	c
	64	c This routine has been modified to take into account 3D, time
	65	c dependent scattering properties of the aerosols.
	66	c---- The look-up tables that contain the scattering parameters
	67	c of a given tracer, for different sizes, are read by SUAER.F90.
	68	c The names of the corresponding ASCII files have to be set in
	69	c this subroutine (file_id variable), and files must be in the
[1918]	70	c directory specified in datafile_mod. Please make sure that the
[38]	71	c ASCII files are correctly written, and that the range
	72	c of particle sizes is consistent with what you would expect.
	73	c---- SUAER.F90 is in charge of reading the ASCII files and averaging
	74	c the scattering parameters in each GCM channel, using the three last
	75	c equations of Forget et al. 1998 (GRL 25, No.7, p.1105-1108).
	76	c---- These look-up tables, loaded during the firstcall, are then
	77	c constantly used by the subroutine "aeroptproperties.F" to compute,
	78	c online, the 3D scattering parameters, based on the size distribution
	79	c (reffrad and nueffrad) of the different tracers, in each grid box.
	80	c These 3D size distributions are loaded by the "updatereffrad.F"
	81	c subroutine. A log-normal distribution is then assumed in
	82	c "aeroptproperties.F", along with a Gauss-Legendre integration.
	83	c---- The optical depth at the visible reference wavelength (set in
	84	c SUAER.F90, after the file_id variable) is then computed by
	85	c the subroutine "aeropacity.F", by using the size and spatial
	86	c distribution of the corresponding tracer. This connection has to
	87	c be implemented in "aeropacity.F" when adding a new tracer. To do so,
	88	c one can use equation 2 of Forget et al. 1998 (Icarus 131, p.302-316).
	89	c---- The resulting variables "aerosol", "QVISsQREF3d", "omegaVIS3d" and
	90	c "gVIS3d" (same in the infrared) are finally used by lwmain.F and
	91	c swmain.F to solve the radiative transfer equation.
	92	c
	93	c changes:
	94	c -------
	95	c
	96	c > SRL 7/2000
	97	c
	98	c This version has been modified to only calculate radiative tendencies
[1047]	99	c over layers 1..NFLEV (set in dimradmars_mod). Returns zero for higher
[38]	100	c layers, if any.
[1266]	101	c In other routines, nlayer -> nflev.
[38]	102	c Routines affected: lwflux, lwi, lwmain, lwxb, lwxd, lwxn.
	103	c
	104	c > J.-B. Madeleine 10W12
	105	c This version uses the variable's splitting, which can be usefull
	106	c when performing very high resolution simulation like LES.
	107	c
	108	c ----------
	109	c Here, solar band#1 is spectral interval between "long1vis" and "long2vis"
[1047]	110	c set in dimradmars_mod
[38]	111	c Here, solar band#2 is spectral interval between "long2vis" and "long3vis"
[1047]	112	c set in dimradmars_mod
[38]	113	c
	114	c input:
	115	c -----
	116	c icount counter of call to subroutine physic by gcm
	117	c ngrid number of gridpoint of horizontal grid
	118	c nlayer Number of layer
	119	c nq Number of tracer
	120	c ls Solar longitude (Ls) , radian
	121	c zday Date (time since Ls=0, in martian days)
	122	c pq(ngrid,nlayer,nq) Advected fields
	123	c
	124	c albedo (ngrid,2) hemispheric surface albedo
	125	c albedo (i,1) : mean albedo for solar band#1
	126	c (see below)
	127	c albedo (i,2) : mean albedo for solar band#2
	128	c (see below)
	129	c emis Thermal IR surface emissivity (no unit)
[1047]	130	c mu0(ngrid) cos of solar zenith angle
[38]	131	c (=1 when sun at zenith)
	132	c pplay(ngrid,nlayer) pressure (Pa) in the middle of each layer
	133	c pplev(ngrid,nlayer+1) pressure (Pa) at boundaries of each layer
	134	c pt(ngrid,nlayer) atmospheric temperature in each layer (K)
	135	c tsurf(ngrid) surface temperature (K)
[1047]	136	c fract(ngrid) day fraction of the time interval
[38]	137	c =1 during the full day ; =0 during the night
	138	c declin latitude of subsolar point
	139	c dist_sol sun-Mars distance (AU)
	140	c igout coordinate of analysed point for debugging
	141	c reffrad(ngrid,nlayer,naerkind) Aerosol effective radius
	142	c nueffrad(ngrid,nlayer,naerkind) Aerosol effective variance
	143
	144	c=======================================================================
	145	c
	146	c Declarations :
	147	c -------------
	148	c
[1974]	149	include "callkeys.h"
[38]	150
	151	c-----------------------------------------------------------------------
	152	c Input/Output
	153	c ------------
[1047]	154	INTEGER,INTENT(IN) :: icount
	155	INTEGER,INTENT(IN) :: ngrid,nlayer,nq
	156	INTEGER,INTENT(IN) :: igout
[38]	157
[1047]	158	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)
[1974]	159	REAL,INTENT(INOUT) :: tauscaling(ngrid) ! Conversion factor for
[358]	160	! qdust and Ndust
[2634]	161	REAL,INTENT(INOUT) :: dust_rad_adjust(ngrid) ! Radiative adjustment
[2417]	162	! factor for dust
[2643]	163	REAL,INTENT(INOUT) :: IRtoVIScoef(ngrid) ! conversion coefficient to apply on
	164	! scenario absorption IR (9.3um) CDOD
	165	! = tau_pref_gcm_VIS / tau_pref_gcm_IR
[1047]	166	REAL,INTENT(IN) :: albedo(ngrid,2),emis(ngrid)
	167	REAL,INTENT(IN) :: ls,zday
[38]	168
[1047]	169	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
	170	REAL,INTENT(IN) :: pt(ngrid,nlayer)
	171	REAL,INTENT(IN) :: tsurf(ngrid)
	172	REAL,INTENT(IN) :: dist_sol,mu0(ngrid),fract(ngrid)
[2685]	173	REAL,INTENT(OUT) :: dtlw(ngrid,nlayer) ! longwave (IR) heating rate (K/s)
	174	REAL,INTENT(OUT) :: dtsw(ngrid,nlayer) ! shortwave (Solar) heating rate (K/s)
	175	REAL,INTENT(OUT) :: fluxsurf_lw(ngrid) ! total LW (thermal IR) downward flux
	176	! (W.m-2) at the surface
	177	REAL,INTENT(OUT) :: fluxtop_lw(ngrid) ! outgoing total LW (thermal IR)
	178	! upward flux (W.m-2) at the top of the atm.
	179	REAL,INTENT(OUT) :: fluxsurf_dn_sw(ngrid,2) ! surface downward SW flux for
	180	! solar bands #1 and #2 (W.m-2)
	181	REAL,INTENT(OUT) :: fluxsurf_up_sw(ngrid,2) ! surface upward SW flux for
	182	! solar bands #1 and #2 (W.m-2)
	183	REAL,INTENT(OUT) :: fluxtop_dn_sw(ngrid,2) ! incoming downward SW flux for
	184	! solar bands #1 and #2 (W.m-2) at top of atm.
	185	REAL,INTENT(OUT) :: fluxtop_up_sw(ngrid,2) ! outgoing upward SW flux for
	186	! solar bands #1 and #2 (W.m-2) at top of atm.
[2415]	187	REAL,INTENT(OUT) :: tau_pref_scenario(ngrid) ! prescribed dust column
	188	! visible opacity at odpref from scenario
	189	REAL,INTENT(OUT) :: tau_pref_gcm(ngrid) ! computed dust column
	190	! visible opacity at odpref in the GCM
[2685]	191	REAL,INTENT(OUT) :: tau(ngrid,naerkind) ! Column visible optical depth
	192	! for each aerosol
[1047]	193	REAL,INTENT(OUT) :: taucloudtes(ngrid)! Cloud opacity at infrared
[520]	194	! reference wavelength using
	195	! Qabs instead of Qext
	196	! (direct comparison with TES)
[2685]	197	REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol extinction
	198	! optical depth at reference wavelength "longrefvis",
	199	! set in dimradmars_h, for each kind of aerosol
[2252]	200	REAL,INTENT(INOUT) :: dsodust(ngrid,nlayer)
[1047]	201	REAL,INTENT(OUT) :: rdust(ngrid,nlayer) ! Dust geometric mean radius (m)
	202	REAL,INTENT(OUT) :: rice(ngrid,nlayer) ! Ice geometric mean radius (m)
	203	REAL,INTENT(OUT) :: nuice(ngrid,nlayer) ! Estimated effective variance
[2459]	204	double precision,INTENT(OUT) :: riceco2(ngrid,nlayer) ! CO2 ice mean radius(m)
[2447]	205	REAL,INTENT(OUT) :: nuiceco2(ngrid,nlayer) ! Effective variance
[1047]	206	REAL,INTENT(IN) :: co2ice(ngrid) ! co2 ice surface layer (kg.m-2)
[1974]	207
	208	c rocket dust storm
	209	LOGICAL,INTENT(IN) :: clearatm ! true for background dust
	210	REAL,INTENT(IN) :: totstormfract(ngrid) ! dust storm mesh fraction
	211	REAL,INTENT(OUT) :: rstormdust(ngrid,nlayer) ! Storm dust geometric mean radius (m)
[2415]	212	REAL,INTENT(OUT) :: dsords(ngrid,nlayer) ! density scaled opacity for rocket dust storm dust
[2246]	213
[2628]	214	c entrainment by mountain top dust flows
	215	LOGICAL, INTENT(IN) :: nohmons ! true for background dust
[2199]	216	REAL,INTENT(OUT) :: rtopdust(ngrid,nlayer) ! Topdust geometric mean radius (m)
[2415]	217	REAL,INTENT(OUT) :: dsotop(ngrid,nlayer) ! density scaled opacity for topmons dust
[2246]	218
[1711]	219	c sub-grid scale water ice clouds
[1974]	220	LOGICAL,INTENT(IN) :: clearsky
	221	REAL,INTENT(IN) :: totcloudfrac(ngrid)
[353]	222
[38]	223	c
	224	c Local variables :
	225	c -----------------
	226
	227	INTEGER j,l,ig,n,ich
	228	INTEGER aer_count,iaer
	229	INTEGER jd,ig0,nd
	230
	231	real cste_mars ! solar constant on Mars (Wm-2)
[1047]	232	REAL ptlev(ngrid,nlayer+1)
[38]	233
[1774]	234	INTEGER :: ndomain
[38]	235
	236	c Thermal IR net radiative budget (W m-2)
	237
	238	real znetrad(ndomainsz,nflev)
	239
	240	real zfluxd_sw(ndomainsz,nflev+1,2)
	241	real zfluxu_sw(ndomainsz,nflev+1,2)
	242
	243	REAL zplev(ndomainsz,nflev+1)
	244	REAL zztlev(ndomainsz,nflev+1)
	245	REAL zplay(ndomainsz,nflev)
	246	REAL zt(ndomainsz,nflev)
	247	REAL zaerosol(ndomainsz,nflev,naerkind)
	248	REAL zalbedo(ndomainsz,2)
	249	REAL zdp(ndomainsz,nflev)
	250	REAL zdt0(ndomainsz)
	251
	252	REAL zzdtlw(ndomainsz,nflev)
	253	REAL zzdtsw(ndomainsz,nflev)
	254	REAL zzflux(ndomainsz,6)
	255	real zrmuz
	256
	257	REAL :: zQVISsQREF3d(ndomainsz,nflev,nsun,naerkind)
	258	REAL :: zomegaVIS3d(ndomainsz,nflev,nsun,naerkind)
	259	REAL :: zgVIS3d(ndomainsz,nflev,nsun,naerkind)
	260
	261	REAL :: zQIRsQREF3d(ndomainsz,nflev,nir,naerkind)
	262	REAL :: zomegaIR3d(ndomainsz,nflev,nir,naerkind)
	263	REAL :: zgIR3d(ndomainsz,nflev,nir,naerkind)
	264
	265	c Aerosol size distribution
	266	REAL :: reffrad(ngrid,nlayer,naerkind)
	267	REAL :: nueffrad(ngrid,nlayer,naerkind)
	268	c Aerosol optical properties
[1047]	269	REAL :: QVISsQREF3d(ngrid,nlayer,nsun,naerkind)
	270	REAL :: omegaVIS3d(ngrid,nlayer,nsun,naerkind)
	271	REAL :: gVIS3d(ngrid,nlayer,nsun,naerkind)
[38]	272
[1047]	273	REAL :: QIRsQREF3d(ngrid,nlayer,nir,naerkind)
	274	REAL :: omegaIR3d(ngrid,nlayer,nir,naerkind)
	275	REAL :: gIR3d(ngrid,nlayer,nir,naerkind)
[38]	276
[1047]	277	REAL :: QREFvis3d(ngrid,nlayer,naerkind)
[2246]	278	! QREFvis3d : Extinction efficiency at the VISible reference wavelength
[1047]	279	REAL :: QREFir3d(ngrid,nlayer,naerkind)
[2246]	280	! QREFir3d : Extinction efficiency at the InfraRed reference wavelength
[38]	281
[1047]	282	REAL :: omegaREFvis3d(ngrid,nlayer,naerkind)
	283	REAL :: omegaREFir3d(ngrid,nlayer,naerkind)
[38]	284
	285	c local saved variables
	286	c ---------------------
	287
	288	real zco2 ! volume fraction of CO2 in Mars atmosphere
[2584]	289	!$OMP THREADPRIVATE(zco2)
[38]	290	DATA zco2/0.95/
	291	SAVE zco2
	292
	293	LOGICAL firstcall
[2584]	294	!$OMP THREADPRIVATE(firstcall)
[38]	295	DATA firstcall/.true./
	296	SAVE firstcall
	297
[2584]	298
	299
[38]	300	c----------------------------------------------------------------------
	301
	302	c Initialisation
	303	c --------------
	304
[1776]	305	! compute ndomain
	306	! AS: moved out of firstcall to allow nesting+evoluting domain
	307	! ------------------------------------------------------------
	308	ndomain= (ngrid-1) / ndomainsz + 1
	309
[38]	310	IF (firstcall) THEN
	311
[1776]	312	write(,) 'Splitting radiative calculations: ',
	313	$ ' ngrid,ndomainsz,ndomain',
	314	$ ngrid,ndomainsz,ndomain
	315
[38]	316	c Assign a number to the different scatterers
	317	c -------------------------------------------
	318
	319	iaer_dust_conrath=0
	320	iaer_dust_doubleq=0
	321	iaer_dust_submicron=0
	322	iaer_h2o_ice=0
[2447]	323	iaer_co2_ice=0
[1974]	324	iaer_stormdust_doubleq=0
[2199]	325	iaer_topdust_doubleq=0
[38]	326
	327	aer_count=0
	328	if (.NOT.active) then
	329	do iaer=1,naerkind
	330	if (name_iaer(iaer).eq."dust_conrath") then
	331	iaer_dust_conrath = iaer
	332	aer_count = aer_count + 1
	333	endif
	334	enddo
	335	endif
	336	if (doubleq.AND.active) then
	337	do iaer=1,naerkind
	338	if (name_iaer(iaer).eq."dust_doubleq") then
	339	iaer_dust_doubleq = iaer
	340	aer_count = aer_count + 1
	341	endif
	342	enddo
	343	endif
	344	if (submicron.AND.active) then
	345	do iaer=1,naerkind
	346	if (name_iaer(iaer).eq."dust_submicron") then
	347	iaer_dust_submicron = iaer
	348	aer_count = aer_count + 1
	349	endif
	350	enddo
	351	endif
	352	if (water.AND.activice) then
	353	do iaer=1,naerkind
	354	if (name_iaer(iaer).eq."h2o_ice") then
	355	iaer_h2o_ice = iaer
	356	aer_count = aer_count + 1
	357	endif
	358	enddo
	359	endif
[2447]	360	if (co2clouds.AND.activeco2ice) then
	361	do iaer=1,naerkind
	362	if (name_iaer(iaer).eq."co2_ice") then
	363	iaer_co2_ice = iaer
	364	aer_count = aer_count + 1
	365	endif
	366	enddo
	367	endif
[1974]	368	if (rdstorm.AND.active) then
	369	do iaer=1,naerkind
	370	if (name_iaer(iaer).eq."stormdust_doubleq") then
	371	iaer_stormdust_doubleq = iaer
	372	aer_count = aer_count + 1
	373	endif
	374	enddo
	375	end if
[2628]	376	if (topflows.AND.active) then
[2199]	377	do iaer=1,naerkind
	378	if (name_iaer(iaer).eq."topdust_doubleq") then
	379	iaer_topdust_doubleq = iaer
	380	aer_count = aer_count + 1
	381	endif
	382	enddo
	383	end if
[38]	384
	385	c Check that we identified all tracers:
	386	if (aer_count.ne.naerkind) then
	387	write(,) "callradite: found only ",aer_count," scatterers"
	388	write(,) " expected ",naerkind
	389	write(,) "please make sure that the number of"
[1047]	390	write(,) "scatterers in scatterers.h, the names"
[38]	391	write(,) "in callradite.F, and the flags in"
	392	write(,) "callphys.def are all consistent!"
	393	do iaer=1,naerkind
	394	write(,)' ',iaer,' ',trim(name_iaer(iaer))
	395	enddo
[2398]	396	call abort_physic("callradite","incoherent scatterers",1)
[38]	397	else
	398	write(,) "callradite: found all scatterers, namely:"
	399	do iaer=1,naerkind
	400	write(,)' ',iaer,' ',trim(name_iaer(iaer))
	401	enddo
	402	endif
	403	c -------------------------------------------
	404
	405	gcp = g/cpp
	406
[2584]	407
[38]	408	c Loading the optical properties in external look-up tables:
[2584]	409
[38]	410	CALL SUAER
[2584]	411
[1047]	412	! CALL SULW ! this step is now done in ini_yomlw_h
[38]	413
[1047]	414	if (ngrid .EQ. 1) then
[38]	415	if (ndomainsz .NE. 1) then
	416	print*
	417	print*,'ATTENTION !!!'
	418	print*,'pour tourner en 1D, '
[1047]	419	print*,'fixer ndomainsz=1 dans phymars/dimradmars_h'
[38]	420	print*
	421	call exit(1)
	422	endif
	423	endif
[1774]	424
[38]	425	firstcall=.false.
	426	END IF
	427
	428	c Computing aerosol optical properties and opacity
	429	c ------------------------------------------------
	430	c Updating aerosol size distributions:
	431	CALL updatereffrad(ngrid,nlayer,
[2199]	432	& rdust,rstormdust,rtopdust,rice,nuice,
[2447]	433	& reffrad,nueffrad, riceco2, nuiceco2,
[2494]	434	& pq,tauscaling,tau,pplay, pt)
[38]	435	c Computing 3D scattering parameters:
[2494]	436	gVIS3d(:,:,:,:) = 0.
[38]	437	CALL aeroptproperties(ngrid,nlayer,reffrad,nueffrad,
	438	& QVISsQREF3d,omegaVIS3d,gVIS3d,
	439	& QIRsQREF3d,omegaIR3d,gIR3d,
	440	& QREFvis3d,QREFir3d,
	441	& omegaREFvis3d,omegaREFir3d)
	442	c Computing aerosol optical depth in each layer:
	443	CALL aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
[2643]	444	& pq,pt,tauscaling,dust_rad_adjust,IRtoVIScoef,
	445	& tau_pref_scenario,tau_pref_gcm,tau,taucloudtes,
	446	& aerosol,dsodust,reffrad,
[1974]	447	& QREFvis3d,QREFir3d,omegaREFir3d,
[2246]	448	& totstormfract,clearatm,dsords,dsotop,
[2628]	449	& nohmons,clearsky,totcloudfrac)
	450
[38]	451	c Starting loop on sub-domain
	452	c ----------------------------
[2494]	453	zgVIS3d(:,:,:,:) = 0.
	454	zfluxd_sw(:,:,:) = 0.
	455	zfluxu_sw(:,:,:) = 0.
	456	zQVISsQREF3d(:,:,:,:) = 0.
	457	zomegaVIS3d(:,:,:,:) = 0.
[38]	458	DO jd=1,ndomain
	459	ig0=(jd-1)*ndomainsz
	460	if (jd.eq.ndomain) then
[1047]	461	nd=ngrid-ig0
[38]	462	else
	463	nd=ndomainsz
	464	endif
	465
	466	c Spliting input variable in sub-domain input variables
	467	c ---------------------------------------------------
	468
	469	do l=1,nlaylte
	470	do ig = 1,nd
	471	do iaer = 1, naerkind
	472	do ich = 1, nsun
	473	zQVISsQREF3d(ig,l,ich,iaer) =
	474	& QVISsQREF3d(ig0+ig,l,ich,iaer)
	475	zomegaVIS3d(ig,l,ich,iaer) =
	476	& omegaVIS3d(ig0+ig,l,ich,iaer)
	477	zgVIS3d(ig,l,ich,iaer) =
	478	& gVIS3d(ig0+ig,l,ich,iaer)
	479	enddo
	480	do ich = 1, nir
	481	zQIRsQREF3d(ig,l,ich,iaer) =
	482	& QIRsQREF3d(ig0+ig,l,ich,iaer)
	483	zomegaIR3d(ig,l,ich,iaer) =
	484	& omegaIR3d(ig0+ig,l,ich,iaer)
	485	zgIR3d(ig,l,ich,iaer) =
	486	& gIR3d(ig0+ig,l,ich,iaer)
	487	enddo
	488	enddo
	489	enddo
	490	enddo
[2494]	491	zplev(:,:) = 0.
[38]	492	do l=1,nlaylte+1
	493	do ig = 1,nd
	494	zplev(ig,l) = pplev(ig0+ig,l)
	495	enddo
	496	enddo
[2494]	497	zdp(:,:) = 0.
	498
[38]	499	do l=1,nlaylte
	500	do ig = 1,nd
	501	zplay(ig,l) = pplay(ig0+ig,l)
	502	zt(ig,l) = pt(ig0+ig,l)
	503	c Thickness of each layer (Pa) :
	504	zdp(ig,l)= pplev(ig0+ig,l) - pplev(ig0+ig,l+1)
	505	enddo
	506	enddo
[2494]	507	zaerosol(:,:,:) = 0.
[38]	508	do n=1,naerkind
	509	do l=1,nlaylte
	510	do ig=1,nd
	511	zaerosol(ig,l,n) = aerosol(ig0+ig,l,n)
	512	enddo
	513	enddo
	514	enddo
[2494]	515	zalbedo(:,:) = 0.
[38]	516	do j=1,2
	517	do ig = 1,nd
	518	zalbedo(ig,j) = albedo(ig0+ig,j)
	519	enddo
	520	enddo
	521
	522	c Intermediate levels: (computing tlev)
	523	c ---------------------------------------
	524	c Extrapolation for the air temperature above the surface
	525	DO ig=1,nd
	526	zztlev(ig,1)=zt(ig,1)+
	527	s (zplev(ig,1)-zplay(ig,1))*
	528	s (zt(ig,1)-zt(ig,2))/(zplay(ig,1)-zplay(ig,2))
	529
	530	zdt0(ig) = tsurf(ig0+ig) - zztlev(ig,1)
	531	ENDDO
	532
	533	DO l=2,nlaylte
	534	DO ig=1, nd
	535	zztlev(ig,l)=0.5*(zt(ig,l-1)+zt(ig,l))
	536	ENDDO
	537	ENDDO
	538
	539	DO ig=1, nd
	540	zztlev(ig,nlaylte+1)=zt(ig,nlaylte)
	541	ENDDO
	542
	543
	544	c Longwave ("lw") radiative transfer (= thermal infrared)
	545	c -------------------------------------------------------
	546	call lwmain (ig0,icount,nd,nflev
	547	. ,zdp,zdt0,emis(ig0+1),zplev,zztlev,zt
	548	. ,zaerosol,zzdtlw
	549	. ,fluxsurf_lw(ig0+1),fluxtop_lw(ig0+1)
	550	. ,znetrad
[353]	551	& ,zQIRsQREF3d,zomegaIR3d,zgIR3d
	552	& ,co2ice(ig0+1))
[38]	553
	554	c Shortwave ("sw") radiative transfer (= solar radiation)
	555	c -------------------------------------------------------
	556	c Mars solar constant (W m-2)
[2678]	557	c flux_1AU = 1370 W.m-2 is the solar constant at 1 AU.
	558	cste_mars=flux_1AU/(dist_sol*dist_sol)
[2494]	559	zzdtsw(:,:) = 0.
[38]	560	call swmain ( nd, nflev,
	561	S cste_mars, zalbedo,
	562	S mu0(ig0+1), zdp, zplev, zaerosol, fract(ig0+1),
	563	S zzdtsw, zfluxd_sw, zfluxu_sw,
	564	& zQVISsQREF3d,zomegaVIS3d,zgVIS3d)
	565	c ------------------------------------------------------------
	566	c Un-spliting output variable from sub-domain input variables
	567	c ------------------------------------------------------------
	568
	569	do l=1,nlaylte
	570	do ig = 1,nd
	571	dtlw(ig0+ig,l) = zzdtlw(ig,l)
	572	dtsw(ig0+ig,l) = zzdtsw(ig,l)
	573	enddo
	574	enddo
	575
[2494]	576	ptlev(:, :) = 0.
[38]	577	do l=1,nlaylte+1
	578	do ig = 1,nd
	579	ptlev(ig0+ig,l) = zztlev(ig,l)
	580	enddo
	581	enddo
	582
[2685]	583	! copy SW fluxes at surface and TOA
[38]	584	do ig = 1,nd
[2685]	585	! surface downward SW flux
	586	fluxsurf_dn_sw(ig0+ig,1) = zfluxd_sw(ig,1,1)
	587	fluxsurf_dn_sw(ig0+ig,2) = zfluxd_sw(ig,1,2)
	588	! surface upward SW flux
	589	fluxsurf_up_sw(ig0+ig,1) = zfluxu_sw(ig,1,1)
	590	fluxsurf_up_sw(ig0+ig,2) = zfluxu_sw(ig,1,2)
	591	! downward SW flux at top of atmosphere
	592	fluxtop_dn_sw(ig0+ig,1) = zfluxd_sw(ig,nlaylte+1,1)
	593	fluxtop_dn_sw(ig0+ig,2) = zfluxd_sw(ig,nlaylte+1,2)
	594	! upward SW flux at top of atmosphere
	595	fluxtop_up_sw(ig0+ig,1) = zfluxu_sw(ig,nlaylte+1,1)
	596	fluxtop_up_sw(ig0+ig,2) = zfluxu_sw(ig,nlaylte+1,2)
[38]	597	enddo
	598
	599	ENDDO ! (boucle jd=1, ndomain)
	600
	601	c Zero tendencies for any remaining layers between nlaylte and nlayer
	602	if (nlayer.gt.nlaylte) then
	603	do l = nlaylte+1, nlayer
	604	do ig = 1, ngrid
	605	dtlw(ig, l) = 0.
	606	dtsw(ig, l) = 0.
	607	enddo
	608	enddo
	609	endif
[3468]	610	c Output for debugging
[38]	611	c --------------------------------
	612	c Write all nlayer layers, even though only nlaylte layers may have
	613	c non-zero tendencies.
	614
[3468]	615	c PRINT*,'Diagnotique for the radiation'
	616	c PRINT*,'albedo, emissiv, mu0,fract,fluxsurf_lw,fluxsurf_sw'
	617	c PRINT*,albedo(igout,1),emis(igout),mu0(igout),
	618	c s fract(igout), fluxsurf_lw(igout),
	619	c $ fluxsurf_dn_sw(igout,1)+fluxsurf_dn_sw(igout,2)
	620	c PRINT*,'Tlay Tlev Play Plev dT/dt SW dT/dt LW (K/s)'
	621	c PRINT*,'daysec',daysec
	622	c DO l=1,nlayer
	623	c PRINT*,pt(igout,l),ptlev(igout,l),
	624	c s pplay(igout,l),pplev(igout,l),
	625	c s dtsw(igout,l),dtlw(igout,l)
	626	c ENDDO
[38]	627
	628
[1711]	629	END SUBROUTINE callradite
	630
	631	END MODULE callradite_mod

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

Download in other formats:

Original Format