source: trunk/LMDZ.MARS/libf/phymars/callradite_mod.F @ 2800

Last change on this file since 2800 was 2685, checked in by emillour, 2 years ago

Mars GCM:
Add possibility to output either upward or downward SW flux at the surface
and top of atmosphere from physiq. Required adding some output arguments
to callradite.
EM

File size: 24.0 KB
Line 
1      MODULE callradite_mod
2
3      IMPLICIT NONE
4
5      CONTAINS
6
7      SUBROUTINE callradite(icount,ngrid,nlayer,nq,zday,ls,pq,albedo,
8     $     emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
9     $     dtlw,dtsw,fluxsurf_lw,fluxsurf_dn_sw,fluxsurf_up_sw,
10     $     fluxtop_lw,fluxtop_dn_sw,fluxtop_up_sw,
11     &     tau_pref_scenario,tau_pref_gcm,
12     &     tau,aerosol,dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef,
13     $     taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,co2ice,
14     $     rstormdust,rtopdust,totstormfract,clearatm,dsords,dsotop,
15     $     nohmons,clearsky,totcloudfrac)
16
17      use aeropacity_mod, only: aeropacity
18      use updatereffrad_mod, only: updatereffrad
19      use dimradmars_mod, only: ndomainsz, nflev, nsun, nir
20      use dimradmars_mod, only: naerkind, name_iaer,
21     &            iaer_dust_conrath,iaer_dust_doubleq,
22     &            iaer_dust_submicron, iaer_h2o_ice, iaer_co2_ice,
23     &            iaer_stormdust_doubleq,iaer_topdust_doubleq,flux_1AU
24      use yomlw_h, only: gcp, nlaylte
25      use comcstfi_h, only: g,cpp
26      use time_phylmdz_mod, only: daysec
27      use lwmain_mod, only: lwmain
28      use swmain_mod, only: swmain
29      use dust_param_mod, only: doubleq, active, submicron
30      IMPLICIT NONE
31c=======================================================================
32c   subject:
33c   --------
34c   Subroutine designed to call the main canonic
35c   radiative transfer subroutine "lwmain" et "swmain"
36c   to compute radiative heating and cooling rate and
37c   radiative fluxes to the surface.
38c
39c   These calculations are only valid on the part of the atmosphere
40c   where Local Thermal Equilibrium (LTE) is verified. In practice
41c   the calculations are only performed for the first "nlaylte"
42c   parameters (nlaylte is calculated by subroutine "nlthermeq"
43c   and stored in module "yomlw_h").
44c
45c   The purpose of this subroutine is to:
46c      1) Make some initial calculation at first call
47c      2) Split the calculation in several sub-grid
48c        ("sub-domain") to save memory and
49c        be able run on a workstation at high resolution
50c        The sub-grid size is defined in dimradmars_mod
51c      3) Compute the 3D scattering parameters depending on the
52c        size distribution of the different tracers (added by JBM)
53c      4) call "lwmain" and "swmain"
54c
55c
56c   authors:   
57c   ------
58c   Francois Forget / Christophe Hourdin / J.-B. Madeleine (2009)
59c
60c
61c   3D scattering scheme user's guide (J.-B. Madeleine)
62c   ---------------------------------
63c
64c   This routine has been modified to take into account 3D, time
65c   dependent scattering properties of the aerosols.
66c---- The look-up tables that contain the scattering parameters
67c   of a given tracer, for different sizes, are read by SUAER.F90.
68c   The names of the corresponding ASCII files have to be set in
69c   this subroutine (file_id variable), and files must be in the
70c   directory specified in datafile_mod. Please make sure that the
71c   ASCII files are correctly written, and that the range
72c   of particle sizes is consistent with what you would expect.
73c---- SUAER.F90 is in charge of reading the ASCII files and averaging
74c   the scattering parameters in each GCM channel, using the three last
75c   equations of Forget et al. 1998 (GRL 25, No.7, p.1105-1108).
76c---- These look-up tables, loaded during the firstcall, are then
77c   constantly used by the subroutine "aeroptproperties.F" to compute,
78c   online, the 3D scattering parameters, based on the size distribution
79c   (reffrad and nueffrad) of the different tracers, in each grid box.
80c   These 3D size distributions are loaded by the "updatereffrad.F"
81c   subroutine. A log-normal distribution is then assumed in
82c   "aeroptproperties.F", along with a Gauss-Legendre integration.
83c---- The optical depth at the visible reference wavelength (set in
84c   SUAER.F90, after the file_id variable) is then computed by
85c   the subroutine "aeropacity.F", by using the size and spatial
86c   distribution of the corresponding tracer. This connection has to
87c   be implemented in "aeropacity.F" when adding a new tracer. To do so,
88c   one can use equation 2 of Forget et al. 1998 (Icarus 131, p.302-316).
89c---- The resulting variables "aerosol", "QVISsQREF3d", "omegaVIS3d" and
90c   "gVIS3d" (same in the infrared) are finally used by lwmain.F and
91c   swmain.F to solve the radiative transfer equation.
92c
93c   changes:
94c   -------
95c
96c   > SRL 7/2000
97c   
98c   This version has been modified to only calculate radiative tendencies
99c   over layers 1..NFLEV (set in dimradmars_mod).  Returns zero for higher
100c   layers, if any.
101c   In other routines, nlayer -> nflev.
102c   Routines affected: lwflux, lwi, lwmain, lwxb, lwxd, lwxn.
103c
104c   > J.-B. Madeleine 10W12
105c   This version uses the variable's splitting, which can be usefull
106c     when performing very high resolution simulation like LES.
107c
108c   ----------
109c   Here, solar band#1 is spectral interval between "long1vis" and "long2vis"
110c   set in dimradmars_mod
111c   Here, solar band#2 is spectral interval between "long2vis" and "long3vis"
112c   set in dimradmars_mod
113c
114c   input:
115c   -----
116c   icount                counter of call to subroutine physic by gcm
117c   ngrid                 number of gridpoint of horizontal grid
118c   nlayer                Number of layer
119c   nq                    Number of tracer
120c   ls                    Solar longitude (Ls) , radian
121c   zday                  Date (time since Ls=0, in martian days)
122c   pq(ngrid,nlayer,nq)   Advected fields
123c
124c   albedo (ngrid,2)      hemispheric surface albedo
125c                         albedo (i,1) : mean albedo for solar band#1
126c                                        (see below)
127c                         albedo (i,2) : mean albedo for solar band#2
128c                                        (see below)
129c   emis                  Thermal IR surface emissivity (no unit)
130c   mu0(ngrid)           cos of solar zenith angle
131c                           (=1 when sun at zenith)
132c   pplay(ngrid,nlayer)    pressure (Pa) in the middle of each layer
133c   pplev(ngrid,nlayer+1)  pressure (Pa) at boundaries of each layer
134c   pt(ngrid,nlayer)       atmospheric temperature in each layer (K)
135c   tsurf(ngrid)           surface temperature (K)
136c   fract(ngrid)         day fraction of the time interval
137c                          =1 during the full day ; =0 during the night
138c   declin                 latitude of subsolar point
139c   dist_sol               sun-Mars distance (AU)
140c   igout                  coordinate of analysed point for debugging
141c   reffrad(ngrid,nlayer,naerkind)  Aerosol effective radius
142c   nueffrad(ngrid,nlayer,naerkind) Aerosol effective variance
143
144c=======================================================================
145c
146c    Declarations :
147c    -------------
148c
149      include "callkeys.h"
150
151c-----------------------------------------------------------------------
152c    Input/Output
153c    ------------
154      INTEGER,INTENT(IN) :: icount       
155      INTEGER,INTENT(IN) :: ngrid,nlayer,nq
156      INTEGER,INTENT(IN) :: igout
157
158      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)
159      REAL,INTENT(INOUT) :: tauscaling(ngrid) ! Conversion factor for
160                               ! qdust and Ndust
161      REAL,INTENT(INOUT) :: dust_rad_adjust(ngrid) ! Radiative adjustment
162                          ! factor for dust
163      REAL,INTENT(INOUT) :: IRtoVIScoef(ngrid) ! conversion coefficient to apply on
164                                               ! scenario absorption IR (9.3um) CDOD
165                                               ! = tau_pref_gcm_VIS / tau_pref_gcm_IR
166      REAL,INTENT(IN) :: albedo(ngrid,2),emis(ngrid)
167      REAL,INTENT(IN) :: ls,zday
168
169      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
170      REAL,INTENT(IN) :: pt(ngrid,nlayer)
171      REAL,INTENT(IN) :: tsurf(ngrid)
172      REAL,INTENT(IN) :: dist_sol,mu0(ngrid),fract(ngrid)
173      REAL,INTENT(OUT) :: dtlw(ngrid,nlayer) ! longwave (IR) heating rate (K/s)
174      REAL,INTENT(OUT) :: dtsw(ngrid,nlayer) ! shortwave (Solar) heating rate (K/s)
175      REAL,INTENT(OUT) :: fluxsurf_lw(ngrid) ! total LW (thermal IR) downward flux
176                                             ! (W.m-2) at the surface
177      REAL,INTENT(OUT) :: fluxtop_lw(ngrid) ! outgoing total LW (thermal IR)
178                                     ! upward flux (W.m-2) at the top of the atm.
179      REAL,INTENT(OUT) :: fluxsurf_dn_sw(ngrid,2) ! surface downward SW flux for
180                                      ! solar bands #1 and #2 (W.m-2)
181      REAL,INTENT(OUT) :: fluxsurf_up_sw(ngrid,2) ! surface upward SW flux for
182                                      ! solar bands #1 and #2 (W.m-2)
183      REAL,INTENT(OUT) :: fluxtop_dn_sw(ngrid,2) ! incoming downward SW flux for
184                                 ! solar bands #1 and #2 (W.m-2) at top of atm.
185      REAL,INTENT(OUT) :: fluxtop_up_sw(ngrid,2) ! outgoing upward SW flux for
186                                 ! solar bands #1 and #2 (W.m-2) at top of atm.
187      REAL,INTENT(OUT) :: tau_pref_scenario(ngrid) ! prescribed dust column
188                          ! visible opacity at odpref from scenario
189      REAL,INTENT(OUT) :: tau_pref_gcm(ngrid) ! computed dust column
190                          ! visible opacity at odpref in the GCM
191      REAL,INTENT(OUT) :: tau(ngrid,naerkind) ! Column visible optical depth
192                          ! for each aerosol
193      REAL,INTENT(OUT) :: taucloudtes(ngrid)! Cloud opacity at infrared
194                               !   reference wavelength using
195                               !   Qabs instead of Qext
196                               !   (direct comparison with TES)
197      REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol extinction
198                          ! optical depth at reference wavelength "longrefvis",
199                          ! set in dimradmars_h, for each kind of aerosol
200      REAL,INTENT(INOUT) :: dsodust(ngrid,nlayer)
201      REAL,INTENT(OUT) :: rdust(ngrid,nlayer)  ! Dust geometric mean radius (m)
202      REAL,INTENT(OUT) :: rice(ngrid,nlayer)   ! Ice geometric mean radius (m)
203      REAL,INTENT(OUT) :: nuice(ngrid,nlayer)  ! Estimated effective variance
204      double precision,INTENT(OUT) :: riceco2(ngrid,nlayer) ! CO2 ice mean radius(m)
205      REAL,INTENT(OUT) :: nuiceco2(ngrid,nlayer) ! Effective variance
206      REAL,INTENT(IN) :: co2ice(ngrid)           ! co2 ice surface layer (kg.m-2)
207
208c     rocket dust storm
209      LOGICAL,INTENT(IN) :: clearatm ! true for background dust
210      REAL,INTENT(IN) :: totstormfract(ngrid) ! dust storm mesh fraction
211      REAL,INTENT(OUT) :: rstormdust(ngrid,nlayer)  ! Storm dust geometric mean radius (m)
212      REAL,INTENT(OUT) :: dsords(ngrid,nlayer) ! density scaled opacity for rocket dust storm dust
213     
214c     entrainment by mountain top dust flows
215      LOGICAL, INTENT(IN) :: nohmons ! true for background dust
216      REAL,INTENT(OUT) :: rtopdust(ngrid,nlayer)  ! Topdust geometric mean radius (m)
217      REAL,INTENT(OUT) :: dsotop(ngrid,nlayer) ! density scaled opacity for topmons dust
218     
219c     sub-grid scale water ice clouds
220      LOGICAL,INTENT(IN) :: clearsky
221      REAL,INTENT(IN) :: totcloudfrac(ngrid)
222
223c
224c    Local variables :
225c    -----------------
226
227      INTEGER j,l,ig,n,ich
228      INTEGER aer_count,iaer
229      INTEGER jd,ig0,nd
230
231      real  cste_mars ! solar constant on Mars (Wm-2)
232      REAL ptlev(ngrid,nlayer+1)
233
234      INTEGER :: ndomain
235
236c     Thermal IR net radiative budget (W m-2)
237
238      real znetrad(ndomainsz,nflev)
239
240      real zfluxd_sw(ndomainsz,nflev+1,2)
241      real zfluxu_sw(ndomainsz,nflev+1,2)
242
243      REAL zplev(ndomainsz,nflev+1)
244      REAL zztlev(ndomainsz,nflev+1)
245      REAL zplay(ndomainsz,nflev)
246      REAL zt(ndomainsz,nflev)
247      REAL zaerosol(ndomainsz,nflev,naerkind)
248      REAL zalbedo(ndomainsz,2)
249      REAL zdp(ndomainsz,nflev)
250      REAL zdt0(ndomainsz)
251
252      REAL zzdtlw(ndomainsz,nflev)
253      REAL zzdtsw(ndomainsz,nflev)
254      REAL zzflux(ndomainsz,6)
255      real zrmuz
256
257      REAL :: zQVISsQREF3d(ndomainsz,nflev,nsun,naerkind)
258      REAL :: zomegaVIS3d(ndomainsz,nflev,nsun,naerkind)
259      REAL :: zgVIS3d(ndomainsz,nflev,nsun,naerkind)
260
261      REAL :: zQIRsQREF3d(ndomainsz,nflev,nir,naerkind)
262      REAL :: zomegaIR3d(ndomainsz,nflev,nir,naerkind)
263      REAL :: zgIR3d(ndomainsz,nflev,nir,naerkind)
264
265c     Aerosol size distribution
266      REAL :: reffrad(ngrid,nlayer,naerkind)
267      REAL :: nueffrad(ngrid,nlayer,naerkind)
268c     Aerosol optical properties
269      REAL :: QVISsQREF3d(ngrid,nlayer,nsun,naerkind)
270      REAL :: omegaVIS3d(ngrid,nlayer,nsun,naerkind)
271      REAL :: gVIS3d(ngrid,nlayer,nsun,naerkind)
272
273      REAL :: QIRsQREF3d(ngrid,nlayer,nir,naerkind)
274      REAL :: omegaIR3d(ngrid,nlayer,nir,naerkind)
275      REAL :: gIR3d(ngrid,nlayer,nir,naerkind)
276
277      REAL :: QREFvis3d(ngrid,nlayer,naerkind)
278      ! QREFvis3d : Extinction efficiency at the VISible reference wavelength
279      REAL :: QREFir3d(ngrid,nlayer,naerkind)
280      ! QREFir3d : Extinction efficiency at the InfraRed reference wavelength
281
282      REAL :: omegaREFvis3d(ngrid,nlayer,naerkind)
283      REAL :: omegaREFir3d(ngrid,nlayer,naerkind)
284
285c   local saved variables
286c   ---------------------
287
288      real zco2   ! volume fraction of CO2 in Mars atmosphere
289!$OMP THREADPRIVATE(zco2)
290      DATA zco2/0.95/
291      SAVE zco2
292
293      LOGICAL firstcall
294!$OMP THREADPRIVATE(firstcall)
295      DATA firstcall/.true./
296      SAVE firstcall
297
298
299
300c----------------------------------------------------------------------
301
302c     Initialisation
303c     --------------
304
305! compute ndomain
306! AS: moved out of firstcall to allow nesting+evoluting domain
307! ------------------------------------------------------------
308      ndomain= (ngrid-1) / ndomainsz + 1
309
310      IF (firstcall) THEN
311
312         write(*,*) 'Splitting radiative calculations: ',
313     $              ' ngrid,ndomainsz,ndomain',
314     $                ngrid,ndomainsz,ndomain
315
316c        Assign a number to the different scatterers
317c        -------------------------------------------
318
319         iaer_dust_conrath=0
320         iaer_dust_doubleq=0
321         iaer_dust_submicron=0
322         iaer_h2o_ice=0
323         iaer_co2_ice=0
324         iaer_stormdust_doubleq=0
325         iaer_topdust_doubleq=0
326
327         aer_count=0
328         if (.NOT.active) then
329           do iaer=1,naerkind
330             if (name_iaer(iaer).eq."dust_conrath") then
331               iaer_dust_conrath = iaer
332               aer_count = aer_count + 1
333             endif
334           enddo
335         endif
336         if (doubleq.AND.active) then
337           do iaer=1,naerkind
338             if (name_iaer(iaer).eq."dust_doubleq") then
339               iaer_dust_doubleq = iaer
340               aer_count = aer_count + 1
341             endif
342           enddo
343         endif
344         if (submicron.AND.active) then
345           do iaer=1,naerkind
346             if (name_iaer(iaer).eq."dust_submicron") then
347               iaer_dust_submicron = iaer
348               aer_count = aer_count + 1
349             endif
350           enddo
351         endif
352         if (water.AND.activice) then
353           do iaer=1,naerkind
354             if (name_iaer(iaer).eq."h2o_ice") then
355               iaer_h2o_ice = iaer
356               aer_count = aer_count + 1
357             endif
358           enddo
359         endif
360         if (co2clouds.AND.activeco2ice) then
361           do iaer=1,naerkind
362             if (name_iaer(iaer).eq."co2_ice") then
363               iaer_co2_ice = iaer
364               aer_count = aer_count + 1
365             endif
366           enddo
367         endif
368         if (rdstorm.AND.active) then
369           do iaer=1,naerkind
370             if (name_iaer(iaer).eq."stormdust_doubleq") then
371               iaer_stormdust_doubleq = iaer
372               aer_count = aer_count + 1
373             endif
374           enddo
375         end if
376         if (topflows.AND.active) then
377           do iaer=1,naerkind
378             if (name_iaer(iaer).eq."topdust_doubleq") then
379               iaer_topdust_doubleq = iaer
380               aer_count = aer_count + 1
381             endif
382           enddo
383         end if
384
385c        Check that we identified all tracers:
386         if (aer_count.ne.naerkind) then
387           write(*,*) "callradite: found only ",aer_count," scatterers"
388           write(*,*) "               expected ",naerkind
389           write(*,*) "please make sure that the number of"
390           write(*,*) "scatterers in scatterers.h, the names"
391           write(*,*) "in callradite.F, and the flags in"
392           write(*,*) "callphys.def are all consistent!"
393           do iaer=1,naerkind
394             write(*,*)'      ',iaer,' ',trim(name_iaer(iaer))
395           enddo
396           call abort_physic("callradite","incoherent scatterers",1)
397         else
398           write(*,*) "callradite: found all scatterers, namely:"
399           do iaer=1,naerkind
400             write(*,*)'      ',iaer,' ',trim(name_iaer(iaer))
401           enddo
402         endif
403c        -------------------------------------------
404
405         gcp = g/cpp
406
407
408c        Loading the optical properties in external look-up tables:
409
410         CALL SUAER
411         
412!         CALL SULW ! this step is now done in ini_yomlw_h
413
414         if (ngrid .EQ. 1) then
415           if (ndomainsz .NE. 1) then
416             print*
417             print*,'ATTENTION !!!'
418             print*,'pour tourner en 1D, '
419             print*,'fixer ndomainsz=1 dans phymars/dimradmars_h'
420             print*
421             call exit(1)
422           endif
423         endif
424
425         firstcall=.false.
426      END IF
427
428c     Computing aerosol optical properties and opacity
429c     ------------------------------------------------
430c     Updating aerosol size distributions:
431      CALL updatereffrad(ngrid,nlayer,
432     &                rdust,rstormdust,rtopdust,rice,nuice,
433     &                reffrad,nueffrad, riceco2, nuiceco2,
434     &                pq,tauscaling,tau,pplay, pt)
435c     Computing 3D scattering parameters:
436      gVIS3d(:,:,:,:) = 0.
437      CALL aeroptproperties(ngrid,nlayer,reffrad,nueffrad,
438     &                      QVISsQREF3d,omegaVIS3d,gVIS3d,
439     &                      QIRsQREF3d,omegaIR3d,gIR3d,
440     &                      QREFvis3d,QREFir3d,
441     &                      omegaREFvis3d,omegaREFir3d)
442c     Computing aerosol optical depth in each layer:
443      CALL aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
444     &    pq,pt,tauscaling,dust_rad_adjust,IRtoVIScoef,
445     &    tau_pref_scenario,tau_pref_gcm,tau,taucloudtes,
446     &    aerosol,dsodust,reffrad,
447     &    QREFvis3d,QREFir3d,omegaREFir3d,
448     &    totstormfract,clearatm,dsords,dsotop,
449     &    nohmons,clearsky,totcloudfrac)
450     
451c     Starting loop on sub-domain
452c     ----------------------------
453      zgVIS3d(:,:,:,:) = 0.
454      zfluxd_sw(:,:,:) = 0.
455      zfluxu_sw(:,:,:) = 0.
456      zQVISsQREF3d(:,:,:,:) = 0.
457      zomegaVIS3d(:,:,:,:) = 0.
458      DO jd=1,ndomain
459        ig0=(jd-1)*ndomainsz
460        if (jd.eq.ndomain) then
461         nd=ngrid-ig0
462        else
463         nd=ndomainsz
464        endif
465
466c       Spliting input variable in sub-domain input variables
467c       ---------------------------------------------------
468
469        do l=1,nlaylte
470         do ig = 1,nd
471           do iaer = 1, naerkind
472             do ich = 1, nsun
473               zQVISsQREF3d(ig,l,ich,iaer) =
474     &                           QVISsQREF3d(ig0+ig,l,ich,iaer)
475               zomegaVIS3d(ig,l,ich,iaer) =
476     &                           omegaVIS3d(ig0+ig,l,ich,iaer)
477               zgVIS3d(ig,l,ich,iaer) =
478     &                           gVIS3d(ig0+ig,l,ich,iaer)
479             enddo
480             do ich = 1, nir
481               zQIRsQREF3d(ig,l,ich,iaer) =
482     &                           QIRsQREF3d(ig0+ig,l,ich,iaer)
483               zomegaIR3d(ig,l,ich,iaer) =
484     &                           omegaIR3d(ig0+ig,l,ich,iaer)
485               zgIR3d(ig,l,ich,iaer) =
486     &                           gIR3d(ig0+ig,l,ich,iaer)
487             enddo
488           enddo
489         enddo
490        enddo
491        zplev(:,:) = 0.
492        do l=1,nlaylte+1
493         do ig = 1,nd
494          zplev(ig,l) = pplev(ig0+ig,l)
495         enddo
496        enddo
497        zdp(:,:) = 0.
498       
499        do l=1,nlaylte
500         do ig = 1,nd
501          zplay(ig,l) = pplay(ig0+ig,l)
502          zt(ig,l) = pt(ig0+ig,l)
503c         Thickness of each layer (Pa) :
504          zdp(ig,l)= pplev(ig0+ig,l) - pplev(ig0+ig,l+1)
505         enddo
506        enddo
507        zaerosol(:,:,:) = 0.
508        do n=1,naerkind
509          do l=1,nlaylte
510            do ig=1,nd
511              zaerosol(ig,l,n) = aerosol(ig0+ig,l,n)
512            enddo
513          enddo
514        enddo
515        zalbedo(:,:) = 0.
516        do j=1,2
517          do ig = 1,nd
518           zalbedo(ig,j) = albedo(ig0+ig,j)
519          enddo
520        enddo
521
522c       Intermediate  levels: (computing tlev)
523c       ---------------------------------------
524c       Extrapolation for the air temperature above the surface
525        DO ig=1,nd
526              zztlev(ig,1)=zt(ig,1)+
527     s        (zplev(ig,1)-zplay(ig,1))*
528     s        (zt(ig,1)-zt(ig,2))/(zplay(ig,1)-zplay(ig,2))
529
530              zdt0(ig) = tsurf(ig0+ig) - zztlev(ig,1)
531        ENDDO
532
533        DO l=2,nlaylte
534         DO ig=1, nd
535               zztlev(ig,l)=0.5*(zt(ig,l-1)+zt(ig,l))
536         ENDDO
537        ENDDO
538
539        DO ig=1, nd
540           zztlev(ig,nlaylte+1)=zt(ig,nlaylte)
541        ENDDO
542
543
544c       Longwave ("lw") radiative transfer (= thermal infrared)
545c       -------------------------------------------------------
546        call lwmain (ig0,icount,nd,nflev
547     .        ,zdp,zdt0,emis(ig0+1),zplev,zztlev,zt
548     .        ,zaerosol,zzdtlw
549     .        ,fluxsurf_lw(ig0+1),fluxtop_lw(ig0+1)
550     .        ,znetrad
551     &        ,zQIRsQREF3d,zomegaIR3d,zgIR3d
552     &        ,co2ice(ig0+1))
553
554c       Shortwave ("sw") radiative transfer (= solar radiation)
555c       -------------------------------------------------------
556c          Mars solar constant (W m-2)
557c          flux_1AU = 1370 W.m-2 is the solar constant at 1 AU.
558           cste_mars=flux_1AU/(dist_sol*dist_sol)
559           zzdtsw(:,:) = 0.
560           call swmain ( nd, nflev,
561     S     cste_mars, zalbedo,
562     S     mu0(ig0+1), zdp, zplev, zaerosol, fract(ig0+1),
563     S     zzdtsw, zfluxd_sw, zfluxu_sw,
564     &     zQVISsQREF3d,zomegaVIS3d,zgVIS3d)
565c       ------------------------------------------------------------
566c       Un-spliting output variable from sub-domain input variables
567c       ------------------------------------------------------------
568
569        do l=1,nlaylte
570         do ig = 1,nd
571          dtlw(ig0+ig,l) = zzdtlw(ig,l)
572          dtsw(ig0+ig,l) = zzdtsw(ig,l)
573         enddo
574        enddo
575
576        ptlev(:, :) = 0.
577        do l=1,nlaylte+1
578         do ig = 1,nd
579          ptlev(ig0+ig,l) = zztlev(ig,l)
580         enddo
581        enddo
582
583        ! copy SW fluxes at surface and TOA
584        do ig = 1,nd
585          ! surface downward SW flux
586          fluxsurf_dn_sw(ig0+ig,1) = zfluxd_sw(ig,1,1)
587          fluxsurf_dn_sw(ig0+ig,2) = zfluxd_sw(ig,1,2)
588          ! surface upward SW flux
589          fluxsurf_up_sw(ig0+ig,1) = zfluxu_sw(ig,1,1)
590          fluxsurf_up_sw(ig0+ig,2) = zfluxu_sw(ig,1,2)
591          ! downward SW flux at top of atmosphere
592          fluxtop_dn_sw(ig0+ig,1) = zfluxd_sw(ig,nlaylte+1,1)
593          fluxtop_dn_sw(ig0+ig,2) = zfluxd_sw(ig,nlaylte+1,2)         
594          ! upward SW flux at top of atmosphere
595          fluxtop_up_sw(ig0+ig,1) = zfluxu_sw(ig,nlaylte+1,1)
596          fluxtop_up_sw(ig0+ig,2) = zfluxu_sw(ig,nlaylte+1,2)
597        enddo
598
599      ENDDO         !   (boucle jd=1, ndomain)
600
601c     Zero tendencies for any remaining layers between nlaylte and nlayer
602      if (nlayer.gt.nlaylte) then
603         do l = nlaylte+1, nlayer
604            do ig = 1, ngrid
605               dtlw(ig, l) = 0.
606               dtsw(ig, l) = 0.
607            enddo
608         enddo
609      endif
610c     Output for debugging if lwrite=T
611c     --------------------------------
612c     Write all nlayer layers, even though only nlaylte layers may have
613c     non-zero tendencies.
614
615         IF(lwrite) THEN
616            PRINT*,'Diagnotique for the radiation'
617            PRINT*,'albedo, emissiv, mu0,fract,fluxsurf_lw,fluxsurf_sw'
618            PRINT*,albedo(igout,1),emis(igout),mu0(igout),
619     s           fract(igout), fluxsurf_lw(igout),
620     $     fluxsurf_dn_sw(igout,1)+fluxsurf_dn_sw(igout,2)
621            PRINT*,'Tlay Tlev Play Plev dT/dt SW dT/dt LW (K/s)'
622            PRINT*,'daysec',daysec
623            DO l=1,nlayer
624               PRINT*,pt(igout,l),ptlev(igout,l),
625     s         pplay(igout,l),pplev(igout,l),
626     s         dtsw(igout,l),dtlw(igout,l)
627            ENDDO
628         ENDIF
629
630
631      END SUBROUTINE callradite
632
633      END MODULE callradite_mod
Note: See TracBrowser for help on using the repository browser.