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