Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 1982)
+++ /trunk/LMDZ.MARS/README	(revision 1983)
@@ -2626,2 +2626,8 @@
 Update xvik.F main program to work with current diagfi outputs.
 
+== 01/08/2018 == MV
+Cosmetic/practical changes:
+- swmain and lwmain become modules swmain_mod, lwmain_mod
+- Addition of the intent in/out characteristics of variables in swmain_mod and lwmain_mod subroutines
+Correction:
+- in callsedim_mod, declaration of variable tau(ngrid,nlay) corrected to tau(ngrid,naerkind)
Index: /trunk/LMDZ.MARS/libf/phymars/callradite_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/callradite_mod.F	(revision 1982)
+++ /trunk/LMDZ.MARS/libf/phymars/callradite_mod.F	(revision 1983)
@@ -23,4 +23,6 @@
       use comcstfi_h, only: g,cpp
       use time_phylmdz_mod, only: daysec
+      use lwmain_mod, only: lwmain
+      use swmain_mod, only: swmain
       IMPLICIT NONE
 c=======================================================================
Index: /trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F	(revision 1982)
+++ /trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F	(revision 1983)
@@ -22,4 +22,5 @@
       USE newsedim_mod, ONLY: newsedim
       USE comcstfi_h, ONLY: g
+      USE dimradmars_mod, only: naerkind
       IMPLICIT NONE
 
@@ -74,5 +75,5 @@
       real,intent(out) :: pdqs_sed(ngrid,nq)    ! flux at surface (kg.m-2.s-1)
       integer,intent(in) :: nq  ! number of tracers
-      real,intent(in) :: tau(ngrid,nlay) ! dust opacity
+      real,intent(in) :: tau(ngrid,naerkind) ! dust opacity
       real,intent(in) :: tauscaling(ngrid)
       
Index: unk/LMDZ.MARS/libf/phymars/lwmain.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/lwmain.F	(revision 1982)
+++ 	(revision )
@@ -1,196 +1,0 @@
-       subroutine lwmain (ig0,icount,kdlon,kflev
-     .                   ,dp,dt0,emis
-     .                   ,plev,tlev,tlay,aerosol,coolrate
-     .                   ,fluxground,fluxtop
-     .                   ,netrad
-     &                   ,QIRsQREF3d,omegaIR3d,gIR3d
-     &                   ,co2ice)
-
-c----------------------------------------------------------------------
-c     LWMAIN     organizes the LTE longwave calculations
-c     for layer 1 to layer "nlaylte" (stored in  "yomlw_h")
-c----------------------------------------------------------------------
-
-      use dimradmars_mod, only: ndlo2, nflev, nir, ndlon, nuco2
-      use dimradmars_mod, only: naerkind
-      use yomlw_h, only: nlaylte, xi
-      implicit none
- 
-#include "callkeys.h"
-#include "comg1d.h"
-
-c----------------------------------------------------------------------
-c         0.1   arguments
-c               ---------
-c                                                            inputs:
-c                                                            -------
-      integer ig0
-      integer icount
-      integer kdlon            ! part of ngrid
-      integer kflev            ! part of nlayer
-
-      real dp (ndlo2,kflev)         ! layer pressure thickness (Pa)
-      real dt0 (ndlo2)              ! surface temperature discontinuity (K)
-      real emis (ndlo2)             ! surface emissivity
-      real plev (ndlo2,kflev+1)     ! level pressure (Pa)
-      real tlev (ndlo2,kflev+1)     ! level temperature (K)
-      real tlay (ndlo2,kflev)       ! layer temperature (K)
-      real aerosol(ndlo2,kflev,naerkind)      !  aerosol extinction optical
-c                         depth at reference wavelength "longrefvis" set
-c                         in dimradmars_mod , in each layer, for one of
-c                         the "naerkind" kind of aerosol optical properties.
-
-
-c                                                            outputs:
-c                                                            --------
-      real coolrate(ndlo2,kflev)      ! cooling rate (K/s)
-      real fluxground(ndlo2)          ! downward ground flux (W/m2)
-      real fluxtop(ndlo2)             ! outgoing upward flux (W/m2) ("OLR")
-      real netrad (ndlo2,kflev)       ! radiative budget (W/m2)
-c     Aerosol optical properties
-      REAL :: QIRsQREF3d(ndlo2,kflev,nir,naerkind)
-      REAL :: omegaIR3d(ndlo2,kflev,nir,naerkind)
-      REAL :: gIR3d(ndlo2,kflev,nir,naerkind)
-
-c----------------------------------------------------------------------
-c         0.2   local arrays
-c               ------------
-
-      real aer_t (ndlon,nuco2,nflev+1)  ! transmission (aer)
-      real co2_u (ndlon,nuco2,nflev+1)  ! absorber amounts (co2)
-      real co2_up (ndlon,nuco2,nflev+1) ! idem scaled by the pressure (co2)
-
-      real bsurf (ndlon,nir)            ! surface spectral planck function
-      real btop (ndlon,nir)             ! top spectral planck function
-      real blev (ndlon,nir,nflev+1)     ! level   spectral planck function
-      real blay (ndlon,nir,nflev)       ! layer   spectral planck function
-      real dblay (ndlon,nir,nflev)      ! layer gradient spectral planck function
-      real dbsublay (ndlon,nir,2*nflev) ! layer gradient spectral planck function
-                                        ! in sub layers
-
-      real tautotal(ndlon,nflev,nir)  ! \   Total single scattering
-      real omegtotal(ndlon,nflev,nir) !  >  properties (Addition of the
-      real gtotal(ndlon,nflev,nir)    ! /   NAERKIND aerosols prop.)
-
-      real newcoolrate(ndlon,nflev) ! cooling rate (K/s) / with implicite scheme
-
-      REAL co2ice(ndlo2)           ! co2 ice surface layer (kg.m-2)
-      REAL emis_gaz(ndlo2)         ! emissivity for gaz computations
-
-      integer jk,jkk,ja,jl
-
-
-c----------------------------------------------------------------------
-c         0.3   Initialisation
-c               --------------
-
-      DO jl=1 , kdlon
-         IF(co2ice(jl) .GT. 20.e-3) THEN
-             emis_gaz(jl)=1. 
-         ELSE
-             emis_gaz(jl)=emis(jl)
-         ENDIF
-      ENDDO
-
-c----------------------------------------------------------------------
-c         1.0   planck function
-c               ---------------
-
-      call lwb ( kdlon, kflev, tlev, tlay, dt0
-     .         , bsurf, btop, blay, blev, dblay, dbsublay)
-
-c----------------------------------------------------------------------
-c         2.0   absorber amounts
-c               ----------------
-
-      call lwu ( kdlon, kflev
-     .         , dp, plev, tlay, aerosol
-     &         , QIRsQREF3d,omegaIR3d,gIR3d
-     .         , aer_t, co2_u, co2_up
-     .         , tautotal,omegtotal,gtotal)
-
-c----------------------------------------------------------------------
-c         3.0   transmission functions / exchange coefficiants
-c               ----------------------------------------------
-
-c                                                                distants
-c                                                                --------
-                    if( mod(icount-1,ilwd).eq.0) then
-
-c     print*, 'CALL of DISTANTS'
-      call lwxd ( ig0, kdlon, kflev, emis_gaz
-     .          , aer_t, co2_u, co2_up)
-
-                    endif
-c                                                              neighbours
-c                                                              ----------
-                    if( mod(icount-1,ilwn).eq.0) then
-
-c     print*, 'CALL of NEIGHBOURS'
-      call lwxn ( ig0, kdlon, kflev
-     .          , dp
-     .          , aer_t, co2_u, co2_up)
-
-                    endif
-c                                                              boundaries
-c                                                              ----------
-                    if( mod(icount-1,ilwb).eq.0) then
-
-c     print*, 'CALL of BOUNDARIES'
-      call lwxb ( ig0, kdlon, kflev, emis_gaz
-     .          , aer_t, co2_u, co2_up)
-
-                    endif
-
-c----------------------------------------------------------------------
-c         4.0   cooling rate
-c               ------------
-
-      call lwflux ( ig0, kdlon, kflev, dp
-     .            , bsurf, btop, blev, blay, dbsublay
-     .            , tlay, tlev, dt0      ! pour sortie dans g2d uniquement
-     .            , emis
-     .            , tautotal,omegtotal,gtotal
-     .            , coolrate, fluxground, fluxtop
-     .            , netrad)
-
-c     do jk = 1, nlaylte
-c       print*,coolrate(1,jk)
-c     enddo
-      
-c       do jkk = 0 , nlaylte+1
-c         do jk = 0 , nlaylte+1
-c           do ja = 1 , nuco2
-c             do jl = 1 , ngrid
-c      if (xi (jl,ja,jk,jkk) .LT. 0
-c    .       .OR. xi (jl,ja,jk,jkk) .GT. 1 ) then
-c                 print*,'xi bande',ja,jk,jkk,xi (jl,ja,jk,jkk)
-c      endif
-c             enddo
-c           enddo
-c         enddo
-c       enddo
-
-c----------------------------------------------------------------------
-c
-c          5.    shema semi-implicite  (lwi) 
-c                ---------------------------
-c
-c
-      call lwi (ig0,kdlon,kflev,netrad,dblay,dp
-     .          , newcoolrate)
-c
-c  Verif que   (X sol,space) + somme(X i,sol) = 1
-c
-      do jkk = 1 , nlaylte
-        do jl = 1 , kdlon
-c     print*,'NEW et OLD coolrate :',jkk,newcoolrate(jl,jkk) 
-c    .  ,coolrate(jl,jkk)
-      coolrate(jl,jkk) = newcoolrate(jl,jkk)
-        enddo
-      enddo
-c
-c----------------------------------------------------------------------
-
-      return
-      end
Index: /trunk/LMDZ.MARS/libf/phymars/lwmain_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/lwmain_mod.F	(revision 1983)
+++ /trunk/LMDZ.MARS/libf/phymars/lwmain_mod.F	(revision 1983)
@@ -0,0 +1,205 @@
+       MODULE lwmain_mod
+
+       IMPLICIT NONE
+
+       CONTAINS
+       
+       subroutine lwmain (ig0,icount,kdlon,kflev
+     .                   ,dp,dt0,emis
+     .                   ,plev,tlev,tlay,aerosol,coolrate
+     .                   ,fluxground,fluxtop
+     .                   ,netrad
+     &                   ,QIRsQREF3d,omegaIR3d,gIR3d
+     &                   ,co2ice)
+
+c----------------------------------------------------------------------
+c     LWMAIN     organizes the LTE longwave calculations
+c     for layer 1 to layer "nlaylte" (stored in  "yomlw_h")
+c----------------------------------------------------------------------
+
+      use dimradmars_mod, only: ndlo2, nflev, nir, ndlon, nuco2
+      use dimradmars_mod, only: naerkind
+      use yomlw_h, only: nlaylte, xi
+      implicit none
+      
+c     declarations
+c     -------------   
+      include "callkeys.h"
+      include "comg1d.h"
+
+c----------------------------------------------------------------------
+c         0.1   arguments
+c               ---------
+c     inputs/outputs:
+c     -------
+      integer, intent(in) :: ig0
+      integer, intent(in) :: icount
+      integer, intent(in) :: kdlon            ! part of ngrid
+      integer, intent(in) :: kflev            ! part of nlayer
+
+      real, intent(in) :: dp (ndlo2,kflev)         ! layer pressure thickness (Pa)
+      real, intent(in) :: dt0 (ndlo2)              ! surface temperature discontinuity (K)
+      real, intent(in) :: emis (ndlo2)             ! surface emissivity
+      real, intent(in) :: plev (ndlo2,kflev+1)     ! level pressure (Pa)
+      real, intent(in) :: tlev (ndlo2,kflev+1)     ! level temperature (K)
+      real, intent(in) :: tlay (ndlo2,kflev)       ! layer temperature (K)
+      real, intent(in) :: aerosol(ndlo2,kflev,naerkind)      !  aerosol extinction optical
+c                         depth at reference wavelength "longrefvis" set
+c                         in dimradmars_mod , in each layer, for one of
+c                         the "naerkind" kind of aerosol optical properties.
+
+
+c                                                            outputs:
+c                                                            --------
+      real, intent(out) :: coolrate(ndlo2,kflev)      ! cooling rate (K/s)
+      real, intent(out) :: fluxground(ndlo2)          ! downward ground flux (W/m2)
+      real, intent(out) :: fluxtop(ndlo2)             ! outgoing upward flux (W/m2) ("OLR")
+      real, intent(out) :: netrad (ndlo2,kflev)       ! radiative budget (W/m2)
+c     Aerosol optical properties
+      real, intent(in) :: QIRsQREF3d(ndlo2,kflev,nir,naerkind)
+      real, intent(in) :: omegaIR3d(ndlo2,kflev,nir,naerkind)
+      real, intent(in) :: gIR3d(ndlo2,kflev,nir,naerkind)
+      real, intent(in) :: co2ice(ndlo2)           ! co2 ice surface layer (kg.m-2)
+c----------------------------------------------------------------------
+c         0.2   local arrays
+c               ------------
+
+      real aer_t (ndlon,nuco2,nflev+1)  ! transmission (aer)
+      real co2_u (ndlon,nuco2,nflev+1)  ! absorber amounts (co2)
+      real co2_up (ndlon,nuco2,nflev+1) ! idem scaled by the pressure (co2)
+
+      real bsurf (ndlon,nir)            ! surface spectral planck function
+      real btop (ndlon,nir)             ! top spectral planck function
+      real blev (ndlon,nir,nflev+1)     ! level   spectral planck function
+      real blay (ndlon,nir,nflev)       ! layer   spectral planck function
+      real dblay (ndlon,nir,nflev)      ! layer gradient spectral planck function
+      real dbsublay (ndlon,nir,2*nflev) ! layer gradient spectral planck function
+                                        ! in sub layers
+
+      real tautotal(ndlon,nflev,nir)  ! \   Total single scattering
+      real omegtotal(ndlon,nflev,nir) !  >  properties (Addition of the
+      real gtotal(ndlon,nflev,nir)    ! /   NAERKIND aerosols prop.)
+
+      real newcoolrate(ndlon,nflev) ! cooling rate (K/s) / with implicite scheme
+
+      real emis_gaz(ndlo2)         ! emissivity for gaz computations
+
+      integer jk,jkk,ja,jl
+
+
+c----------------------------------------------------------------------
+c         0.3   Initialisation
+c               --------------
+
+      DO jl=1 , kdlon
+         IF(co2ice(jl) .GT. 20.e-3) THEN
+             emis_gaz(jl)=1. 
+         ELSE
+             emis_gaz(jl)=emis(jl)
+         ENDIF
+      ENDDO
+
+c----------------------------------------------------------------------
+c         1.0   planck function
+c               ---------------
+
+      call lwb ( kdlon, kflev, tlev, tlay, dt0
+     .         , bsurf, btop, blay, blev, dblay, dbsublay)
+
+c----------------------------------------------------------------------
+c         2.0   absorber amounts
+c               ----------------
+
+      call lwu ( kdlon, kflev
+     .         , dp, plev, tlay, aerosol
+     &         , QIRsQREF3d,omegaIR3d,gIR3d
+     .         , aer_t, co2_u, co2_up
+     .         , tautotal,omegtotal,gtotal)
+
+c----------------------------------------------------------------------
+c         3.0   transmission functions / exchange coefficiants
+c               ----------------------------------------------
+
+c                                                                distants
+c                                                                --------
+                    if( mod(icount-1,ilwd).eq.0) then
+
+c     print*, 'CALL of DISTANTS'
+      call lwxd ( ig0, kdlon, kflev, emis_gaz
+     .          , aer_t, co2_u, co2_up)
+
+                    endif
+c                                                              neighbours
+c                                                              ----------
+                    if( mod(icount-1,ilwn).eq.0) then
+
+c     print*, 'CALL of NEIGHBOURS'
+      call lwxn ( ig0, kdlon, kflev
+     .          , dp
+     .          , aer_t, co2_u, co2_up)
+
+                    endif
+c                                                              boundaries
+c                                                              ----------
+                    if( mod(icount-1,ilwb).eq.0) then
+
+c     print*, 'CALL of BOUNDARIES'
+      call lwxb ( ig0, kdlon, kflev, emis_gaz
+     .          , aer_t, co2_u, co2_up)
+
+                    endif
+
+c----------------------------------------------------------------------
+c         4.0   cooling rate
+c               ------------
+
+      call lwflux ( ig0, kdlon, kflev, dp
+     .            , bsurf, btop, blev, blay, dbsublay
+     .            , tlay, tlev, dt0      ! pour sortie dans g2d uniquement
+     .            , emis
+     .            , tautotal,omegtotal,gtotal
+     .            , coolrate, fluxground, fluxtop
+     .            , netrad)
+
+c     do jk = 1, nlaylte
+c       print*,coolrate(1,jk)
+c     enddo
+      
+c       do jkk = 0 , nlaylte+1
+c         do jk = 0 , nlaylte+1
+c           do ja = 1 , nuco2
+c             do jl = 1 , ngrid
+c      if (xi (jl,ja,jk,jkk) .LT. 0
+c    .       .OR. xi (jl,ja,jk,jkk) .GT. 1 ) then
+c                 print*,'xi bande',ja,jk,jkk,xi (jl,ja,jk,jkk)
+c      endif
+c             enddo
+c           enddo
+c         enddo
+c       enddo
+
+c----------------------------------------------------------------------
+c
+c          5.    shema semi-implicite  (lwi) 
+c                ---------------------------
+c
+c
+      call lwi (ig0,kdlon,kflev,netrad,dblay,dp
+     .          , newcoolrate)
+c
+c  Verif que   (X sol,space) + somme(X i,sol) = 1
+c
+      do jkk = 1 , nlaylte
+        do jl = 1 , kdlon
+c     print*,'NEW et OLD coolrate :',jkk,newcoolrate(jl,jkk) 
+c    .  ,coolrate(jl,jkk)
+      coolrate(jl,jkk) = newcoolrate(jl,jkk)
+        enddo
+      enddo
+c
+c----------------------------------------------------------------------
+
+      END SUBROUTINE lwmain
+
+      END MODULE lwmain_mod 
+
Index: unk/LMDZ.MARS/libf/phymars/swmain.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/swmain.F	(revision 1982)
+++ 	(revision )
@@ -1,186 +1,0 @@
-      SUBROUTINE SWMAIN ( KDLON, KFLEV,
-     $                PCST, albedo,
-     $                PRMU0, PDP, PPLEV, aerosol,PFRACT,
-     $                PHEAT, PFLUXD,PFLUXU,
-     &                QVISsQREF3d,omegaVIS3d,gVIS3d)
-
-      use dimradmars_mod, only: ndlo2, ndlon, nflev, 
-     &                          nsun,naerkind
-      use yomlw_h, only: nlaylte, gcp
-      IMPLICIT NONE
-      
-#include "callkeys.h"
-c     
-c     PURPOSE.
-c     --------
-c     
-c     This routine computes the shortwave (solar wavelength)
-c     radiation fluxes in two spectral intervals 
-c     and heating rate on the first "nlaylte" layers.
-C     
-c     Francois Forget (2000), adapted from 
-C     Fouquart and Bonnel's ECMWF program
-c 
-C     IMPLICIT ARGUMENTS :
-C     --------------------
-C     
-C     ==== INPUTS ===
-c
-c    KDLON :  number of horizontal grid points
-c    PCST  :  Solar constant on Mars (W.m-2)
-c    albedo   hemispheric surface albedo
-c                         albedo (i,1) : mean albedo for solar band#1
-c                                        (see below)
-c                         albedo (i,2) : mean albedo for solar band#2
-c                                        (see below)
-c    PRMU0 :         cos of solar zenith angle (=1 when sun at zenith)
-c    PDP :           Layer thickness (Pa)
-c    PPLEV           pressure (Pa) at boundaries of each layer 
-c   aerosol               aerosol extinction optical depth
-c                         at reference wavelength "longrefvis" set
-c                         in dimradmars_mod , in each layer, for one of
-c                         the "naerkind" kind of aerosol optical properties.
-c    Pfract :        day fraction of the time interval
-c                          =1 during the full day ; =0 during the night
-c    QVISsQREF3d,omegaVIS3d,gVIS3d Aerosol optical properties
-c
-C     ==== OUTPUTS ===
-c    PHEAT :         Heating rate (K/s)
-c    PFLUXD :        SW downward flux at boundaries of each layer (W.m-2)
-c    PFLUXU :        SW upward flux at boundaries of each layer (W.m-2)
-C     
-C     ----------
-C     
-C-----------------------------------------------------------------------
-C     
-C     
-C-----------------------------------------------------------------------
-C     
-C     ARGUMENTS
-C     ---------
-      
-      INTEGER KDLON, KFLEV
-      REAL ZPSOL(NDLO2), aerosol(NDLO2,KFLEV,naerkind),PRMU0(NDLO2)
-      real PCST
-      REAL albedo(NDLO2,2)
-      REAL PDP(NDLO2,KFLEV)
-      REAL PPLEV(NDLO2,KFLEV+1)
-      REAL PHEAT(NDLO2,KFLEV)
-      REAL PFRACT(NDLO2)
-      real PFLUXD(NDLON,NFLEV+1,2)
-      real PFLUXU(NDLON,NFLEV+1,2)
-      REAL :: QVISsQREF3d(NDLO2,KFLEV,nsun,naerkind)
-      REAL :: omegaVIS3d(NDLO2,KFLEV,nsun,naerkind)
-      REAL :: gVIS3d(NDLO2,KFLEV,nsun,naerkind)
-      
-C     LOCAL ARRAYS
-C     ------------
-      
-      REAL ZDSIG(NDLON,NFLEV), ZFACT(NDLON)
-     S     ,  ZFD(NDLON,NFLEV+1)
-     S     ,  ZFU(NDLON,NFLEV+1)
-     S     ,  ZRMU(NDLON), ZSEC(NDLON)
-     S     ,  ZUD(NDLON,3,NFLEV+1), ZUM(NDLON,NFLEV+1)
-      REAL ZSIGN(NDLON),  ZSIGO(NDLON)
-
-c following line has been changed, kflev--->nflev (to avoid error message 
-c when compiling on NASA Ames Sun)
-      REAL ZFDOWN(NDLO2,NFLEV+1),ZFUP(NDLO2,NFLEV+1)
-
-      integer jl, jk, jkp1, jkl
-      integer INU 
-      real  zdfnet
-
-C     ------------------------------------------------------------------
-C     Initializations :
-C     ------------------------------------------------------------------
-
-c     Incident Solar flux and corrected angle in the atmosphere 
-c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      DO JL = 1 , KDLON
-c        Incident Flux at the top of the atmosphere
-         ZFACT(JL)= PRMU0(JL) * PCST * PFRACT(JL)
-
-c        Cos of solar zenith angle CORRECTED for high zenith angle
-       if (PRMU0(JL).GT.0) then
-         ZRMU(JL)=SQRT(1224.* PRMU0(JL) * PRMU0(JL) + 1.) / 35.
-       else
-         ZRMU(JL)= 1. / 35.
-       endif
-         ZSEC(JL)=1./ZRMU(JL)
-      END DO
-
-c     Calcul of ZDSIG  (thickness of layers in sigma coordinates)
-c     ~~~~~~~~~~~~~~~
-      DO JL = 1 , KDLON
-         ZSIGO(JL) = 1.0
-         ZPSOL(JL) = PPLEV(JL,1)
-      END DO
-      DO JK = 1 , nlaylte
-         JKP1 = JK + 1
-         JKL = nlaylte+1 - JK
-         DO  JL = 1 , KDLON
-            ZSIGN(JL) = PPLEV(JL,JKP1) / PPLEV(JL,1)
-            ZDSIG(JL,JK) = ZSIGO(JL) - ZSIGN(JL)
-            ZSIGO(JL) = ZSIGN(JL)
-         END DO
-      END DO
-
-C------------------------------------------------------------------
-C          LOOP ON SPECTRAL INTERVAL in solar spectrum
-C------------------------------------------------------------------
-c  2 spectral interval in solar spectrum :
-c  - INU=1: between wavelength "long1vis" and "long2vis" set in dimradmars_mod
-c  - INU=2: between wavelength "long2vis" and "long3vis" set in dimradmars_mod
-
-      DO INU = 1,2
-
-! NB: swrtype is set in callkeys.h
-        if (swrtype.eq.1) then ! Fouquart
-          CALL SWR_FOUQUART( KDLON, kflev, INU
-     &     ,  aerosol,QVISsQREF3d,omegaVIS3d,gVIS3d
-     &     ,  albedo,ZDSIG,ZPSOL,ZRMU,ZSEC
-     &     ,  ZFD,ZFU )
-        else
-          if (swrtype.eq.2) then ! Toon
-            CALL SWR_TOON( KDLON, kflev, INU
-     &       ,  aerosol,QVISsQREF3d,omegaVIS3d,gVIS3d
-     &       ,  albedo,ZDSIG,ZPSOL,ZRMU,ZSEC
-     &       ,  ZFD,ZFU )
-          else
-            write(*,*) "swmain: invalid swrtype value !!"
-            stop
-          endif ! of if (swrtype.eq.2)
-        endif ! of if (swrtype.eq.1)
-        
-         DO JK = 1 , nlaylte+1
-            DO  JL = 1 , KDLON
-              PFLUXD(JL,JK,INU)=ZFD(JL,JK)*ZFACT(JL)
-              PFLUXU(JL,JK,INU)=ZFU(JL,JK)*ZFACT(JL)
-            END DO
-         END DO
-      END DO ! of DO INU=1,2
-
-C     ------------------------------------------------------
-C     HEATING RATES
-C     ------------------------------------------------------
-
-      DO JK = 1 , nlaylte+1
-         DO  JL = 1 , KDLON
-c           wavelength integrated flux at every level:
-            ZFUP(JL,JK)= (PFLUXU(JL,JK,1)+ PFLUXU(JL,JK,2))
-            ZFDOWN(JL,JK)= (PFLUXD(JL,JK,1)+ PFLUXD(JL,JK,2))
-         END DO
-      END DO
-     
-      DO JK = 1 , nlaylte
-         DO  JL = 1 , KDLON
-            ZDFNET = ZFUP (JL,JK  ) - ZFDOWN(JL,JK  )
-     S              -ZFUP (JL,JK+1) + ZFDOWN(JL,JK+1)
-c           Heating rate
-            PHEAT(JL,JK) = gcp * ZDFNET / PDP(JL,JK)
-         END DO
-      END DO
-
-      RETURN
-      END
Index: /trunk/LMDZ.MARS/libf/phymars/swmain_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/swmain_mod.F	(revision 1983)
+++ /trunk/LMDZ.MARS/libf/phymars/swmain_mod.F	(revision 1983)
@@ -0,0 +1,197 @@
+      MODULE swmain_mod
+      
+      IMPLICIT NONE
+
+      CONTAINS
+
+      SUBROUTINE SWMAIN ( KDLON, KFLEV,
+     $                PCST, albedo,
+     $                PRMU0, PDP, PPLEV, aerosol,PFRACT,
+     $                PHEAT, PFLUXD,PFLUXU,
+     &                QVISsQREF3d,omegaVIS3d,gVIS3d)
+
+      use dimradmars_mod, only: ndlo2, ndlon, nflev, 
+     &                          nsun,naerkind
+      use yomlw_h, only: nlaylte, gcp
+      IMPLICIT NONE
+
+c     DECLARATIONS
+c     -------------      
+      include "callkeys.h"
+c     
+c     PURPOSE.
+c     --------
+c     
+c     This routine computes the shortwave (solar wavelength)
+c     radiation fluxes in two spectral intervals 
+c     and heating rate on the first "nlaylte" layers.
+C     
+c     Francois Forget (2000), adapted from 
+C     Fouquart and Bonnel's ECMWF program
+c 
+C     IMPLICIT ARGUMENTS :
+C     --------------------
+C     
+C     ==== INPUTS ===
+c
+c    KDLON :  number of horizontal grid points
+c    PCST  :  Solar constant on Mars (W.m-2)
+c    albedo   hemispheric surface albedo
+c                         albedo (i,1) : mean albedo for solar band#1
+c                                        (see below)
+c                         albedo (i,2) : mean albedo for solar band#2
+c                                        (see below)
+c    PRMU0 :         cos of solar zenith angle (=1 when sun at zenith)
+c    PDP :           Layer thickness (Pa)
+c    PPLEV           pressure (Pa) at boundaries of each layer 
+c   aerosol               aerosol extinction optical depth
+c                         at reference wavelength "longrefvis" set
+c                         in dimradmars_mod , in each layer, for one of
+c                         the "naerkind" kind of aerosol optical properties.
+c    Pfract :        day fraction of the time interval
+c                          =1 during the full day ; =0 during the night
+c    QVISsQREF3d,omegaVIS3d,gVIS3d Aerosol optical properties
+c
+C     ==== OUTPUTS ===
+c    PHEAT :         Heating rate (K/s)
+c    PFLUXD :        SW downward flux at boundaries of each layer (W.m-2)
+c    PFLUXU :        SW upward flux at boundaries of each layer (W.m-2)
+C     
+C     ----------
+C     
+C-----------------------------------------------------------------------
+C     
+C     
+C-----------------------------------------------------------------------
+C     
+C     ARGUMENTS
+C     ---------
+c     INPUTS/OUTPUTS:
+c     ---------
+      
+      INTEGER, iNTENT(IN) :: KDLON, KFLEV
+      REAL, iNTENT(IN) :: aerosol(NDLO2,KFLEV,naerkind),PRMU0(NDLO2)
+      REAL, iNTENT(IN) :: PCST
+      REAL, iNTENT(IN) :: albedo(NDLO2,2)
+      REAL, iNTENT(IN) :: PDP(NDLO2,KFLEV)
+      REAL, iNTENT(IN) :: PPLEV(NDLO2,KFLEV+1)
+      REAL, iNTENT(OUT) :: PHEAT(NDLO2,KFLEV)
+      REAL, iNTENT(IN) :: PFRACT(NDLO2)
+      REAL, iNTENT(OUT) :: PFLUXD(NDLON,NFLEV+1,2)
+      REAL, iNTENT(OUT) :: PFLUXU(NDLON,NFLEV+1,2)
+      REAL, iNTENT(IN) :: QVISsQREF3d(NDLO2,KFLEV,nsun,naerkind)
+      REAL, iNTENT(IN) :: omegaVIS3d(NDLO2,KFLEV,nsun,naerkind)
+      REAL, iNTENT(IN) :: gVIS3d(NDLO2,KFLEV,nsun,naerkind)
+      
+C     LOCAL ARRAYS
+C     ------------
+      REAL ZPSOL(NDLO2)
+      REAL ZDSIG(NDLON,NFLEV), ZFACT(NDLON)
+     S     ,  ZFD(NDLON,NFLEV+1)
+     S     ,  ZFU(NDLON,NFLEV+1)
+     S     ,  ZRMU(NDLON), ZSEC(NDLON)
+     S     ,  ZUD(NDLON,3,NFLEV+1), ZUM(NDLON,NFLEV+1)
+      REAL ZSIGN(NDLON),  ZSIGO(NDLON)
+
+c following line has been changed, kflev--->nflev (to avoid error message 
+c when compiling on NASA Ames Sun)
+      REAL ZFDOWN(NDLO2,NFLEV+1),ZFUP(NDLO2,NFLEV+1)
+
+      integer jl, jk, jkp1, jkl
+      integer INU 
+      real  zdfnet
+
+C     ------------------------------------------------------------------
+C     Initializations :
+C     ------------------------------------------------------------------
+
+c     Incident Solar flux and corrected angle in the atmosphere 
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      DO JL = 1 , KDLON
+c        Incident Flux at the top of the atmosphere
+         ZFACT(JL)= PRMU0(JL) * PCST * PFRACT(JL)
+
+c        Cos of solar zenith angle CORRECTED for high zenith angle
+       if (PRMU0(JL).GT.0) then
+         ZRMU(JL)=SQRT(1224.* PRMU0(JL) * PRMU0(JL) + 1.) / 35.
+       else
+         ZRMU(JL)= 1. / 35.
+       endif
+         ZSEC(JL)=1./ZRMU(JL)
+      END DO
+
+c     Calcul of ZDSIG  (thickness of layers in sigma coordinates)
+c     ~~~~~~~~~~~~~~~
+      DO JL = 1 , KDLON
+         ZSIGO(JL) = 1.0
+         ZPSOL(JL) = PPLEV(JL,1)
+      END DO
+      DO JK = 1 , nlaylte
+         JKP1 = JK + 1
+         JKL = nlaylte+1 - JK
+         DO  JL = 1 , KDLON
+            ZSIGN(JL) = PPLEV(JL,JKP1) / PPLEV(JL,1)
+            ZDSIG(JL,JK) = ZSIGO(JL) - ZSIGN(JL)
+            ZSIGO(JL) = ZSIGN(JL)
+         END DO
+      END DO
+
+C------------------------------------------------------------------
+C          LOOP ON SPECTRAL INTERVAL in solar spectrum
+C------------------------------------------------------------------
+c  2 spectral interval in solar spectrum :
+c  - INU=1: between wavelength "long1vis" and "long2vis" set in dimradmars_mod
+c  - INU=2: between wavelength "long2vis" and "long3vis" set in dimradmars_mod
+
+      DO INU = 1,2
+
+! NB: swrtype is set in callkeys.h
+        if (swrtype.eq.1) then ! Fouquart
+          CALL SWR_FOUQUART( KDLON, kflev, INU
+     &     ,  aerosol,QVISsQREF3d,omegaVIS3d,gVIS3d
+     &     ,  albedo,ZDSIG,ZPSOL,ZRMU,ZSEC
+     &     ,  ZFD,ZFU )
+        else
+          if (swrtype.eq.2) then ! Toon
+            CALL SWR_TOON( KDLON, kflev, INU
+     &       ,  aerosol,QVISsQREF3d,omegaVIS3d,gVIS3d
+     &       ,  albedo,ZDSIG,ZPSOL,ZRMU,ZSEC
+     &       ,  ZFD,ZFU )
+          else
+            write(*,*) "swmain: invalid swrtype value !!"
+            stop
+          endif ! of if (swrtype.eq.2)
+        endif ! of if (swrtype.eq.1)
+        
+         DO JK = 1 , nlaylte+1
+            DO  JL = 1 , KDLON
+              PFLUXD(JL,JK,INU)=ZFD(JL,JK)*ZFACT(JL)
+              PFLUXU(JL,JK,INU)=ZFU(JL,JK)*ZFACT(JL)
+            END DO
+         END DO
+      END DO ! of DO INU=1,2
+
+C     ------------------------------------------------------
+C     HEATING RATES
+C     ------------------------------------------------------
+
+      DO JK = 1 , nlaylte+1
+         DO  JL = 1 , KDLON
+c           wavelength integrated flux at every level:
+            ZFUP(JL,JK)= (PFLUXU(JL,JK,1)+ PFLUXU(JL,JK,2))
+            ZFDOWN(JL,JK)= (PFLUXD(JL,JK,1)+ PFLUXD(JL,JK,2))
+         END DO
+      END DO
+     
+      DO JK = 1 , nlaylte
+         DO  JL = 1 , KDLON
+            ZDFNET = ZFUP (JL,JK  ) - ZFDOWN(JL,JK  )
+     S              -ZFUP (JL,JK+1) + ZFDOWN(JL,JK+1)
+c           Heating rate
+            PHEAT(JL,JK) = gcp * ZDFNET / PDP(JL,JK)
+         END DO
+      END DO
+
+      END SUBROUTINE SWMAIN
+
+      END MODULE swmain_mod
