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, ok_aie, &
       tau_ae, piz_ae, cg_ae, topswad, solswad, cldtaupi, topswai, solswai)

    ! 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 and not ok_aie
    ! both are zero

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

    ! not ok_ade and ok_aie
    ! aerosol indirect forcing is F_{AI} = topsw - topswai
    ! direct is zero

    ! ok_ade and ok_aie
    ! aerosol indirect forcing is F_{AI} = topsw - topswai
    ! aerosol direct forcing is F_{AD} = topswai - topswad

    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)
    real, intent(out):: cool(klon, klev) ! refroidissement dans l'IR (K/jour)
    real, intent(out):: cool0(klon, klev)

    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
    logical, intent(in):: ok_aie ! apply the Aerosol Indirect Effect

    ! aerosol optical properties (calculated in aeropt.F):
    real, intent(in):: tau_ae(klon, klev, 2), piz_ae(klon, klev, 2)
    real, intent(in):: cg_ae(klon, klev, 2)

    real, intent(out):: topswad(klon), solswad(klon)
    ! aerosol direct forcing at TOA and surface
    ! ray. solaire net absorbe
    
    real, intent(in):: cldtaupi(klon, klev)
    ! cloud visible optical thickness for pre-industrial aerosol concentrations
    ! i.e. with smaller droplet concentration, thus larger droplets,
    ! thus generally cdltaupi cldtaupd it is needed for the
    ! diagnostics of the aerosol indirect radiative forcing

    real, intent(out):: topswai(klon), solswai(klon)
    ! aerosol indirect forcing at TOA and surface
    ! ray. solaire net absorbe

    ! Local:

    double precision tauae(kdlon, klev, 2) ! aer opt properties
    double precision pizae(kdlon, klev, 2)
    double precision cgae(kdlon, klev, 2)

    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

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

    tauae = 0.
    pizae = 0.
    cgae = 0.

    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.0 
          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.0 - 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), 1.0e-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(cldtaupi(iof+i, k), 1.0e-05)
             ! (1e-12 serait instable)
             PTAU(i, 2, k) = MAX(cldtaupi(iof+i, k), 1.0e-05)
             ! (pour 32-bit machines)
             POMEGA(i, 1, k) = 0.9999 - 5.0e-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), 1.0e-05)
             ! (1e-12 serait instable)
             PTAUA(i, 2, k) = MAX(cldtaupd(iof+i, k), 1.0e-05)
             ! (pour 32-bit machines)
             POMEGAA(i, 1, k) = 0.9999 - 5.0e-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.0
          ENDDO
       ENDDO

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

       DO k = 1, klev
          DO i = 1, kdlon
             tauae(i, k, 1) = tau_ae(iof+i, k, 1)
             pizae(i, k, 1) = piz_ae(iof+i, k, 1)
             cgae(i, k, 1) =cg_ae(iof+i, k, 1)
             tauae(i, k, 2) = tau_ae(iof+i, k, 2)
             pizae(i, k, 2) = piz_ae(iof+i, k, 2)
             cgae(i, k, 2) =cg_ae(iof+i, k, 2)
          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, tauae, pizae, cgae, PTAUA, POMEGAA, ztopswad, zsolswad, &
            ztopswai, zsolswai, ok_ade, ok_aie)

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