source: trunk/LMDZ.MARS/libf/phymars/swmain_mod.F @ 2311

Last change on this file since 2311 was 2311, checked in by emillour, 5 years ago

Mars GCM:
Code tidying: use getin_p() instead of getin() and use "call abort_physic"
instead of "stop" or "call abort".
EM

File size: 6.7 KB
Line 
1      MODULE swmain_mod
2     
3      IMPLICIT NONE
4
5      CONTAINS
6
7      SUBROUTINE SWMAIN ( KDLON, KFLEV,
8     $                PCST, albedo,
9     $                PRMU0, PDP, PPLEV, aerosol,PFRACT,
10     $                PHEAT, PFLUXD,PFLUXU,
11     &                QVISsQREF3d,omegaVIS3d,gVIS3d)
12
13      use dimradmars_mod, only: ndlo2, ndlon, nflev,
14     &                          nsun,naerkind
15      use yomlw_h, only: nlaylte, gcp
16      IMPLICIT NONE
17
18c     DECLARATIONS
19c     -------------     
20      include "callkeys.h"
21c     
22c     PURPOSE.
23c     --------
24c     
25c     This routine computes the shortwave (solar wavelength)
26c     radiation fluxes in two spectral intervals
27c     and heating rate on the first "nlaylte" layers.
28C     
29c     Francois Forget (2000), adapted from
30C     Fouquart and Bonnel's ECMWF program
31c
32C     IMPLICIT ARGUMENTS :
33C     --------------------
34C     
35C     ==== INPUTS ===
36c
37c    KDLON :  number of horizontal grid points
38c    PCST  :  Solar constant on Mars (W.m-2)
39c    albedo   hemispheric surface albedo
40c                         albedo (i,1) : mean albedo for solar band#1
41c                                        (see below)
42c                         albedo (i,2) : mean albedo for solar band#2
43c                                        (see below)
44c    PRMU0 :         cos of solar zenith angle (=1 when sun at zenith)
45c    PDP :           Layer thickness (Pa)
46c    PPLEV           pressure (Pa) at boundaries of each layer
47c   aerosol               aerosol extinction optical depth
48c                         at reference wavelength "longrefvis" set
49c                         in dimradmars_mod , in each layer, for one of
50c                         the "naerkind" kind of aerosol optical properties.
51c    Pfract :        day fraction of the time interval
52c                          =1 during the full day ; =0 during the night
53c    QVISsQREF3d,omegaVIS3d,gVIS3d Aerosol optical properties
54c
55C     ==== OUTPUTS ===
56c    PHEAT :         Heating rate (K/s)
57c    PFLUXD :        SW downward flux at boundaries of each layer (W.m-2)
58c    PFLUXU :        SW upward flux at boundaries of each layer (W.m-2)
59C     
60C     ----------
61C     
62C-----------------------------------------------------------------------
63C     
64C     
65C-----------------------------------------------------------------------
66C     
67C     ARGUMENTS
68C     ---------
69c     INPUTS/OUTPUTS:
70c     ---------
71     
72      INTEGER, iNTENT(IN) :: KDLON, KFLEV
73      REAL, iNTENT(IN) :: aerosol(NDLO2,KFLEV,naerkind),PRMU0(NDLO2)
74      REAL, iNTENT(IN) :: PCST
75      REAL, iNTENT(IN) :: albedo(NDLO2,2)
76      REAL, iNTENT(IN) :: PDP(NDLO2,KFLEV)
77      REAL, iNTENT(IN) :: PPLEV(NDLO2,KFLEV+1)
78      REAL, iNTENT(OUT) :: PHEAT(NDLO2,KFLEV)
79      REAL, iNTENT(IN) :: PFRACT(NDLO2)
80      REAL, iNTENT(OUT) :: PFLUXD(NDLON,NFLEV+1,2)
81      REAL, iNTENT(OUT) :: PFLUXU(NDLON,NFLEV+1,2)
82      REAL, iNTENT(IN) :: QVISsQREF3d(NDLO2,KFLEV,nsun,naerkind)
83      REAL, iNTENT(IN) :: omegaVIS3d(NDLO2,KFLEV,nsun,naerkind)
84      REAL, iNTENT(IN) :: gVIS3d(NDLO2,KFLEV,nsun,naerkind)
85     
86C     LOCAL ARRAYS
87C     ------------
88      REAL ZPSOL(NDLO2)
89      REAL ZDSIG(NDLON,NFLEV), ZFACT(NDLON)
90     S     ,  ZFD(NDLON,NFLEV+1)
91     S     ,  ZFU(NDLON,NFLEV+1)
92     S     ,  ZRMU(NDLON), ZSEC(NDLON)
93     S     ,  ZUD(NDLON,3,NFLEV+1), ZUM(NDLON,NFLEV+1)
94      REAL ZSIGN(NDLON),  ZSIGO(NDLON)
95
96c following line has been changed, kflev--->nflev (to avoid error message
97c when compiling on NASA Ames Sun)
98      REAL ZFDOWN(NDLO2,NFLEV+1),ZFUP(NDLO2,NFLEV+1)
99
100      integer jl, jk, jkp1, jkl
101      integer INU
102      real  zdfnet
103
104C     ------------------------------------------------------------------
105C     Initializations :
106C     ------------------------------------------------------------------
107
108c     Incident Solar flux and corrected angle in the atmosphere
109c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
110      DO JL = 1 , KDLON
111c        Incident Flux at the top of the atmosphere
112         ZFACT(JL)= PRMU0(JL) * PCST * PFRACT(JL)
113
114c        Cos of solar zenith angle CORRECTED for high zenith angle
115       if (PRMU0(JL).GT.0) then
116         ZRMU(JL)=SQRT(1224.* PRMU0(JL) * PRMU0(JL) + 1.) / 35.
117       else
118         ZRMU(JL)= 1. / 35.
119       endif
120         ZSEC(JL)=1./ZRMU(JL)
121      END DO
122
123c     Calcul of ZDSIG  (thickness of layers in sigma coordinates)
124c     ~~~~~~~~~~~~~~~
125      DO JL = 1 , KDLON
126         ZSIGO(JL) = 1.0
127         ZPSOL(JL) = PPLEV(JL,1)
128      END DO
129      DO JK = 1 , nlaylte
130         JKP1 = JK + 1
131         JKL = nlaylte+1 - JK
132         DO  JL = 1 , KDLON
133            ZSIGN(JL) = PPLEV(JL,JKP1) / PPLEV(JL,1)
134            ZDSIG(JL,JK) = ZSIGO(JL) - ZSIGN(JL)
135            ZSIGO(JL) = ZSIGN(JL)
136         END DO
137      END DO
138
139C------------------------------------------------------------------
140C          LOOP ON SPECTRAL INTERVAL in solar spectrum
141C------------------------------------------------------------------
142c  2 spectral interval in solar spectrum :
143c  - INU=1: between wavelength "long1vis" and "long2vis" set in dimradmars_mod
144c  - INU=2: between wavelength "long2vis" and "long3vis" set in dimradmars_mod
145
146      DO INU = 1,2
147
148! NB: swrtype is set in callkeys.h
149        if (swrtype.eq.1) then ! Fouquart
150          CALL SWR_FOUQUART( KDLON, kflev, INU
151     &     ,  aerosol,QVISsQREF3d,omegaVIS3d,gVIS3d
152     &     ,  albedo,ZDSIG,ZPSOL,ZRMU,ZSEC
153     &     ,  ZFD,ZFU )
154        else
155          if (swrtype.eq.2) then ! Toon
156            CALL SWR_TOON( KDLON, kflev, INU
157     &       ,  aerosol,QVISsQREF3d,omegaVIS3d,gVIS3d
158     &       ,  albedo,ZDSIG,ZPSOL,ZRMU,ZSEC
159     &       ,  ZFD,ZFU )
160          else
161            write(*,*) "swmain: invalid swrtype value !!"
162            call abort_physic("swmain","invalid swrtype",1)
163          endif ! of if (swrtype.eq.2)
164        endif ! of if (swrtype.eq.1)
165       
166         DO JK = 1 , nlaylte+1
167            DO  JL = 1 , KDLON
168              PFLUXD(JL,JK,INU)=ZFD(JL,JK)*ZFACT(JL)
169              PFLUXU(JL,JK,INU)=ZFU(JL,JK)*ZFACT(JL)
170            END DO
171         END DO
172      END DO ! of DO INU=1,2
173
174C     ------------------------------------------------------
175C     HEATING RATES
176C     ------------------------------------------------------
177
178      DO JK = 1 , nlaylte+1
179         DO  JL = 1 , KDLON
180c           wavelength integrated flux at every level:
181            ZFUP(JL,JK)= (PFLUXU(JL,JK,1)+ PFLUXU(JL,JK,2))
182            ZFDOWN(JL,JK)= (PFLUXD(JL,JK,1)+ PFLUXD(JL,JK,2))
183         END DO
184      END DO
185     
186      DO JK = 1 , nlaylte
187         DO  JL = 1 , KDLON
188            ZDFNET = ZFUP (JL,JK  ) - ZFDOWN(JL,JK  )
189     S              -ZFUP (JL,JK+1) + ZFDOWN(JL,JK+1)
190c           Heating rate
191            PHEAT(JL,JK) = gcp * ZDFNET / PDP(JL,JK)
192         END DO
193      END DO
194
195      END SUBROUTINE SWMAIN
196
197      END MODULE swmain_mod
Note: See TracBrowser for help on using the repository browser.