Context Navigation

radlwsw_m.F90 @ 2095

Last change on this file since 2095 was 2043, checked in by Laurent Fairhead, 11 years ago

Pour compilation avec ifort

For ifort compilation

Baek

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:keywords set to Author Date Id Revi

File size: 41.2 KB

Rev	Line
[2003]	1	!
	2	! $Id: radlwsw_m.F90 2043 2014-05-14 11:09:20Z fairhead $
	3	!
[1687]	4	module radlwsw_m
	5
	6	IMPLICIT NONE
	7
	8	contains
	9
	10	SUBROUTINE radlwsw( &
	11	dist, rmu0, fract, &
	12	paprs, pplay,tsol,alb1, alb2, &
	13	t,q,wo,&
	14	cldfra, cldemi, cldtaupd,&
	15	ok_ade, ok_aie, flag_aerosol,&
[1764]	16	flag_aerosol_strat,&
[1687]	17	tau_aero, piz_aero, cg_aero,&
[2003]	18	tau_aero_rrtm, piz_aero_rrtm, cg_aero_rrtm,& ! rajoute par OB pour RRTM
[1687]	19	cldtaupi, new_aod, &
	20	qsat, flwc, fiwc, &
[1989]	21	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
[1687]	22	heat,heat0,cool,cool0,radsol,albpla,&
	23	topsw,toplw,solsw,sollw,&
	24	sollwdown,&
	25	topsw0,toplw0,solsw0,sollw0,&
	26	lwdn0, lwdn, lwup0, lwup,&
	27	swdn0, swdn, swup0, swup,&
	28	topswad_aero, solswad_aero,&
	29	topswai_aero, solswai_aero, &
	30	topswad0_aero, solswad0_aero,&
	31	topsw_aero, topsw0_aero,&
	32	solsw_aero, solsw0_aero, &
[1989]	33	topswcf_aero, solswcf_aero,&
	34	ZLWFT0_i, ZFLDN0, ZFLUP0,&
	35	ZSWFT0_i, ZFSDN0, ZFSUP0)
[1687]	36
	37
	38
	39	USE DIMPHY
	40	USE assert_m, ONLY : assert
	41	USE infotrac, ONLY : type_trac
[1989]	42	USE write_field_phy
[1687]	43	#ifdef REPROBUS
	44	USE CHEM_REP, ONLY : solaireTIME, ok_SUNTIME, ndimozon
	45	#endif
[1989]	46	#ifdef CPP_RRTM
	47	! modules necessaires au rayonnement
	48	! -----------------------------------------
	49	! USE YOMCST , ONLY : RG ,RD ,RTT ,RPI
	50	! USE YOERAD , ONLY : NSW ,LRRTM ,LINHOM , LCCNL,LCCNO,
	51	! USE YOERAD , ONLY : NSW ,LRRTM ,LCCNL ,LCCNO ,&
	52	! NSW mis dans .def MPL 20140211
	53	USE YOERAD , ONLY : LRRTM ,LCCNL ,LCCNO ,&
	54	NRADIP , NRADLP , NICEOPT, NLIQOPT ,RCCNLND , RCCNSEA
	55	USE YOELW , ONLY : NSIL ,NTRA ,NUA ,TSTAND ,XP
	56	USE YOESW , ONLY : RYFWCA ,RYFWCB ,RYFWCC ,RYFWCD,&
	57	RYFWCE ,RYFWCF ,REBCUA ,REBCUB ,REBCUC,&
	58	REBCUD ,REBCUE ,REBCUF ,REBCUI ,REBCUJ,&
	59	REBCUG ,REBCUH ,RHSAVI ,RFULIO ,RFLAA0,&
	60	RFLAA1 ,RFLBB0 ,RFLBB1 ,RFLBB2 ,RFLBB3,&
	61	RFLCC0 ,RFLCC1 ,RFLCC2 ,RFLCC3 ,RFLDD0,&
	62	RFLDD1 ,RFLDD2 ,RFLDD3 ,RFUETA ,RASWCA,&
	63	RASWCB ,RASWCC ,RASWCD ,RASWCE ,RASWCF
	64	! & RASWCB ,RASWCC ,RASWCD ,RASWCE ,RASWCF, RLINLI
	65	USE YOERDU , ONLY : NUAER ,NTRAER ,REPLOG ,REPSC ,REPSCW ,DIFF
[2043]	66	! USE YOETHF , ONLY : RTICE
[1989]	67	USE YOERRTWN , ONLY : DELWAVE ,TOTPLNK
	68	USE YOMPHY3 , ONLY : RII0
	69	#endif
[1687]	70
	71	!======================================================================
	72	! Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
	73	! Objet: interface entre le modele et les rayonnements
	74	! Arguments:
	75	! dist-----input-R- distance astronomique terre-soleil
	76	! rmu0-----input-R- cosinus de l'angle zenithal
	77	! fract----input-R- duree d'ensoleillement normalisee
	78	! co2_ppm--input-R- concentration du gaz carbonique (en ppm)
	79	! paprs----input-R- pression a inter-couche (Pa)
	80	! pplay----input-R- pression au milieu de couche (Pa)
	81	! tsol-----input-R- temperature du sol (en K)
	82	! alb1-----input-R- albedo du sol(entre 0 et 1) dans l'interval visible
	83	! alb2-----input-R- albedo du sol(entre 0 et 1) dans l'interval proche infra-rouge
	84	! t--------input-R- temperature (K)
	85	! q--------input-R- vapeur d'eau (en kg/kg)
	86	! cldfra---input-R- fraction nuageuse (entre 0 et 1)
	87	! cldtaupd---input-R- epaisseur optique des nuages dans le visible (present-day value)
	88	! cldemi---input-R- emissivite des nuages dans l'IR (entre 0 et 1)
	89	! ok_ade---input-L- apply the Aerosol Direct Effect or not?
	90	! ok_aie---input-L- apply the Aerosol Indirect Effect or not?
	91	! flag_aerosol-input-I- aerosol flag from 0 to 6
[1764]	92	! flag_aerosol_strat-input-I- use stratospheric aerosols flag (T/F)
[1687]	93	! tau_ae, piz_ae, cg_ae-input-R- aerosol optical properties (calculated in aeropt.F)
	94	! cldtaupi-input-R- epaisseur optique des nuages dans le visible
	95	! calculated for pre-industrial (pi) aerosol concentrations, i.e. with smaller
	96	! droplet concentration, thus larger droplets, thus generally cdltaupi cldtaupd
	97	! it is needed for the diagnostics of the aerosol indirect radiative forcing
	98	!
	99	! heat-----output-R- echauffement atmospherique (visible) (K/jour)
	100	! cool-----output-R- refroidissement dans l'IR (K/jour)
	101	! radsol---output-R- bilan radiatif net au sol (W/m**2) (+ vers le bas)
	102	! albpla---output-R- albedo planetaire (entre 0 et 1)
	103	! topsw----output-R- flux solaire net au sommet de l'atm.
	104	! toplw----output-R- ray. IR montant au sommet de l'atmosphere
	105	! solsw----output-R- flux solaire net a la surface
	106	! sollw----output-R- ray. IR montant a la surface
	107	! solswad---output-R- ray. solaire net absorbe a la surface (aerosol dir)
	108	! topswad---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol dir)
	109	! solswai---output-R- ray. solaire net absorbe a la surface (aerosol ind)
	110	! topswai---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol ind)
	111	!
	112	! ATTENTION: swai and swad have to be interpreted in the following manner:
	113	! ---------
	114	! ok_ade=F & ok_aie=F -both are zero
	115	! ok_ade=T & ok_aie=F -aerosol direct forcing is F_{AD} = topsw-topswad
	116	! indirect is zero
	117	! ok_ade=F & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
	118	! direct is zero
	119	! ok_ade=T & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
	120	! aerosol direct forcing is F_{AD} = topswai-topswad
	121	!
[1989]	122	! --------- RRTM: output RECMWFL
	123	! ZEMTD (KPROMA,KLEV+1) ; TOTAL DOWNWARD LONGWAVE EMISSIVITY
	124	! ZEMTU (KPROMA,KLEV+1) ; TOTAL UPWARD LONGWAVE EMISSIVITY
	125	! ZTRSO (KPROMA,KLEV+1) ; TOTAL SHORTWAVE TRANSMISSIVITY
	126	! ZTH (KPROMA,KLEV+1) ; HALF LEVEL TEMPERATURE
	127	! ZCTRSO(KPROMA,2) ; CLEAR-SKY SHORTWAVE TRANSMISSIVITY
	128	! ZCEMTR(KPROMA,2) ; CLEAR-SKY NET LONGWAVE EMISSIVITY
	129	! ZTRSOD(KPROMA) ; TOTAL-SKY SURFACE SW TRANSMISSITY
	130	! ZLWFC (KPROMA,2) ; CLEAR-SKY LONGWAVE FLUXES
	131	! ZLWFT (KPROMA,KLEV+1) ; TOTAL-SKY LONGWAVE FLUXES
	132	! ZLWFT0(KPROMA,KLEV+1) ; CLEAR-SKY LONGWAVE FLUXES ! added by MPL 090109
	133	! ZSWFC (KPROMA,2) ; CLEAR-SKY SHORTWAVE FLUXES
	134	! ZSWFT (KPROMA,KLEV+1) ; TOTAL-SKY SHORTWAVE FLUXES
	135	! ZSWFT0(KPROMA,KLEV+1) ; CLEAR-SKY SHORTWAVE FLUXES ! added by MPL 090109
	136	! ZFLUX (KLON,2,KLEV+1) ; TOTAL LW FLUXES 1=up, 2=DWN ! added by MPL 080411
	137	! ZFLUC (KLON,2,KLEV+1) ; CLEAR SKY LW FLUXES ! added by MPL 080411
	138	! ZFSDWN(klon,KLEV+1) ; TOTAL SW DWN FLUXES ! added by MPL 080411
	139	! ZFCDWN(klon,KLEV+1) ; CLEAR SKY SW DWN FLUXES ! added by MPL 080411
	140	! ZFSUP (klon,KLEV+1) ; TOTAL SW UP FLUXES ! added by MPL 080411
	141	! ZFCUP (klon,KLEV+1) ; CLEAR SKY SW UP FLUXES ! added by MPL 080411
[1687]	142
	143	!======================================================================
	144
	145	! ====================================================================
	146	! Adapte au modele de chimie INCA par Celine Deandreis & Anne Cozic -- 2009
	147	! 1 = ZERO
	148	! 2 = AER total
	149	! 3 = NAT
	150	! 4 = BC
	151	! 5 = SO4
	152	! 6 = POM
	153	! 7 = DUST
	154	! 8 = SS
	155	! 9 = NO3
	156	!
	157	! ====================================================================
	158	include "YOETHF.h"
	159	include "YOMCST.h"
	160	include "clesphys.h"
	161	include "iniprint.h"
	162
	163	! Input arguments
	164	REAL, INTENT(in) :: dist
	165	REAL, INTENT(in) :: rmu0(KLON), fract(KLON)
	166	REAL, INTENT(in) :: paprs(KLON,KLEV+1), pplay(KLON,KLEV)
	167	REAL, INTENT(in) :: alb1(KLON), alb2(KLON), tsol(KLON)
	168	REAL, INTENT(in) :: t(KLON,KLEV), q(KLON,KLEV)
	169
	170	REAL, INTENT(in):: wo(:, :, :) ! dimension(KLON,KLEV, 1 or 2)
	171	! column-density of ozone in a layer, in kilo-Dobsons
	172	! "wo(:, :, 1)" is for the average day-night field,
	173	! "wo(:, :, 2)" is for daylight time.
	174
	175	LOGICAL, INTENT(in) :: ok_ade, ok_aie ! switches whether to use aerosol direct (indirect) effects or not
[1989]	176	LOGICAL :: lldebug
[1687]	177	INTEGER, INTENT(in) :: flag_aerosol ! takes value 0 (no aerosol) or 1 to 6 (aerosols)
[1764]	178	LOGICAL, INTENT(in) :: flag_aerosol_strat ! use stratospheric aerosols
[1687]	179	REAL, INTENT(in) :: cldfra(KLON,KLEV), cldemi(KLON,KLEV), cldtaupd(KLON,KLEV)
	180	REAL, INTENT(in) :: tau_aero(KLON,KLEV,9,2) ! aerosol optical properties (see aeropt.F)
	181	REAL, INTENT(in) :: piz_aero(KLON,KLEV,9,2) ! aerosol optical properties (see aeropt.F)
	182	REAL, INTENT(in) :: cg_aero(KLON,KLEV,9,2) ! aerosol optical properties (see aeropt.F)
[2003]	183	!--OB
	184	REAL, INTENT(in) :: tau_aero_rrtm(KLON,KLEV,2,NSW) ! aerosol optical properties RRTM
	185	REAL, INTENT(in) :: piz_aero_rrtm(KLON,KLEV,2,NSW) ! aerosol optical properties RRTM
	186	REAL, INTENT(in) :: cg_aero_rrtm(KLON,KLEV,2,NSW) ! aerosol optical properties RRTM
	187	!--OB fin
[1687]	188	REAL, INTENT(in) :: cldtaupi(KLON,KLEV) ! cloud optical thickness for pre-industrial aerosol concentrations
	189	LOGICAL, INTENT(in) :: new_aod ! flag pour retrouver les resultats exacts de l'AR4 dans le cas ou l'on ne travaille qu'avec les sulfates
	190	REAL, INTENT(in) :: qsat(klon,klev) ! Variable pour iflag_rrtm=1
	191	REAL, INTENT(in) :: flwc(klon,klev) ! Variable pour iflag_rrtm=1
	192	REAL, INTENT(in) :: fiwc(klon,klev) ! Variable pour iflag_rrtm=1
[1989]	193	REAL, INTENT(in) :: ref_liq(klon,klev) ! cloud droplet radius present-day from newmicro
	194	REAL, INTENT(in) :: ref_ice(klon,klev) ! ice crystal radius present-day from newmicro
	195	REAL, INTENT(in) :: ref_liq_pi(klon,klev) ! cloud droplet radius pre-industrial from newmicro
	196	REAL, INTENT(in) :: ref_ice_pi(klon,klev) ! ice crystal radius pre-industrial from newmicro
[1687]	197
	198	! Output arguments
	199	REAL, INTENT(out) :: heat(KLON,KLEV), cool(KLON,KLEV)
	200	REAL, INTENT(out) :: heat0(KLON,KLEV), cool0(KLON,KLEV)
	201	REAL, INTENT(out) :: radsol(KLON), topsw(KLON), toplw(KLON)
	202	REAL, INTENT(out) :: solsw(KLON), sollw(KLON), albpla(KLON)
	203	REAL, INTENT(out) :: topsw0(KLON), toplw0(KLON), solsw0(KLON), sollw0(KLON)
	204	REAL, INTENT(out) :: sollwdown(KLON)
	205	REAL, INTENT(out) :: swdn(KLON,kflev+1),swdn0(KLON,kflev+1)
	206	REAL, INTENT(out) :: swup(KLON,kflev+1),swup0(KLON,kflev+1)
	207	REAL, INTENT(out) :: lwdn(KLON,kflev+1),lwdn0(KLON,kflev+1)
	208	REAL, INTENT(out) :: lwup(KLON,kflev+1),lwup0(KLON,kflev+1)
	209	REAL, INTENT(out) :: topswad_aero(KLON), solswad_aero(KLON) ! output: aerosol direct forcing at TOA and surface
	210	REAL, INTENT(out) :: topswai_aero(KLON), solswai_aero(KLON) ! output: aerosol indirect forcing atTOA and surface
	211	REAL, DIMENSION(klon), INTENT(out) :: topswad0_aero
	212	REAL, DIMENSION(klon), INTENT(out) :: solswad0_aero
	213	REAL, DIMENSION(kdlon,9), INTENT(out) :: topsw_aero
	214	REAL, DIMENSION(kdlon,9), INTENT(out) :: topsw0_aero
	215	REAL, DIMENSION(kdlon,9), INTENT(out) :: solsw_aero
	216	REAL, DIMENSION(kdlon,9), INTENT(out) :: solsw0_aero
	217	REAL, DIMENSION(kdlon,3), INTENT(out) :: topswcf_aero
	218	REAL, DIMENSION(kdlon,3), INTENT(out) :: solswcf_aero
[1989]	219	REAL, DIMENSION(kdlon,kflev+1), INTENT(out) :: ZSWFT0_i
	220	REAL, DIMENSION(kdlon,kflev+1), INTENT(out) :: ZLWFT0_i
[1687]	221
	222	! Local variables
	223	REAL(KIND=8) ZFSUP(KDLON,KFLEV+1)
	224	REAL(KIND=8) ZFSDN(KDLON,KFLEV+1)
	225	REAL(KIND=8) ZFSUP0(KDLON,KFLEV+1)
	226	REAL(KIND=8) ZFSDN0(KDLON,KFLEV+1)
	227	REAL(KIND=8) ZFLUP(KDLON,KFLEV+1)
	228	REAL(KIND=8) ZFLDN(KDLON,KFLEV+1)
	229	REAL(KIND=8) ZFLUP0(KDLON,KFLEV+1)
	230	REAL(KIND=8) ZFLDN0(KDLON,KFLEV+1)
	231	REAL(KIND=8) zx_alpha1, zx_alpha2
	232	INTEGER k, kk, i, j, iof, nb_gr
[1989]	233	INTEGER ist,iend,ktdia,kmode
[1687]	234	REAL(KIND=8) PSCT
	235	REAL(KIND=8) PALBD(kdlon,2), PALBP(kdlon,2)
[1989]	236	! MPL 06.01.09: pour RRTM, creation de PALBD_NEW et PALBP_NEW
	237	! avec NSW en deuxieme dimension
	238	REAL(KIND=8) PALBD_NEW(kdlon,NSW), PALBP_NEW(kdlon,NSW)
[1687]	239	REAL(KIND=8) PEMIS(kdlon), PDT0(kdlon), PVIEW(kdlon)
	240	REAL(KIND=8) PPSOL(kdlon), PDP(kdlon,KLEV)
	241	REAL(KIND=8) PTL(kdlon,kflev+1), PPMB(kdlon,kflev+1)
	242	REAL(KIND=8) PTAVE(kdlon,kflev)
	243	REAL(KIND=8) PWV(kdlon,kflev), PQS(kdlon,kflev)
	244
	245	real(kind=8) POZON(kdlon, kflev, size(wo, 3)) ! mass fraction of ozone
	246	! "POZON(:, :, 1)" is for the average day-night field,
	247	! "POZON(:, :, 2)" is for daylight time.
[1989]	248	!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6
	249	REAL(KIND=8) PAER(kdlon,kflev,6)
[1687]	250	REAL(KIND=8) PCLDLD(kdlon,kflev)
	251	REAL(KIND=8) PCLDLU(kdlon,kflev)
	252	REAL(KIND=8) PCLDSW(kdlon,kflev)
	253	REAL(KIND=8) PTAU(kdlon,2,kflev)
	254	REAL(KIND=8) POMEGA(kdlon,2,kflev)
	255	REAL(KIND=8) PCG(kdlon,2,kflev)
	256	REAL(KIND=8) zfract(kdlon), zrmu0(kdlon), zdist
	257	REAL(KIND=8) zheat(kdlon,kflev), zcool(kdlon,kflev)
	258	REAL(KIND=8) zheat0(kdlon,kflev), zcool0(kdlon,kflev)
	259	REAL(KIND=8) ztopsw(kdlon), ztoplw(kdlon)
	260	REAL(KIND=8) zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon)
	261	REAL(KIND=8) zsollwdown(kdlon)
	262	REAL(KIND=8) ztopsw0(kdlon), ztoplw0(kdlon)
	263	REAL(KIND=8) zsolsw0(kdlon), zsollw0(kdlon)
	264	REAL(KIND=8) zznormcp
	265	REAL(KIND=8) tauaero(kdlon,kflev,9,2) ! aer opt properties
	266	REAL(KIND=8) pizaero(kdlon,kflev,9,2)
	267	REAL(KIND=8) cgaero(kdlon,kflev,9,2)
	268	REAL(KIND=8) PTAUA(kdlon,2,kflev) ! present-day value of cloud opt thickness (PTAU is pre-industrial value), local use
	269	REAL(KIND=8) POMEGAA(kdlon,2,kflev) ! dito for single scatt albedo
	270	REAL(KIND=8) ztopswadaero(kdlon), zsolswadaero(kdlon) ! Aerosol direct forcing at TOAand surface
	271	REAL(KIND=8) ztopswad0aero(kdlon), zsolswad0aero(kdlon) ! Aerosol direct forcing at TOAand surface
	272	REAL(KIND=8) ztopswaiaero(kdlon), zsolswaiaero(kdlon) ! dito, indirect
	273	REAL(KIND=8) ztopsw_aero(kdlon,9), ztopsw0_aero(kdlon,9)
	274	REAL(KIND=8) zsolsw_aero(kdlon,9), zsolsw0_aero(kdlon,9)
	275	REAL(KIND=8) ztopswcf_aero(kdlon,3), zsolswcf_aero(kdlon,3)
[1989]	276	! real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2 deje declare dans physiq.F MPL 20130618
	277	!MPL input supplementaires pour RECMWFL
	278	! flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
	279	REAL(KIND=8) GEMU(klon)
	280	!MPL input RECMWFL:
	281	! Tableaux aux niveaux inverses pour respecter convention Arpege
	282	REAL(KIND=8) ref_liq_i(klon,klev) ! cloud droplet radius present-day from newmicro (inverted)
	283	REAL(KIND=8) ref_ice_i(klon,klev) ! ice crystal radius present-day from newmicro (inverted)
[2003]	284	!--OB
	285	REAL(KIND=8) ref_liq_pi_i(klon,klev) ! cloud droplet radius pre-industrial from newmicro (inverted)
	286	REAL(KIND=8) ref_ice_pi_i(klon,klev) ! ice crystal radius pre-industrial from newmicro (inverted)
	287	!--end OB
[1989]	288	REAL(KIND=8) paprs_i(klon,klev+1)
	289	REAL(KIND=8) pplay_i(klon,klev)
	290	REAL(KIND=8) cldfra_i(klon,klev)
	291	REAL(KIND=8) POZON_i(kdlon,kflev, size(wo, 3)) ! mass fraction of ozone
	292	! "POZON(:, :, 1)" is for the average day-night field,
	293	! "POZON(:, :, 2)" is for daylight time.
	294	!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6
	295	REAL(KIND=8) PAER_i(kdlon,kflev,6)
	296	REAL(KIND=8) PDP_i(klon,klev)
	297	REAL(KIND=8) t_i(klon,klev),q_i(klon,klev),qsat_i(klon,klev)
	298	REAL(KIND=8) flwc_i(klon,klev),fiwc_i(klon,klev)
	299	!MPL output RECMWFL:
	300	REAL(KIND=8) ZEMTD (klon,klev+1),ZEMTD_i (klon,klev+1)
	301	REAL(KIND=8) ZEMTU (klon,klev+1),ZEMTU_i (klon,klev+1)
	302	REAL(KIND=8) ZTRSO (klon,klev+1),ZTRSO_i (klon,klev+1)
	303	REAL(KIND=8) ZTH (klon,klev+1),ZTH_i (klon,klev+1)
	304	REAL(KIND=8) ZCTRSO(klon,2)
	305	REAL(KIND=8) ZCEMTR(klon,2)
	306	REAL(KIND=8) ZTRSOD(klon)
	307	REAL(KIND=8) ZLWFC (klon,2)
	308	REAL(KIND=8) ZLWFT (klon,klev+1),ZLWFT_i (klon,klev+1)
	309	REAL(KIND=8) ZSWFC (klon,2)
	310	REAL(KIND=8) ZSWFT (klon,klev+1),ZSWFT_i (klon,klev+1)
	311	REAL(KIND=8) ZFLUCDWN_i(klon,klev+1),ZFLUCUP_i(klon,klev+1)
[2003]	312	REAL(KIND=8) PPIZA_TOT(klon,klev,NSW)
	313	REAL(KIND=8) PCGA_TOT(klon,klev,NSW)
	314	REAL(KIND=8) PTAU_TOT(klon,klev,NSW)
	315	REAL(KIND=8) PPIZA_NAT(klon,klev,NSW)
	316	REAL(KIND=8) PCGA_NAT(klon,klev,NSW)
	317	REAL(KIND=8) PTAU_NAT(klon,klev,NSW)
[1989]	318	REAL(KIND=8) PSFSWDIR(klon,NSW)
	319	REAL(KIND=8) PSFSWDIF(klon,NSW)
	320	REAL(KIND=8) PFSDNN(klon)
	321	REAL(KIND=8) PFSDNV(klon)
	322	!MPL On ne redefinit pas les tableaux ZFLUX,ZFLUC,
	323	!MPL ZFSDWN,ZFCDWN,ZFSUP,ZFCUP car ils existent deja
	324	!MPL sous les noms de ZFLDN,ZFLDN0,ZFLUP,ZFLUP0,
	325	!MPL ZFSDN,ZFSDN0,ZFSUP,ZFSUP0
	326	REAL(KIND=8) ZFLUX_i (klon,2,klev+1)
	327	REAL(KIND=8) ZFLUC_i (klon,2,klev+1)
	328	REAL(KIND=8) ZFSDWN_i (klon,klev+1)
	329	REAL(KIND=8) ZFCDWN_i (klon,klev+1)
	330	REAL(KIND=8) ZFSUP_i (klon,klev+1)
	331	REAL(KIND=8) ZFCUP_i (klon,klev+1)
	332	! 3 lignes suivantes a activer pour CCMVAL (MPL 20100412)
	333	! REAL(KIND=8) RSUN(3,2)
	334	! REAL(KIND=8) SUN(3)
	335	! REAL(KIND=8) SUN_FRACT(2)
[1687]	336	real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
[2003]	337	CHARACTER (LEN=80) :: abort_message
	338	CHARACTER (LEN=80) :: modname='radlwsw_m'
[1687]	339
	340	call assert(size(wo, 1) == klon, size(wo, 2) == klev, "radlwsw wo")
	341	! initialisation
[1989]	342	ist=1
	343	iend=klon
	344	ktdia=1
	345	kmode=ist
[1687]	346	tauaero(:,:,:,:)=0.
	347	pizaero(:,:,:,:)=0.
	348	cgaero(:,:,:,:)=0.
[1989]	349	lldebug=.FALSE.
[1687]	350
	351	!
	352	!-------------------------------------------
	353	nb_gr = KLON / kdlon
	354	IF (nb_gr*kdlon .NE. KLON) THEN
	355	PRINT*, "kdlon mauvais:", KLON, kdlon, nb_gr
[1931]	356	call abort_gcm("radlwsw", "", 1)
[1687]	357	ENDIF
	358	IF (kflev .NE. KLEV) THEN
	359	PRINT*, "kflev differe de KLEV, kflev, KLEV"
[1931]	360	call abort_gcm("radlwsw", "", 1)
[1687]	361	ENDIF
	362	!-------------------------------------------
	363	DO k = 1, KLEV
	364	DO i = 1, KLON
	365	heat(i,k)=0.
	366	cool(i,k)=0.
	367	heat0(i,k)=0.
	368	cool0(i,k)=0.
	369	ENDDO
	370	ENDDO
	371	!
	372	zdist = dist
	373	!
	374	PSCT = solaire/zdist/zdist
	375
	376	IF (type_trac == 'repr') THEN
	377	#ifdef REPROBUS
	378	if(ok_SUNTIME) PSCT = solaireTIME/zdist/zdist
	379	print,'Constante solaire: ',PSCTzdist*zdist
	380	#endif
	381	END IF
	382
	383	DO j = 1, nb_gr
	384	iof = kdlon*(j-1)
	385	DO i = 1, kdlon
	386	zfract(i) = fract(iof+i)
[1989]	387	! zfract(i) = 1. !!!!!! essai MPL 19052010
[1687]	388	zrmu0(i) = rmu0(iof+i)
	389	PALBD(i,1) = alb1(iof+i)
	390	PALBD(i,2) = alb2(iof+i)
[1989]	391	!
	392	PALBD_NEW(i,1) = alb1(iof+i) !!!!! A REVOIR (MPL) PALBD_NEW en fonction bdes SW
	393	do kk=2,NSW
	394	PALBD_NEW(i,kk) = alb2(iof+i)
	395	enddo
[1687]	396	PALBP(i,1) = alb1(iof+i)
	397	PALBP(i,2) = alb2(iof+i)
[1989]	398	!
	399	PALBP_NEW(i,1) = alb1(iof+i) !!!!! A REVOIR (MPL) PALBP_NEW en fonction bdes SW
	400	do kk=2,NSW
	401	PALBP_NEW(i,kk) = alb2(iof+i)
	402	enddo
	403	PEMIS(i) = 1.0 !!!!! A REVOIR (MPL)
[1687]	404	PVIEW(i) = 1.66
	405	PPSOL(i) = paprs(iof+i,1)
	406	zx_alpha1 = (paprs(iof+i,1)-pplay(iof+i,2))/(pplay(iof+i,1)-pplay(iof+i,2))
	407	zx_alpha2 = 1.0 - zx_alpha1
	408	PTL(i,1) = t(iof+i,1) * zx_alpha1 + t(iof+i,2) * zx_alpha2
	409	PTL(i,KLEV+1) = t(iof+i,KLEV)
	410	PDT0(i) = tsol(iof+i) - PTL(i,1)
	411	ENDDO
	412	DO k = 2, kflev
	413	DO i = 1, kdlon
	414	PTL(i,k) = (t(iof+i,k)+t(iof+i,k-1))*0.5
	415	ENDDO
	416	ENDDO
	417	DO k = 1, kflev
	418	DO i = 1, kdlon
	419	PDP(i,k) = paprs(iof+i,k)-paprs(iof+i,k+1)
	420	PTAVE(i,k) = t(iof+i,k)
	421	PWV(i,k) = MAX (q(iof+i,k), 1.0e-12)
	422	PQS(i,k) = PWV(i,k)
	423	POZON(i,k, :) = wo(iof+i, k, :) * RG * dobson_u * 1e3 &
	424	/ (paprs(iof+i, k) - paprs(iof+i, k+1))
[1989]	425	! A activer pour CCMVAL on prend l'ozone impose (MPL 07042010)
	426	! POZON(i,k,:) = wo(i,k,:)
	427	! print *,'RADLWSW: POZON',k, POZON(i,k,1)
[1687]	428	PCLDLD(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
	429	PCLDLU(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
	430	PCLDSW(i,k) = cldfra(iof+i,k)
	431	PTAU(i,1,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! 1e-12 serait instable
	432	PTAU(i,2,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! pour 32-bit machines
	433	POMEGA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAU(i,1,k))
	434	POMEGA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i,2,k))
	435	PCG(i,1,k) = 0.865
	436	PCG(i,2,k) = 0.910
	437	!-
	438	! Introduced for aerosol indirect forcings.
	439	! The following values use the cloud optical thickness calculated from
	440	! present-day aerosol concentrations whereas the quantities without the
	441	! "A" at the end are for pre-industial (natural-only) aerosol concentrations
	442	!
	443	PTAUA(i,1,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! 1e-12 serait instable
	444	PTAUA(i,2,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! pour 32-bit machines
	445	POMEGAA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAUA(i,1,k))
	446	POMEGAA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i,2,k))
	447	ENDDO
	448	ENDDO
	449
	450	IF (type_trac == 'repr') THEN
	451	#ifdef REPROBUS
	452	ndimozon = size(wo, 3)
	453	CALL RAD_INTERACTIF(POZON,iof)
	454	#endif
	455	END IF
	456
	457	!
	458	DO k = 1, kflev+1
	459	DO i = 1, kdlon
	460	PPMB(i,k) = paprs(iof+i,k)/100.0
	461	ENDDO
	462	ENDDO
	463	!
[1989]	464	!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6
	465	DO kk = 1, 6
[1687]	466	DO k = 1, kflev
	467	DO i = 1, kdlon
[1989]	468	PAER(i,k,kk) = 1.0E-15 !!!!! A REVOIR (MPL)
[1687]	469	ENDDO
	470	ENDDO
	471	ENDDO
	472	DO k = 1, kflev
	473	DO i = 1, kdlon
	474	tauaero(i,k,:,1)=tau_aero(iof+i,k,:,1)
	475	pizaero(i,k,:,1)=piz_aero(iof+i,k,:,1)
	476	cgaero(i,k,:,1) =cg_aero(iof+i,k,:,1)
	477	tauaero(i,k,:,2)=tau_aero(iof+i,k,:,2)
	478	pizaero(i,k,:,2)=piz_aero(iof+i,k,:,2)
	479	cgaero(i,k,:,2) =cg_aero(iof+i,k,:,2)
	480	ENDDO
	481	ENDDO
	482
	483	!
	484	!===== iflag_rrtm ================================================
	485	!
[1989]	486	IF (iflag_rrtm == 0) THEN !!!! remettre 0 juste pour tester l'ancien rayt via rrtm
	487	!--- Mise a zero des tableaux output du rayonnement LW-AR4 ----------
	488	DO k = 1, kflev+1
	489	DO i = 1, kdlon
	490	! print *,'RADLWSW: boucle mise a zero i k',i,k
	491	ZFLUP(i,k)=0.
	492	ZFLDN(i,k)=0.
	493	ZFLUP0(i,k)=0.
	494	ZFLDN0(i,k)=0.
	495	ZLWFT0_i(i,k)=0.
	496	ZFLUCUP_i(i,k)=0.
	497	ZFLUCDWN_i(i,k)=0.
	498	ENDDO
	499	ENDDO
	500	DO k = 1, kflev
	501	DO i = 1, kdlon
	502	zcool(i,k)=0.
	503	zcool0(i,k)=0.
	504	ENDDO
	505	ENDDO
	506	DO i = 1, kdlon
	507	ztoplw(i)=0.
	508	zsollw(i)=0.
	509	ztoplw0(i)=0.
	510	zsollw0(i)=0.
	511	zsollwdown(i)=0.
	512	ENDDO
[1687]	513	! Old radiation scheme, used for AR4 runs
	514	! average day-night ozone for longwave
	515	CALL LW_LMDAR4(&
	516	PPMB, PDP,&
	517	PPSOL,PDT0,PEMIS,&
	518	PTL, PTAVE, PWV, POZON(:, :, 1), PAER,&
	519	PCLDLD,PCLDLU,&
	520	PVIEW,&
	521	zcool, zcool0,&
	522	ztoplw,zsollw,ztoplw0,zsollw0,&
	523	zsollwdown,&
	524	ZFLUP, ZFLDN, ZFLUP0,ZFLDN0)
[1989]	525	!----- Mise a zero des tableaux output du rayonnement SW-AR4
	526	DO k = 1, kflev+1
	527	DO i = 1, kdlon
	528	ZFSUP(i,k)=0.
	529	ZFSDN(i,k)=0.
	530	ZFSUP0(i,k)=0.
	531	ZFSDN0(i,k)=0.
	532	ZSWFT0_i(i,k)=0.
	533	ZFCUP_i(i,k)=0.
	534	ZFCDWN_i(i,k)=0.
	535	ENDDO
	536	ENDDO
	537	DO k = 1, kflev
	538	DO i = 1, kdlon
	539	zheat(i,k)=0.
	540	zheat0(i,k)=0.
	541	ENDDO
	542	ENDDO
	543	DO i = 1, kdlon
	544	zalbpla(i)=0.
	545	ztopsw(i)=0.
	546	zsolsw(i)=0.
	547	ztopsw0(i)=0.
	548	zsolsw0(i)=0.
	549	ztopswadaero(i)=0.
	550	zsolswadaero(i)=0.
	551	ztopswaiaero(i)=0.
	552	zsolswaiaero(i)=0.
	553	ENDDO
	554	! print *,'Avant SW_LMDAR4: PSCT zrmu0 zfract',PSCT, zrmu0, zfract
[1687]	555	! daylight ozone, if we have it, for short wave
	556	IF (.NOT. new_aod) THEN
	557	! use old version
	558	CALL SW_LMDAR4(PSCT, zrmu0, zfract,&
	559	PPMB, PDP, &
	560	PPSOL, PALBD, PALBP,&
	561	PTAVE, PWV, PQS, POZON(:, :, size(wo, 3)), PAER,&
	562	PCLDSW, PTAU, POMEGA, PCG,&
	563	zheat, zheat0,&
	564	zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,&
	565	ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,&
	566	tauaero(:,:,5,:), pizaero(:,:,5,:), cgaero(:,:,5,:),&
	567	PTAUA, POMEGAA,&
	568	ztopswadaero,zsolswadaero,&
	569	ztopswaiaero,zsolswaiaero,&
[1764]	570	ok_ade, ok_aie)
[1687]	571
	572	ELSE ! new_aod=T
	573	CALL SW_AEROAR4(PSCT, zrmu0, zfract,&
	574	PPMB, PDP,&
	575	PPSOL, PALBD, PALBP,&
	576	PTAVE, PWV, PQS, POZON(:, :, size(wo, 3)), PAER,&
	577	PCLDSW, PTAU, POMEGA, PCG,&
	578	zheat, zheat0,&
	579	zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,&
	580	ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,&
	581	tauaero, pizaero, cgaero, &
	582	PTAUA, POMEGAA,&
	583	ztopswadaero,zsolswadaero,&
	584	ztopswad0aero,zsolswad0aero,&
	585	ztopswaiaero,zsolswaiaero, &
	586	ztopsw_aero,ztopsw0_aero,&
	587	zsolsw_aero,zsolsw0_aero,&
	588	ztopswcf_aero,zsolswcf_aero, &
[1764]	589	ok_ade, ok_aie, flag_aerosol,flag_aerosol_strat)
[1687]	590	ENDIF
	591
[1989]	592
	593	DO i=1,kdlon
	594	DO k=1,kflev+1
	595	ZSWFT0_i(1:klon,k) = ZFSDN0(1:klon,k)-ZFSUP0(1:klon,k)
	596	ZLWFT0_i(1:klon,k)=-ZFLDN0(1:klon,k)-ZFLUP0(1:klon,k)
	597	! print *,'iof i k klon klev=',iof,i,k,klon,klev
	598	lwdn0 ( iof+i,k) = ZFLDN0 ( i,k)
	599	lwdn ( iof+i,k) = ZFLDN ( i,k)
	600	lwup0 ( iof+i,k) = ZFLUP0 ( i,k)
	601	lwup ( iof+i,k) = ZFLUP ( i,k)
	602	swdn0 ( iof+i,k) = ZFSDN0 ( i,k)
	603	swdn ( iof+i,k) = ZFSDN ( i,k)
	604	swup0 ( iof+i,k) = ZFSUP0 ( i,k)
	605	swup ( iof+i,k) = ZFSUP ( i,k)
	606	ENDDO
	607	ENDDO
	608	! print*,'SW_AR4 ZFSDN0 1 , klev:',ZFSDN0(1:klon,1),ZFSDN0(1:klon,klev)
	609	! print*,'SW_AR4 swdn0 1 , klev:',swdn0(1:klon,1),swdn0(1:klon,klev)
	610	! print*,'SW_AR4 ZFSUP0 1 , klev:',ZFSUP0(1:klon,1),ZFSUP0(1:klon,klev)
	611	! print*,'SW_AR4 swup0 1 , klev:',swup0(1:klon,1),swup0(1:klon,klev)
	612	! print*,'SW_AR4 ZFSDN 1 , klev:',ZFSDN(1:klon,1) ,ZFSDN(1:klon,klev)
	613	! print*,'SW_AR4 ZFSUP 1 , klev:',ZFSUP(1:klon,1) ,ZFSUP(1:klon,klev)
[1687]	614	ELSE
[1989]	615	#ifdef CPP_RRTM
	616	! if (prt_level.gt.10)write(lunout,*)'CPP_RRTM=.T.'
[1687]	617	!===== iflag_rrtm=1, on passe dans SW via RECMWFL ===============
	618
[1989]	619	DO k = 1, kflev+1
	620	DO i = 1, kdlon
	621	ZEMTD_i(i,k)=0.
	622	ZEMTU_i(i,k)=0.
	623	ZTRSO_i(i,k)=0.
	624	ZTH_i(i,k)=0.
	625	ZLWFT_i(i,k)=0.
	626	ZSWFT_i(i,k)=0.
	627	ZFLUX_i(i,1,k)=0.
	628	ZFLUX_i(i,2,k)=0.
	629	ZFLUC_i(i,1,k)=0.
	630	ZFLUC_i(i,2,k)=0.
	631	ZFSDWN_i(i,k)=0.
	632	ZFCDWN_i(i,k)=0.
	633	ZFSUP_i(i,k)=0.
	634	ZFCUP_i(i,k)=0.
	635	ENDDO
	636	ENDDO
[2003]	637	!
	638	!--OB
	639	!--aerosol TOT - anthropogenic+natural
	640	!--aerosol NAT - natural only
	641	!
[1989]	642	DO i = 1, kdlon
	643	DO k = 1, kflev
	644	DO kk=1, NSW
[2003]	645	!
	646	PTAU_TOT(i,kflev+1-k,kk)=tau_aero_rrtm(i,k,2,kk)
	647	PPIZA_TOT(i,kflev+1-k,kk)=piz_aero_rrtm(i,k,2,kk)
	648	PCGA_TOT(i,kflev+1-k,kk)=cg_aero_rrtm(i,k,2,kk)
	649	!
	650	PTAU_NAT(i,kflev+1-k,kk)=tau_aero_rrtm(i,k,1,kk)
	651	PPIZA_NAT(i,kflev+1-k,kk)=piz_aero_rrtm(i,k,1,kk)
	652	PCGA_NAT(i,kflev+1-k,kk)=cg_aero_rrtm(i,k,1,kk)
	653	!
[1989]	654	ENDDO
	655	ENDDO
	656	ENDDO
[2003]	657	!-end OB
[1989]	658	!
	659	!
	660	DO i = 1, kdlon
	661	ZCTRSO(i,1)=0.
	662	ZCTRSO(i,2)=0.
	663	ZCEMTR(i,1)=0.
	664	ZCEMTR(i,2)=0.
	665	ZTRSOD(i)=0.
	666	ZLWFC(i,1)=0.
	667	ZLWFC(i,2)=0.
	668	ZSWFC(i,1)=0.
	669	ZSWFC(i,2)=0.
	670	PFSDNN(i)=0.
	671	PFSDNV(i)=0.
	672	DO kk = 1, NSW
	673	PSFSWDIR(i,kk)=0.
	674	PSFSWDIF(i,kk)=0.
	675	ENDDO
	676	ENDDO
	677	!----- Fin des mises a zero des tableaux output de RECMWF -------------------
	678	! GEMU(1:klon)=sin(rlatd(1:klon))
	679	! On met les donnees dans l'ordre des niveaux arpege
	680	paprs_i(:,1)=paprs(:,klev+1)
	681	do k=1,klev
	682	paprs_i(1:klon,k+1) =paprs(1:klon,klev+1-k)
	683	pplay_i(1:klon,k) =pplay(1:klon,klev+1-k)
	684	cldfra_i(1:klon,k) =cldfra(1:klon,klev+1-k)
	685	PDP_i(1:klon,k) =PDP(1:klon,klev+1-k)
	686	t_i(1:klon,k) =t(1:klon,klev+1-k)
	687	q_i(1:klon,k) =q(1:klon,klev+1-k)
	688	qsat_i(1:klon,k) =qsat(1:klon,klev+1-k)
	689	flwc_i(1:klon,k) =flwc(1:klon,klev+1-k)
	690	fiwc_i(1:klon,k) =fiwc(1:klon,klev+1-k)
	691	ref_liq_i(1:klon,k) =ref_liq(1:klon,klev+1-k)
	692	ref_ice_i(1:klon,k) =ref_ice(1:klon,klev+1-k)
[2003]	693	!-OB
	694	ref_liq_pi_i(1:klon,k) =ref_liq_pi(1:klon,klev+1-k)
	695	ref_ice_pi_i(1:klon,k) =ref_ice_pi(1:klon,klev+1-k)
[1989]	696	enddo
	697	do k=1,kflev
	698	POZON_i(1:klon,k,:)=POZON(1:klon,kflev+1-k,:)
	699	!!! POZON_i(1:klon,k)=POZON(1:klon,k) !!! on laisse 1=sol et klev=top
	700	! print *,'Juste avant RECMWFL: k tsol temp',k,tsol,t(1,k)
	701	!!!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6
	702	do i=1,6
	703	PAER_i(1:klon,k,i)=PAER(1:klon,kflev+1-k,i)
	704	enddo
	705	enddo
	706	! print *,'RADLWSW: avant RECMWFL, RI0,rmu0=',solaire,rmu0
	707
	708	! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	709	! La version ARPEGE1D utilise differentes valeurs de la constante
	710	! solaire suivant le rayonnement utilise.
	711	! A controler ...
	712	! SOLAR FLUX AT THE TOP (/YOMPHY3/)
	713	! introduce season correction
	714	!--------------------------------------
	715	! RII0 = RIP0
	716	! IF(LRAYFM)
	717	! RII0 = RIP0M ! =rip0m if Morcrette non-each time step call.
	718	! IF(LRAYFM15)
	719	! RII0 = RIP0M15 ! =rip0m if Morcrette non-each time step call.
	720	RII0=solaire/zdist/zdist
	721	print*,'+++ radlwsw: solaire ,RII0',solaire,RII0
	722	! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	723	! Ancien appel a RECMWF (celui du cy25)
	724	! CALL RECMWF (ist , iend, klon , ktdia , klev , kmode ,
	725	! s PALBD , PALBP , paprs_i , pplay_i , RCO2 , cldfra_i,
	726	! s POZON_i , PAER_i , PDP_i , PEMIS , GEMU , rmu0,
	727	! s q_i , qsat_i , fiwc_i , flwc_i , zmasq , t_i ,tsol,
	728	! s ZEMTD_i , ZEMTU_i , ZTRSO_i ,
	729	! s ZTH_i , ZCTRSO , ZCEMTR , ZTRSOD ,
	730	! s ZLWFC , ZLWFT_i , ZSWFC , ZSWFT_i ,
	731	! s ZFLUX_i , ZFLUC_i , ZFSDWN_i, ZFSUP_i , ZFCDWN_i,ZFCUP_i)
	732	! s 'RECMWF ')
	733	!
	734	if(lldebug) then
	735	CALL writefield_phy('paprs_i',paprs_i,klev+1)
	736	CALL writefield_phy('pplay_i',pplay_i,klev)
	737	CALL writefield_phy('cldfra_i',cldfra_i,klev)
	738	CALL writefield_phy('pozon_i',POZON_i,klev)
	739	CALL writefield_phy('paer_i',PAER_i,klev)
	740	CALL writefield_phy('pdp_i',PDP_i,klev)
	741	CALL writefield_phy('q_i',q_i,klev)
	742	CALL writefield_phy('qsat_i',qsat_i,klev)
	743	CALL writefield_phy('fiwc_i',fiwc_i,klev)
	744	CALL writefield_phy('flwc_i',flwc_i,klev)
	745	CALL writefield_phy('t_i',t_i,klev)
	746	CALL writefield_phy('palbd_new',PALBD_NEW,NSW)
	747	CALL writefield_phy('palbp_new',PALBP_NEW,NSW)
	748	endif
	749
	750	! Nouvel appel a RECMWF (celui du cy32t0)
[2003]	751	CALL RECMWF_AERO (ist , iend, klon , ktdia , klev , kmode ,&
[1989]	752	PALBD_NEW,PALBP_NEW, paprs_i , pplay_i , RCO2 , cldfra_i,&
	753	POZON_i , PAER_i , PDP_i , PEMIS , rmu0 ,&
	754	q_i , qsat_i , fiwc_i , flwc_i , zmasq , t_i ,tsol,&
	755	ref_liq_i, ref_ice_i, &
[2003]	756	ref_liq_pi_i, ref_ice_pi_i, & ! rajoute par OB pour diagnostiquer effet indirect
[1989]	757	ZEMTD_i , ZEMTU_i , ZTRSO_i ,&
	758	ZTH_i , ZCTRSO , ZCEMTR , ZTRSOD ,&
	759	ZLWFC , ZLWFT_i , ZSWFC , ZSWFT_i ,&
	760	PSFSWDIR , PSFSWDIF, PFSDNN , PFSDNV ,&
[2003]	761	PPIZA_TOT, PCGA_TOT,PTAU_TOT,&
	762	PPIZA_NAT, PCGA_NAT,PTAU_NAT, & ! rajoute par OB pour diagnostiquer effet direct
	763	ZFLUX_i , ZFLUC_i ,&
	764	ZFSDWN_i , ZFSUP_i , ZFCDWN_i, ZFCUP_i,&
	765	ZTOPSWADAERO,ZSOLSWADAERO,& ! rajoute par OB pour diagnostics
	766	ZTOPSWAD0AERO,ZSOLSWAD0AERO,&
	767	ZTOPSWAIAERO,ZSOLSWAIAERO, &
	768	ZTOPSWCF_AERO,ZSOLSWCF_AERO, &
	769	ok_ade, ok_aie, flag_aerosol,flag_aerosol_strat) ! flags aerosols
[1989]	770
	771	print *,'RADLWSW: apres RECMWF'
	772	if(lldebug) then
	773	CALL writefield_phy('zemtd_i',ZEMTD_i,klev+1)
	774	CALL writefield_phy('zemtu_i',ZEMTU_i,klev+1)
	775	CALL writefield_phy('ztrso_i',ZTRSO_i,klev+1)
	776	CALL writefield_phy('zth_i',ZTH_i,klev+1)
	777	CALL writefield_phy('zctrso',ZCTRSO,2)
	778	CALL writefield_phy('zcemtr',ZCEMTR,2)
	779	CALL writefield_phy('ztrsod',ZTRSOD,1)
	780	CALL writefield_phy('zlwfc',ZLWFC,2)
	781	CALL writefield_phy('zlwft_i',ZLWFT_i,klev+1)
	782	CALL writefield_phy('zswfc',ZSWFC,2)
	783	CALL writefield_phy('zswft_i',ZSWFT_i,klev+1)
	784	CALL writefield_phy('psfswdir',PSFSWDIR,6)
	785	CALL writefield_phy('psfswdif',PSFSWDIF,6)
	786	CALL writefield_phy('pfsdnn',PFSDNN,1)
	787	CALL writefield_phy('pfsdnv',PFSDNV,1)
[2003]	788	CALL writefield_phy('ppiza_dst',PPIZA_TOT,klev)
	789	CALL writefield_phy('pcga_dst',PCGA_TOT,klev)
	790	CALL writefield_phy('ptaurel_dst',PTAU_TOT,klev)
[1989]	791	CALL writefield_phy('zflux_i',ZFLUX_i,klev+1)
	792	CALL writefield_phy('zfluc_i',ZFLUC_i,klev+1)
	793	CALL writefield_phy('zfsdwn_i',ZFSDWN_i,klev+1)
	794	CALL writefield_phy('zfsup_i',ZFSUP_i,klev+1)
	795	CALL writefield_phy('zfcdwn_i',ZFCDWN_i,klev+1)
	796	CALL writefield_phy('zfcup_i',ZFCUP_i,klev+1)
	797	endif
	798	! --------- output RECMWFL
	799	! ZEMTD (KPROMA,KLEV+1) ; TOTAL DOWNWARD LONGWAVE EMISSIVITY
	800	! ZEMTU (KPROMA,KLEV+1) ; TOTAL UPWARD LONGWAVE EMISSIVITY
	801	! ZTRSO (KPROMA,KLEV+1) ; TOTAL SHORTWAVE TRANSMISSIVITY
	802	! ZTH (KPROMA,KLEV+1) ; HALF LEVEL TEMPERATURE
	803	! ZCTRSO (KPROMA,2) ; CLEAR-SKY SHORTWAVE TRANSMISSIVITY
	804	! ZCEMTR (KPROMA,2) ; CLEAR-SKY NET LONGWAVE EMISSIVITY
	805	! ZTRSOD (KPROMA) ; TOTAL-SKY SURFACE SW TRANSMISSITY
	806	! ZLWFC (KPROMA,2) ; CLEAR-SKY LONGWAVE FLUXES
	807	! ZLWFT (KPROMA,KLEV+1) ; TOTAL-SKY LONGWAVE FLUXES
	808	! ZSWFC (KPROMA,2) ; CLEAR-SKY SHORTWAVE FLUXES
	809	! ZSWFT (KPROMA,KLEV+1) ; TOTAL-SKY SHORTWAVE FLUXES
[2003]	810	! PPIZA_TOT (KPROMA,KLEV,NSW); Single scattering albedo of total aerosols
	811	! PCGA_TOT (KPROMA,KLEV,NSW); Assymetry factor for total aerosols
	812	! PTAU_TOT (KPROMA,KLEV,NSW); Optical depth of total aerosols
	813	! PPIZA_NAT (KPROMA,KLEV,NSW); Single scattering albedo of natural aerosols
	814	! PCGA_NAT (KPROMA,KLEV,NSW); Assymetry factor for natural aerosols
	815	! PTAU_NAT (KPROMA,KLEV,NSW); Optical depth of natiral aerosols
[1989]	816	! PSFSWDIR (KPROMA,NSW) ;
	817	! PSFSWDIF (KPROMA,NSW) ;
	818	! PFSDNN (KPROMA) ;
	819	! PFSDNV (KPROMA) ;
	820	! ---------
	821	! ---------
	822	! On retablit l'ordre des niveaux lmd pour les tableaux de sortie
	823	! D autre part, on multiplie les resultats SW par fract pour etre coherent
	824	! avec l ancien rayonnement AR4. Si nuit, fract=0 donc pas de
	825	! rayonnement SW. (MPL 260609)
	826	DO k=0,klev
	827	DO i=1,klon
	828	ZEMTD(i,k+1) = ZEMTD_i(i,k+1)
	829	ZEMTU(i,k+1) = ZEMTU_i(i,k+1)
	830	ZTRSO(i,k+1) = ZTRSO_i(i,k+1)
	831	ZTH(i,k+1) = ZTH_i(i,k+1)
	832	! ZLWFT(i,k+1) = ZLWFT_i(i,klev+1-k)
	833	! ZSWFT(i,k+1) = ZSWFT_i(i,klev+1-k)
	834	ZFLUP(i,k+1) = ZFLUX_i(i,1,k+1)
	835	ZFLDN(i,k+1) = ZFLUX_i(i,2,k+1)
	836	ZFLUP0(i,k+1) = ZFLUC_i(i,1,k+1)
	837	ZFLDN0(i,k+1) = ZFLUC_i(i,2,k+1)
	838	ZFSDN(i,k+1) = ZFSDWN_i(i,k+1)*fract(i)
	839	ZFSDN0(i,k+1) = ZFCDWN_i(i,k+1)*fract(i)
	840	ZFSUP (i,k+1) = ZFSUP_i(i,k+1)*fract(i)
	841	ZFSUP0(i,k+1) = ZFCUP_i(i,k+1)*fract(i)
	842	! Nouveau calcul car visiblement ZSWFT et ZSWFC sont nuls dans RRTM cy32
	843	! en sortie de radlsw.F90 - MPL 7.01.09
	844	ZSWFT(i,k+1) = (ZFSDWN_i(i,k+1)-ZFSUP_i(i,k+1))*fract(i)
	845	ZSWFT0_i(i,k+1) = (ZFCDWN_i(i,k+1)-ZFCUP_i(i,k+1))*fract(i)
	846	! WRITE(*,'("FSDN FSUP FCDN FCUP: ",4E12.5)') ZFSDWN_i(i,k+1),&
	847	! ZFSUP_i(i,k+1),ZFCDWN_i(i,k+1),ZFCUP_i(i,k+1)
	848	ZLWFT(i,k+1) =-ZFLUX_i(i,2,k+1)-ZFLUX_i(i,1,k+1)
	849	ZLWFT0_i(i,k+1)=-ZFLUC_i(i,2,k+1)-ZFLUC_i(i,1,k+1)
	850	! print *,'FLUX2 FLUX1 FLUC2 FLUC1',ZFLUX_i(i,2,k+1),&
	851	! & ZFLUX_i(i,1,k+1),ZFLUC_i(i,2,k+1),ZFLUC_i(i,1,k+1)
	852	ENDDO
	853	ENDDO
	854
[2003]	855	!--ajout OB
	856	ZTOPSWADAERO(:) =ZTOPSWADAERO(:) *fract(:)
	857	ZSOLSWADAERO(:) =ZSOLSWADAERO(:) *fract(:)
	858	ZTOPSWAD0AERO(:)=ZTOPSWAD0AERO(:)*fract(:)
	859	ZSOLSWAD0AERO(:)=ZSOLSWAD0AERO(:)*fract(:)
	860	ZTOPSWAIAERO(:) =ZTOPSWAIAERO(:) *fract(:)
	861	ZSOLSWAIAERO(:) =ZSOLSWAIAERO(:) *fract(:)
	862	ZTOPSWCF_AERO(:,1)=ZTOPSWCF_AERO(:,1)*fract(:)
	863	ZTOPSWCF_AERO(:,2)=ZTOPSWCF_AERO(:,2)*fract(:)
	864	ZTOPSWCF_AERO(:,3)=ZTOPSWCF_AERO(:,3)*fract(:)
	865	ZSOLSWCF_AERO(:,1)=ZSOLSWCF_AERO(:,1)*fract(:)
	866	ZSOLSWCF_AERO(:,2)=ZSOLSWCF_AERO(:,2)*fract(:)
	867	ZSOLSWCF_AERO(:,3)=ZSOLSWCF_AERO(:,3)*fract(:)
	868
[1989]	869	! print*,'SW_RRTM ZFSDN0 1 , klev:',ZFSDN0(1:klon,1),ZFSDN0(1:klon,klev)
	870	! print*,'SW_RRTM ZFSUP0 1 , klev:',ZFSUP0(1:klon,1),ZFSUP0(1:klon,klev)
	871	! print*,'SW_RRTM ZFSDN 1 , klev:',ZFSDN(1:klon,1),ZFSDN(1:klon,klev)
	872	! print*,'SW_RRTM ZFSUP 1 , klev:',ZFSUP(1:klon,1),ZFSUP(1:klon,klev)
	873	! print*,'OK1'
	874	! ---------
	875	! ---------
	876	! On renseigne les champs LMDz, pour avoir la meme chose qu'en sortie de
	877	! LW_LMDAR4 et SW_LMDAR4
	878	DO i = 1, kdlon
	879	zsolsw(i) = ZSWFT(i,1)
	880	zsolsw0(i) = ZSWFT0_i(i,1)
	881	! zsolsw0(i) = ZFSDN0(i,1) -ZFSUP0(i,1)
	882	ztopsw(i) = ZSWFT(i,klev+1)
	883	ztopsw0(i) = ZSWFT0_i(i,klev+1)
	884	! ztopsw0(i) = ZFSDN0(i,klev+1)-ZFSUP0(i,klev+1)
	885	!
	886	! zsollw(i) = ZFLDN(i,1) -ZFLUP(i,1)
	887	! zsollw0(i) = ZFLDN0(i,1) -ZFLUP0(i,1)
	888	! ztoplw(i) = ZFLDN(i,klev+1) -ZFLUP(i,klev+1)
	889	! ztoplw0(i) = ZFLDN0(i,klev+1)-ZFLUP0(i,klev+1)
	890	zsollw(i) = ZLWFT(i,1)
	891	zsollw0(i) = ZLWFT0_i(i,1)
	892	ztoplw(i) = ZLWFT(i,klev+1)*(-1)
	893	ztoplw0(i) = ZLWFT0_i(i,klev+1)*(-1)
	894	!
	895	IF (fract(i) == 0.) THEN
	896	!!!!! A REVOIR MPL (20090630) ca n a pas de sens quand fract=0
	897	! pas plus que dans le sw_AR4
	898	zalbpla(i) = 1.0e+39
	899	ELSE
	900	zalbpla(i) = ZFSUP(i,klev+1)/ZFSDN(i,klev+1)
	901	ENDIF
	902	zsollwdown(i)= ZFLDN(i,1)
	903	ENDDO
	904	print*,'OK2'
	905
	906	! extrait de SW_AR4
	907	! DO k = 1, KFLEV
	908	! kpl1 = k+1
	909	! DO i = 1, KDLON
	910	! PHEAT(i,k) = -(ZFSUP(i,kpl1)-ZFSUP(i,k)) -(ZFSDN(i,k)-ZFSDN(i,kpl1))
	911	! PHEAT(i,k) = PHEAT(i,k) * RDAY*RG/RCPD / PDP(i,k)
	912	! ZLWFT(klon,k),ZSWFT
	913
	914	do k=1,kflev
	915	do i=1,kdlon
	916	zheat(i,k)=(ZSWFT(i,k+1)-ZSWFT(i,k))RDAYRG/RCPD/PDP(i,k)
	917	zheat0(i,k)=(ZSWFT0_i(i,k+1)-ZSWFT0_i(i,k))RDAYRG/RCPD/PDP(i,k)
	918	zcool(i,k)=(ZLWFT(i,k)-ZLWFT(i,k+1))RDAYRG/RCPD/PDP(i,k)
	919	zcool0(i,k)=(ZLWFT0_i(i,k)-ZLWFT0_i(i,k+1))RDAYRG/RCPD/PDP(i,k)
	920	! print *,'heat cool heat0 cool0 ',zheat(i,k),zcool(i,k),zheat0(i,k),zcool0(i,k)
	921	! ZFLUCUP_i(i,k)=ZFLUC_i(i,1,k)
	922	! ZFLUCDWN_i(i,k)=ZFLUC_i(i,2,k)
	923	enddo
	924	enddo
	925	#else
[1991]	926	abort_message="You should compile with -rrtm if running with iflag_rrtm=1"
[1989]	927	call abort_gcm(modname, abort_message, 1)
	928	#endif
[1687]	929	ENDIF ! iflag_rrtm
	930	!======================================================================
	931
	932	DO i = 1, kdlon
	933	radsol(iof+i) = zsolsw(i) + zsollw(i)
	934	topsw(iof+i) = ztopsw(i)
	935	toplw(iof+i) = ztoplw(i)
	936	solsw(iof+i) = zsolsw(i)
	937	sollw(iof+i) = zsollw(i)
	938	sollwdown(iof+i) = zsollwdown(i)
	939	DO k = 1, kflev+1
	940	lwdn0 ( iof+i,k) = ZFLDN0 ( i,k)
	941	lwdn ( iof+i,k) = ZFLDN ( i,k)
	942	lwup0 ( iof+i,k) = ZFLUP0 ( i,k)
	943	lwup ( iof+i,k) = ZFLUP ( i,k)
	944	ENDDO
	945	topsw0(iof+i) = ztopsw0(i)
	946	toplw0(iof+i) = ztoplw0(i)
	947	solsw0(iof+i) = zsolsw0(i)
	948	sollw0(iof+i) = zsollw0(i)
	949	albpla(iof+i) = zalbpla(i)
	950
	951	DO k = 1, kflev+1
	952	swdn0 ( iof+i,k) = ZFSDN0 ( i,k)
	953	swdn ( iof+i,k) = ZFSDN ( i,k)
	954	swup0 ( iof+i,k) = ZFSUP0 ( i,k)
	955	swup ( iof+i,k) = ZFSUP ( i,k)
	956	ENDDO
	957	ENDDO
	958	!-transform the aerosol forcings, if they have
	959	! to be calculated
	960	IF (ok_ade) THEN
	961	DO i = 1, kdlon
	962	topswad_aero(iof+i) = ztopswadaero(i)
	963	topswad0_aero(iof+i) = ztopswad0aero(i)
	964	solswad_aero(iof+i) = zsolswadaero(i)
	965	solswad0_aero(iof+i) = zsolswad0aero(i)
	966	! MS the following lines seem to be wrong, why is iof on right hand side???
	967	! topsw_aero(iof+i,:) = ztopsw_aero(iof+i,:)
	968	! topsw0_aero(iof+i,:) = ztopsw0_aero(iof+i,:)
	969	! solsw_aero(iof+i,:) = zsolsw_aero(iof+i,:)
	970	! solsw0_aero(iof+i,:) = zsolsw0_aero(iof+i,:)
	971	topsw_aero(iof+i,:) = ztopsw_aero(i,:)
	972	topsw0_aero(iof+i,:) = ztopsw0_aero(i,:)
	973	solsw_aero(iof+i,:) = zsolsw_aero(i,:)
	974	solsw0_aero(iof+i,:) = zsolsw0_aero(i,:)
	975	topswcf_aero(iof+i,:) = ztopswcf_aero(i,:)
	976	solswcf_aero(iof+i,:) = zsolswcf_aero(i,:)
	977	ENDDO
	978	ELSE
	979	DO i = 1, kdlon
	980	topswad_aero(iof+i) = 0.0
	981	solswad_aero(iof+i) = 0.0
	982	topswad0_aero(iof+i) = 0.0
	983	solswad0_aero(iof+i) = 0.0
	984	topsw_aero(iof+i,:) = 0.
	985	topsw0_aero(iof+i,:) =0.
	986	solsw_aero(iof+i,:) = 0.
	987	solsw0_aero(iof+i,:) = 0.
	988	ENDDO
	989	ENDIF
	990	IF (ok_aie) THEN
	991	DO i = 1, kdlon
	992	topswai_aero(iof+i) = ztopswaiaero(i)
	993	solswai_aero(iof+i) = zsolswaiaero(i)
	994	ENDDO
	995	ELSE
	996	DO i = 1, kdlon
	997	topswai_aero(iof+i) = 0.0
	998	solswai_aero(iof+i) = 0.0
	999	ENDDO
	1000	ENDIF
	1001	DO k = 1, kflev
	1002	DO i = 1, kdlon
	1003	! scale factor to take into account the difference between
	1004	! dry air and watter vapour scpecifi! heat capacity
	1005	zznormcp=1.0+RVTMP2*PWV(i,k)
	1006	heat(iof+i,k) = zheat(i,k)/zznormcp
	1007	cool(iof+i,k) = zcool(i,k)/zznormcp
	1008	heat0(iof+i,k) = zheat0(i,k)/zznormcp
	1009	cool0(iof+i,k) = zcool0(i,k)/zznormcp
	1010	ENDDO
	1011	ENDDO
	1012
	1013	ENDDO ! j = 1, nb_gr
	1014
	1015	END SUBROUTINE radlwsw
	1016
	1017	end module radlwsw_m

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

Download in other formats:

Original Format