phylmd/Radlwsw/radlwsw.f

module radlwsw_m

  IMPLICIT none

contains

  SUBROUTINE radlwsw(dist, mu0, fract, paprs, play, tsol, albedo, t, q, wo, &
       cldfra, cldemi, cldtaupd, heat, heat0, cool, cool0, radsol, albpla, &
       topsw, toplw, solsw, sollw, sollwdown, topsw0, toplw0, solsw0, sollw0, &
       lwdn0, lwdn, lwup0, lwup, swdn0, swdn, swup0, swup, ok_ade, topswad, &
       solswad)

    ! From LMDZ4/libf/phylmd/radlwsw.F, version 1.4, 2005/06/06 13:16:33
    ! Author: Z. X. Li (LMD/CNRS) 
    ! Date: 1996/07/19

    ! Objet : interface entre le modèle et les rayonnements solaire et
    ! infrarouge

    ! ATTENTION: swai and swad have to be interpreted in the following manner:

    ! not ok_ade
    ! both are zero

    ! ok_ade
    ! aerosol direct forcing is F_{AD} = topsw - topswad
    ! indirect is zero

    USE clesphys, ONLY: solaire
    USE dimphy, ONLY: klev, klon
    use lw_m, only: lw
    USE raddim, ONLY: kdlon
    USE suphec_m, ONLY: rg
    use sw_m, only: sw
    USE yoethf_m, ONLY: rvtmp2

    real, intent(in):: dist ! distance astronomique terre-soleil
    real, intent(in):: mu0(klon) ! cosinus de l'angle zenithal
    real, intent(in):: fract(klon) ! duree d'ensoleillement normalisee
    real, intent(in):: paprs(klon, klev+1) ! pression a inter-couche (Pa)
    real, intent(in):: play(klon, klev) ! pression au milieu de couche (Pa)
    real, intent(in):: tsol(klon) ! temperature du sol (en K)
    real, intent(in):: albedo(klon) ! albedo du sol (entre 0 et 1)
    real, intent(in):: t(klon, klev) ! temperature (K)
    real, intent(in):: q(klon, klev) ! vapeur d'eau (en kg/kg)

    real, intent(in):: wo(klon, klev)
    ! column-density of ozone in a layer, in kilo-Dobsons

    real, intent(in):: cldfra(klon, klev) ! fraction nuageuse (entre 0 et 1)

    real, intent(in):: cldemi(klon, klev)
    ! emissivite des nuages dans l'IR (entre 0 et 1)

    real, intent(in):: cldtaupd(klon, klev)
    ! epaisseur optique des nuages dans le visible (present-day value)

    real, intent(out):: heat(klon, klev)
    ! échauffement atmosphérique (visible) (K/jour)

    real, intent(out):: heat0(klon, klev) ! chauffage solaire ciel clair
    real, intent(out):: cool(klon, klev) ! refroidissement dans l'IR (K/jour)

    real, intent(out):: cool0(klon, klev)
    ! refroidissement infrarouge ciel clair

    real, intent(out):: radsol(klon)
    ! bilan radiatif net au sol (W/m**2) (+ vers le bas)

    real, intent(out):: albpla(klon) ! albedo planetaire (entre 0 et 1)
    real, intent(out):: topsw(klon) ! flux solaire net au sommet de l'atm.

    real, intent(out):: toplw(klon)
    ! rayonnement infrarouge montant au sommet de l'atmosphère

    real, intent(out):: solsw(klon) ! flux solaire net à la surface

    real, intent(out):: sollw(klon)
    ! rayonnement infrarouge montant à la surface

    real, intent(out):: sollwdown(klon)
    real, intent(out):: topsw0(klon)
    real, intent(out):: toplw0(klon)
    real, intent(out):: solsw0(klon), sollw0(klon)
    REAL, intent(out):: lwdn0(klon, klev+1), lwdn(klon, klev+1)
    REAL, intent(out):: lwup0(klon, klev+1), lwup(klon, klev+1)
    REAL, intent(out):: swdn0(klon, klev+1), swdn(klon, klev+1)
    REAL, intent(out):: swup0(klon, klev+1), swup(klon, klev+1)

    logical, intent(in):: ok_ade ! apply the Aerosol Direct Effect

    real, intent(out):: topswad(klon), solswad(klon)
    ! aerosol direct forcing at TOA and surface
    ! rayonnement solaire net absorb\'e

    ! Local:

    DOUBLE PRECISION ZFSUP(KDLON, KLEV+1)
    DOUBLE PRECISION ZFSDN(KDLON, KLEV+1)
    DOUBLE PRECISION ZFSUP0(KDLON, KLEV+1)
    DOUBLE PRECISION ZFSDN0(KDLON, KLEV+1)

    DOUBLE PRECISION ZFLUP(KDLON, KLEV+1)
    DOUBLE PRECISION ZFLDN(KDLON, KLEV+1)
    DOUBLE PRECISION ZFLUP0(KDLON, KLEV+1)
    DOUBLE PRECISION ZFLDN0(KDLON, KLEV+1)

    DOUBLE PRECISION zx_alpha1, zx_alpha2
    INTEGER k, kk, i, iof, nb_gr
    DOUBLE PRECISION PSCT

    DOUBLE PRECISION PALBD(kdlon, 2), PALBP(kdlon, 2)
    DOUBLE PRECISION PEMIS(kdlon), PDT0(kdlon), PVIEW(kdlon)
    DOUBLE PRECISION PPSOL(kdlon), PDP(kdlon, klev)
    DOUBLE PRECISION PTL(kdlon, klev+1), PPMB(kdlon, klev+1)
    DOUBLE PRECISION PTAVE(kdlon, klev)
    DOUBLE PRECISION PWV(kdlon, klev), PQS(kdlon, klev)
    DOUBLE PRECISION POZON(kdlon, klev) ! mass fraction of ozone
    DOUBLE PRECISION PAER(kdlon, klev, 5) ! AEROSOLS' OPTICAL THICKNESS
    DOUBLE PRECISION PCLDLD(kdlon, klev)
    DOUBLE PRECISION PCLDLU(kdlon, klev)
    DOUBLE PRECISION PCLDSW(kdlon, klev)
    DOUBLE PRECISION PTAU(kdlon, 2, klev)
    DOUBLE PRECISION POMEGA(kdlon, 2, klev)
    DOUBLE PRECISION PCG(kdlon, 2, klev)

    DOUBLE PRECISION zfract(kdlon), zrmu0(kdlon)

    DOUBLE PRECISION zheat(kdlon, klev), zcool(kdlon, klev)
    DOUBLE PRECISION zheat0(kdlon, klev), zcool0(kdlon, klev)
    DOUBLE PRECISION ztopsw(kdlon), ztoplw(kdlon)
    DOUBLE PRECISION zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon)
    DOUBLE PRECISION zsollwdown(kdlon)

    DOUBLE PRECISION ztopsw0(kdlon), ztoplw0(kdlon)
    DOUBLE PRECISION zsolsw0(kdlon), zsollw0(kdlon)
    DOUBLE PRECISION zznormcp

    !jq the following quantities are needed for the aerosol radiative forcings

    DOUBLE PRECISION PTAUA(kdlon, 2, klev)
    ! present-day value of cloud opt thickness (PTAU is pre-industrial
    ! value), local use

    DOUBLE PRECISION POMEGAA(kdlon, 2, klev) ! dito for single scatt albedo

    DOUBLE PRECISION ztopswad(kdlon), zsolswad(kdlon) 
    ! Aerosol direct forcing at TOAand surface

    DOUBLE PRECISION ztopswai(kdlon), zsolswai(kdlon) ! dito, indirect
    real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2

    !----------------------------------------------------------------------

    nb_gr = klon / kdlon
    IF (nb_gr * kdlon /= klon) THEN
       PRINT *, "kdlon mauvais :", klon, kdlon, nb_gr
       stop 1
    ENDIF

    heat = 0.
    cool = 0.
    heat0 = 0.
    cool0 = 0.
    PSCT = solaire / dist**2

    loop_iof: DO iof = 0, klon - kdlon, kdlon
       DO i = 1, kdlon
          zfract(i) = fract(iof+i)
          zrmu0(i) = mu0(iof+i)
          PALBD(i, 1) = albedo(iof+i)
          PALBD(i, 2) = albedo(iof+i)
          PALBP(i, 1) = albedo(iof+i)
          PALBP(i, 2) = albedo(iof+i)
          ! cf. JLD pour etre en accord avec ORCHIDEE il faut mettre
          ! PEMIS(i) = 0.96
          PEMIS(i) = 1. 
          PVIEW(i) = 1.66
          PPSOL(i) = paprs(iof+i, 1)
          zx_alpha1 = (paprs(iof+i, 1)-play(iof+i, 2))  &
               / (play(iof+i, 1)-play(iof+i, 2))
          zx_alpha2 = 1. - zx_alpha1
          PTL(i, 1) = t(iof+i, 1) * zx_alpha1 + t(iof+i, 2) * zx_alpha2
          PTL(i, klev+1) = t(iof+i, klev)
          PDT0(i) = tsol(iof+i) - PTL(i, 1)
       ENDDO
       DO k = 2, klev
          DO i = 1, kdlon
             PTL(i, k) = (t(iof+i, k)+t(iof+i, k-1))*0.5
          ENDDO
       ENDDO
       DO k = 1, klev
          DO i = 1, kdlon
             PDP(i, k) = paprs(iof+i, k)-paprs(iof+i, k+1)
             PTAVE(i, k) = t(iof+i, k)
             PWV(i, k) = MAX (q(iof+i, k), 1e-12)
             PQS(i, k) = PWV(i, k)
             POZON(i, k) = wo(iof+i, k) * RG * dobson_u * 1e3 &
                  / (paprs(iof+i, k) - paprs(iof+i, k+1))
             PCLDLD(i, k) = cldfra(iof+i, k)*cldemi(iof+i, k)
             PCLDLU(i, k) = cldfra(iof+i, k)*cldemi(iof+i, k)
             PCLDSW(i, k) = cldfra(iof+i, k)
             PTAU(i, 1, k) = MAX(cldtaupd(iof+i, k), 1e-05)
             ! (1e-12 serait instable)
             PTAU(i, 2, k) = MAX(cldtaupd(iof+i, k), 1e-05)
             ! (pour 32-bit machines)
             POMEGA(i, 1, k) = 0.9999 - 5e-04 * EXP(-0.5 * PTAU(i, 1, k))
             POMEGA(i, 2, k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i, 2, k))
             PCG(i, 1, k) = 0.865
             PCG(i, 2, k) = 0.910

             ! Introduced for aerosol indirect forcings.  The
             ! following values use the cloud optical thickness
             ! calculated from present-day aerosol concentrations
             ! whereas the quantities without the "A" at the end are
             ! for pre-industial (natural-only) aerosol concentrations
             PTAUA(i, 1, k) = MAX(cldtaupd(iof+i, k), 1e-05)
             ! (1e-12 serait instable)
             PTAUA(i, 2, k) = MAX(cldtaupd(iof+i, k), 1e-05)
             ! (pour 32-bit machines)
             POMEGAA(i, 1, k) = 0.9999 - 5e-04 * EXP(-0.5 * PTAUA(i, 1, k))
             POMEGAA(i, 2, k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i, 2, k))
             !jq-end
          ENDDO
       ENDDO

       DO k = 1, klev+1
          DO i = 1, kdlon
             PPMB(i, k) = paprs(iof+i, k)/100.
          ENDDO
       ENDDO

       DO kk = 1, 5
          DO k = 1, klev
             DO i = 1, kdlon
                PAER(i, k, kk) = 1E-15
             ENDDO
          ENDDO
       ENDDO

       CALL LW(PPMB, PDP, PDT0, PEMIS, PTL, PTAVE, PWV, POZON, PAER, PCLDLD, &
            PCLDLU, PVIEW, zcool, zcool0, ztoplw, zsollw, ztoplw0, zsollw0, &
            zsollwdown, ZFLUP, ZFLDN, ZFLUP0, ZFLDN0)
       CALL SW(PSCT, zrmu0, zfract, PPMB, PDP, PPSOL, PALBD, PALBP, PTAVE, &
            PWV, PQS, POZON, PCLDSW, PTAU, POMEGA, PCG, zheat, zheat0, &
            zalbpla, ztopsw, zsolsw, ztopsw0, zsolsw0, ZFSUP, ZFSDN, ZFSUP0, &
            ZFSDN0, ztopswad, zsolswad, ztopswai, zsolswai, ok_ade)

       DO i = 1, kdlon
          radsol(iof+i) = zsolsw(i) + zsollw(i)
          topsw(iof+i) = ztopsw(i)
          toplw(iof+i) = ztoplw(i)
          solsw(iof+i) = zsolsw(i)
          sollw(iof+i) = zsollw(i)
          sollwdown(iof+i) = zsollwdown(i)

          DO k = 1, klev+1
             lwdn0 ( iof+i, k)   = ZFLDN0 ( i, k)
             lwdn  ( iof+i, k)   = ZFLDN  ( i, k)
             lwup0 ( iof+i, k)   = ZFLUP0 ( i, k)
             lwup  ( iof+i, k)   = ZFLUP  ( i, k)
          ENDDO

          topsw0(iof+i) = ztopsw0(i)
          toplw0(iof+i) = ztoplw0(i)
          solsw0(iof+i) = zsolsw0(i)
          sollw0(iof+i) = zsollw0(i)
          albpla(iof+i) = zalbpla(i)

          DO k = 1, klev+1
             swdn0 ( iof+i, k)   = ZFSDN0 ( i, k)
             swdn  ( iof+i, k)   = ZFSDN  ( i, k)
             swup0 ( iof+i, k)   = ZFSUP0 ( i, k)
             swup  ( iof+i, k)   = ZFSUP  ( i, k)
          ENDDO
       ENDDO
       ! transform the aerosol forcings, if they have to be calculated
       IF (ok_ade) THEN
          DO i = 1, kdlon
             topswad(iof+i) = ztopswad(i)
             solswad(iof+i) = zsolswad(i)
          ENDDO
       ELSE
          DO i = 1, kdlon
             topswad(iof+i) = 0.
             solswad(iof+i) = 0.
          ENDDO
       ENDIF

       DO k = 1, klev
          DO i = 1, kdlon
             ! scale factor to take into account the difference
             ! between dry air and water vapour specific heat capacity
             zznormcp = 1. + RVTMP2 * PWV(i, k)
             heat(iof+i, k) = zheat(i, k) / zznormcp
             cool(iof+i, k) = zcool(i, k)/zznormcp
             heat0(iof+i, k) = zheat0(i, k)/zznormcp
             cool0(iof+i, k) = zcool0(i, k)/zznormcp
          ENDDO
       ENDDO
    end DO loop_iof

  END SUBROUTINE radlwsw

end module radlwsw_m
1	guez	53	module radlwsw_m
2	guez	3
3	guez	53	IMPLICIT none
4	guez	3
5	guez	53	contains
6
7	guez	212	SUBROUTINE radlwsw(dist, mu0, fract, paprs, play, tsol, albedo, t, q, wo, &
8			cldfra, cldemi, cldtaupd, heat, heat0, cool, cool0, radsol, albpla, &
9			topsw, toplw, solsw, sollw, sollwdown, topsw0, toplw0, solsw0, sollw0, &
10	guez	217	lwdn0, lwdn, lwup0, lwup, swdn0, swdn, swup0, swup, ok_ade, topswad, &
11			solswad)
12	guez	53
13	guez	217	! From LMDZ4/libf/phylmd/radlwsw.F, version 1.4, 2005/06/06 13:16:33
14	guez	69	! Author: Z. X. Li (LMD/CNRS)
15			! Date: 1996/07/19
16	guez	53
17	guez	69	! Objet : interface entre le modèle et les rayonnements solaire et
18			! infrarouge
19
20			! ATTENTION: swai and swad have to be interpreted in the following manner:
21
22	guez	217	! not ok_ade
23	guez	69	! both are zero
24
25	guez	217	! ok_ade
26	guez	69	! aerosol direct forcing is F_{AD} = topsw - topswad
27			! indirect is zero
28
29	guez	118	USE clesphys, ONLY: solaire
30	guez	53	USE dimphy, ONLY: klev, klon
31	guez	71	use lw_m, only: lw
32			USE raddim, ONLY: kdlon
33	guez	53	USE suphec_m, ONLY: rg
34	guez	71	use sw_m, only: sw
35	guez	53	USE yoethf_m, ONLY: rvtmp2
36	guez	217
37	guez	118	real, intent(in):: dist ! distance astronomique terre-soleil
38			real, intent(in):: mu0(klon) ! cosinus de l'angle zenithal
39			real, intent(in):: fract(klon) ! duree d'ensoleillement normalisee
40			real, intent(in):: paprs(klon, klev+1) ! pression a inter-couche (Pa)
41			real, intent(in):: play(klon, klev) ! pression au milieu de couche (Pa)
42	guez	155	real, intent(in):: tsol(klon) ! temperature du sol (en K)
43			real, intent(in):: albedo(klon) ! albedo du sol (entre 0 et 1)
44			real, intent(in):: t(klon, klev) ! temperature (K)
45	guez	212	real, intent(in):: q(klon, klev) ! vapeur d'eau (en kg/kg)
46	guez	118
47	guez	69	real, intent(in):: wo(klon, klev)
48	guez	118	! column-density of ozone in a layer, in kilo-Dobsons
49
50	guez	212	real, intent(in):: cldfra(klon, klev) ! fraction nuageuse (entre 0 et 1)
51	guez	53
52	guez	212	real, intent(in):: cldemi(klon, klev)
53			! emissivite des nuages dans l'IR (entre 0 et 1)
54	guez	53
55	guez	212	real, intent(in):: cldtaupd(klon, klev)
56			! epaisseur optique des nuages dans le visible (present-day value)
57
58	guez	53	real, intent(out):: heat(klon, klev)
59			! échauffement atmosphérique (visible) (K/jour)
60
61	guez	213	real, intent(out):: heat0(klon, klev) ! chauffage solaire ciel clair
62	guez	212	real, intent(out):: cool(klon, klev) ! refroidissement dans l'IR (K/jour)
63	guez	213
64	guez	212	real, intent(out):: cool0(klon, klev)
65	guez	213	! refroidissement infrarouge ciel clair
66	guez	212
67			real, intent(out):: radsol(klon)
68			! bilan radiatif net au sol (W/m**2) (+ vers le bas)
69
70	guez	155	real, intent(out):: albpla(klon) ! albedo planetaire (entre 0 et 1)
71	guez	212	real, intent(out):: topsw(klon) ! flux solaire net au sommet de l'atm.
72	guez	62
73			real, intent(out):: toplw(klon)
74			! rayonnement infrarouge montant au sommet de l'atmosphère
75
76	guez	72	real, intent(out):: solsw(klon) ! flux solaire net à la surface
77
78			real, intent(out):: sollw(klon)
79			! rayonnement infrarouge montant à la surface
80
81			real, intent(out):: sollwdown(klon)
82	guez	212	real, intent(out):: topsw0(klon)
83	guez	62	real, intent(out):: toplw0(klon)
84	guez	212	real, intent(out):: solsw0(klon), sollw0(klon)
85			REAL, intent(out):: lwdn0(klon, klev+1), lwdn(klon, klev+1)
86			REAL, intent(out):: lwup0(klon, klev+1), lwup(klon, klev+1)
87			REAL, intent(out):: swdn0(klon, klev+1), swdn(klon, klev+1)
88			REAL, intent(out):: swup0(klon, klev+1), swup(klon, klev+1)
89	guez	72
90	guez	212	logical, intent(in):: ok_ade ! apply the Aerosol Direct Effect
91	guez	72
92	guez	212	real, intent(out):: topswad(klon), solswad(klon)
93			! aerosol direct forcing at TOA and surface
94	guez	217	! rayonnement solaire net absorb\'e
95	guez	72
96			! Local:
97
98	guez	53	DOUBLE PRECISION ZFSUP(KDLON, KLEV+1)
99			DOUBLE PRECISION ZFSDN(KDLON, KLEV+1)
100			DOUBLE PRECISION ZFSUP0(KDLON, KLEV+1)
101			DOUBLE PRECISION ZFSDN0(KDLON, KLEV+1)
102
103			DOUBLE PRECISION ZFLUP(KDLON, KLEV+1)
104			DOUBLE PRECISION ZFLDN(KDLON, KLEV+1)
105			DOUBLE PRECISION ZFLUP0(KDLON, KLEV+1)
106			DOUBLE PRECISION ZFLDN0(KDLON, KLEV+1)
107
108			DOUBLE PRECISION zx_alpha1, zx_alpha2
109	guez	62	INTEGER k, kk, i, iof, nb_gr
110	guez	53	DOUBLE PRECISION PSCT
111
112			DOUBLE PRECISION PALBD(kdlon, 2), PALBP(kdlon, 2)
113			DOUBLE PRECISION PEMIS(kdlon), PDT0(kdlon), PVIEW(kdlon)
114			DOUBLE PRECISION PPSOL(kdlon), PDP(kdlon, klev)
115			DOUBLE PRECISION PTL(kdlon, klev+1), PPMB(kdlon, klev+1)
116			DOUBLE PRECISION PTAVE(kdlon, klev)
117	guez	118	DOUBLE PRECISION PWV(kdlon, klev), PQS(kdlon, klev)
118			DOUBLE PRECISION POZON(kdlon, klev) ! mass fraction of ozone
119	guez	178	DOUBLE PRECISION PAER(kdlon, klev, 5) ! AEROSOLS' OPTICAL THICKNESS
120	guez	53	DOUBLE PRECISION PCLDLD(kdlon, klev)
121			DOUBLE PRECISION PCLDLU(kdlon, klev)
122			DOUBLE PRECISION PCLDSW(kdlon, klev)
123			DOUBLE PRECISION PTAU(kdlon, 2, klev)
124			DOUBLE PRECISION POMEGA(kdlon, 2, klev)
125			DOUBLE PRECISION PCG(kdlon, 2, klev)
126
127	guez	118	DOUBLE PRECISION zfract(kdlon), zrmu0(kdlon)
128	guez	53
129			DOUBLE PRECISION zheat(kdlon, klev), zcool(kdlon, klev)
130			DOUBLE PRECISION zheat0(kdlon, klev), zcool0(kdlon, klev)
131			DOUBLE PRECISION ztopsw(kdlon), ztoplw(kdlon)
132			DOUBLE PRECISION zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon)
133			DOUBLE PRECISION zsollwdown(kdlon)
134
135			DOUBLE PRECISION ztopsw0(kdlon), ztoplw0(kdlon)
136			DOUBLE PRECISION zsolsw0(kdlon), zsollw0(kdlon)
137			DOUBLE PRECISION zznormcp
138
139			!jq the following quantities are needed for the aerosol radiative forcings
140
141			DOUBLE PRECISION PTAUA(kdlon, 2, klev)
142			! present-day value of cloud opt thickness (PTAU is pre-industrial
143			! value), local use
144
145			DOUBLE PRECISION POMEGAA(kdlon, 2, klev) ! dito for single scatt albedo
146
147			DOUBLE PRECISION ztopswad(kdlon), zsolswad(kdlon)
148			! Aerosol direct forcing at TOAand surface
149
150			DOUBLE PRECISION ztopswai(kdlon), zsolswai(kdlon) ! dito, indirect
151	guez	118	real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
152	guez	53
153			!----------------------------------------------------------------------
154
155			nb_gr = klon / kdlon
156			IF (nb_gr * kdlon /= klon) THEN
157			PRINT *, "kdlon mauvais :", klon, kdlon, nb_gr
158			stop 1
159			ENDIF
160	guez	217
161	guez	53	heat = 0.
162			cool = 0.
163			heat0 = 0.
164			cool0 = 0.
165	guez	118	PSCT = solaire / dist**2
166	guez	53
167	guez	62	loop_iof: DO iof = 0, klon - kdlon, kdlon
168	guez	53	DO i = 1, kdlon
169			zfract(i) = fract(iof+i)
170	guez	118	zrmu0(i) = mu0(iof+i)
171	guez	53	PALBD(i, 1) = albedo(iof+i)
172	guez	155	PALBD(i, 2) = albedo(iof+i)
173	guez	53	PALBP(i, 1) = albedo(iof+i)
174	guez	155	PALBP(i, 2) = albedo(iof+i)
175	guez	53	! cf. JLD pour etre en accord avec ORCHIDEE il faut mettre
176			! PEMIS(i) = 0.96
177	guez	217	PEMIS(i) = 1.
178	guez	53	PVIEW(i) = 1.66
179			PPSOL(i) = paprs(iof+i, 1)
180	guez	72	zx_alpha1 = (paprs(iof+i, 1)-play(iof+i, 2)) &
181			/ (play(iof+i, 1)-play(iof+i, 2))
182	guez	217	zx_alpha2 = 1. - zx_alpha1
183	guez	53	PTL(i, 1) = t(iof+i, 1) * zx_alpha1 + t(iof+i, 2) * zx_alpha2
184			PTL(i, klev+1) = t(iof+i, klev)
185			PDT0(i) = tsol(iof+i) - PTL(i, 1)
186			ENDDO
187			DO k = 2, klev
188			DO i = 1, kdlon
189			PTL(i, k) = (t(iof+i, k)+t(iof+i, k-1))*0.5
190			ENDDO
191			ENDDO
192			DO k = 1, klev
193			DO i = 1, kdlon
194			PDP(i, k) = paprs(iof+i, k)-paprs(iof+i, k+1)
195			PTAVE(i, k) = t(iof+i, k)
196	guez	217	PWV(i, k) = MAX (q(iof+i, k), 1e-12)
197	guez	53	PQS(i, k) = PWV(i, k)
198	guez	118	POZON(i, k) = wo(iof+i, k) * RG * dobson_u * 1e3 &
199			/ (paprs(iof+i, k) - paprs(iof+i, k+1))
200	guez	53	PCLDLD(i, k) = cldfra(iof+i, k)*cldemi(iof+i, k)
201			PCLDLU(i, k) = cldfra(iof+i, k)*cldemi(iof+i, k)
202			PCLDSW(i, k) = cldfra(iof+i, k)
203	guez	217	PTAU(i, 1, k) = MAX(cldtaupd(iof+i, k), 1e-05)
204	guez	53	! (1e-12 serait instable)
205	guez	217	PTAU(i, 2, k) = MAX(cldtaupd(iof+i, k), 1e-05)
206	guez	53	! (pour 32-bit machines)
207	guez	217	POMEGA(i, 1, k) = 0.9999 - 5e-04 * EXP(-0.5 * PTAU(i, 1, k))
208	guez	53	POMEGA(i, 2, k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i, 2, k))
209			PCG(i, 1, k) = 0.865
210			PCG(i, 2, k) = 0.910
211
212			! Introduced for aerosol indirect forcings. The
213			! following values use the cloud optical thickness
214			! calculated from present-day aerosol concentrations
215			! whereas the quantities without the "A" at the end are
216			! for pre-industial (natural-only) aerosol concentrations
217	guez	217	PTAUA(i, 1, k) = MAX(cldtaupd(iof+i, k), 1e-05)
218	guez	53	! (1e-12 serait instable)
219	guez	217	PTAUA(i, 2, k) = MAX(cldtaupd(iof+i, k), 1e-05)
220	guez	53	! (pour 32-bit machines)
221	guez	217	POMEGAA(i, 1, k) = 0.9999 - 5e-04 * EXP(-0.5 * PTAUA(i, 1, k))
222	guez	53	POMEGAA(i, 2, k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i, 2, k))
223			!jq-end
224			ENDDO
225			ENDDO
226
227			DO k = 1, klev+1
228			DO i = 1, kdlon
229	guez	217	PPMB(i, k) = paprs(iof+i, k)/100.
230	guez	53	ENDDO
231			ENDDO
232
233			DO kk = 1, 5
234			DO k = 1, klev
235			DO i = 1, kdlon
236	guez	217	PAER(i, k, kk) = 1E-15
237	guez	53	ENDDO
238			ENDDO
239			ENDDO
240
241	guez	145	CALL LW(PPMB, PDP, PDT0, PEMIS, PTL, PTAVE, PWV, POZON, PAER, PCLDLD, &
242			PCLDLU, PVIEW, zcool, zcool0, ztoplw, zsollw, ztoplw0, zsollw0, &
243			zsollwdown, ZFLUP, ZFLDN, ZFLUP0, ZFLDN0)
244	guez	53	CALL SW(PSCT, zrmu0, zfract, PPMB, PDP, PPSOL, PALBD, PALBP, PTAVE, &
245	guez	178	PWV, PQS, POZON, PCLDSW, PTAU, POMEGA, PCG, zheat, zheat0, &
246	guez	53	zalbpla, ztopsw, zsolsw, ztopsw0, zsolsw0, ZFSUP, ZFSDN, ZFSUP0, &
247	guez	217	ZFSDN0, ztopswad, zsolswad, ztopswai, zsolswai, ok_ade)
248	guez	53
249			DO i = 1, kdlon
250			radsol(iof+i) = zsolsw(i) + zsollw(i)
251			topsw(iof+i) = ztopsw(i)
252			toplw(iof+i) = ztoplw(i)
253			solsw(iof+i) = zsolsw(i)
254			sollw(iof+i) = zsollw(i)
255			sollwdown(iof+i) = zsollwdown(i)
256
257			DO k = 1, klev+1
258			lwdn0 ( iof+i, k) = ZFLDN0 ( i, k)
259			lwdn ( iof+i, k) = ZFLDN ( i, k)
260			lwup0 ( iof+i, k) = ZFLUP0 ( i, k)
261			lwup ( iof+i, k) = ZFLUP ( i, k)
262			ENDDO
263
264			topsw0(iof+i) = ztopsw0(i)
265			toplw0(iof+i) = ztoplw0(i)
266			solsw0(iof+i) = zsolsw0(i)
267			sollw0(iof+i) = zsollw0(i)
268			albpla(iof+i) = zalbpla(i)
269
270			DO k = 1, klev+1
271			swdn0 ( iof+i, k) = ZFSDN0 ( i, k)
272			swdn ( iof+i, k) = ZFSDN ( i, k)
273			swup0 ( iof+i, k) = ZFSUP0 ( i, k)
274			swup ( iof+i, k) = ZFSUP ( i, k)
275			ENDDO
276			ENDDO
277			! transform the aerosol forcings, if they have to be calculated
278			IF (ok_ade) THEN
279			DO i = 1, kdlon
280			topswad(iof+i) = ztopswad(i)
281			solswad(iof+i) = zsolswad(i)
282			ENDDO
283			ELSE
284			DO i = 1, kdlon
285	guez	217	topswad(iof+i) = 0.
286			solswad(iof+i) = 0.
287	guez	53	ENDDO
288			ENDIF
289
290			DO k = 1, klev
291			DO i = 1, kdlon
292			! scale factor to take into account the difference
293			! between dry air and water vapour specific heat capacity
294			zznormcp = 1. + RVTMP2 * PWV(i, k)
295			heat(iof+i, k) = zheat(i, k) / zznormcp
296			cool(iof+i, k) = zcool(i, k)/zznormcp
297			heat0(iof+i, k) = zheat0(i, k)/zznormcp
298			cool0(iof+i, k) = zcool0(i, k)/zznormcp
299			ENDDO
300			ENDDO
301	guez	62	end DO loop_iof
302	guez	53
303			END SUBROUTINE radlwsw
304
305			end module radlwsw_m