source: trunk/LMDZ.GENERIC/libf/phystd/callsedim.F @ 848

Last change on this file since 848 was 787, checked in by aslmd, 12 years ago

LMDZ.GENERIC. (Sorry for long text but this is a quite major commit)

Paving the path for parallel computations. And moving towards a more flexible code.

Automatic allocation is used within all routines in phystd. No further mention to ngridmx and nqmx.

  1. ngridmx and nqmx are still used in LMDZ.GENERIC in the dyn3d part
  2. if the LMDZ4/LMDZ5 dynamical core is used, there is no more fixed dimensions ngridmx and nqmx --> a fully flexible parallel implementation is now possible (e.g. no need to recompile when changing numbers of processors)

The important stuff :

  • Compilation checked with ifort. OK with and without debug mode. No errors. Checked for: gcm, newstart, rcm1d, kcm1d
  • RUN GCM: Running an Earth test case. Comparison with previous revision --> debug mode : perfect match. bit by bit (diff command). checked with plots --> O1 mode : close match (checked with plots) --> O2 mode : sometimes up to 0.5 K departure.... BUT in this new version O2 and O1 are quite close while in previous version O1 and O2 differed by about, well, typically 0.5 K (pictures available on request)
  • RUN NEWSTART : perfect match (bit-by-bit) in either debug or normal mode.
  • RUN RCM1D : perfect match in normal mode.
  • RUN KCM1D : not tested (I don't know what is the use of kcm1d)

List of main changes :

  • Additional arguments to some subroutines (ngrid and nq)
  • F77 include strategy is obsolete and replaced by F90 module strategy In this new strategy arrays are allocatable and allocated once at first use This has to be done for all common featuring arrays defined with ngridmx or nqmx

surfdat.h >> surfdat_h.F90
tracer.h >> tracer_h.F90
comsaison.h >> comsaison_h.F90
comgeomfi.h >> comgeomfi_h.F90
comsoil.h >> comsoil_h.F90
comdiurn.h >> comdiurn_h.F90
fisice.h >> DELETED. was not used. probably a fossil.
watercap.h >> DELETED. variable put in surfdat_h.F90

  • F77 'save' strategy is obsolete and replaced by F90 'allocatable save' strategy (see previous point and e.g. new version of physiq.F90)
  • Suppressing any mention to advtrac.h which is a common in the dynamics and needs nqmx This was easily solved by adding an argument with tracer names, coming from the dynamics This is probably not a definitive solution, ... but this allows for generic physics to work easily with either LMDZ.GENERIC or LMDZ dynamical cores
  • Removing consistency tests between nq and nqmx ; and ngrid and ngridmx. No use now!
  • Adaptation of rcm1d, kcm1d, newstart given above-mentioned changes

A note on phyetat0 and soil_setting:

  • Now written so that a slice of horizontal size 'ngrid' starting at grid point 'cursor' is read in startfi.nc 'cursor' is defined in dimphys.h and initialized by inifis (or in newstart) this is useful for parallel computations. default behavior is the usual one : sequential runs, cursor is 1, size ngrid is the whole global domain

A note on an additional change :

  • nueffrad is now an argument to callcorrk as is the case for reffrad both are saved in physiq this is for consistency and lisibility (previously nueffrad was saved in callcorrk) ... but there is a call to a function which modifies nueffrad in physiq ... previously this was not modifying nueffrad (although it was quite cumbersome to detect this) ... to be conservative I kept this behaviour and highlighted it with an array nueffrad_dummy ... I added a comment because someone might want to change this
File size: 4.2 KB
Line 
1      SUBROUTINE callsedim(ngrid,nlay, ptimestep,
2     $                pplev,zlev, pt,
3     &                pq, pdqfi, pdqsed,pdqs_sed,nq,rfall)
4
5      USE tracer_h
6
7      IMPLICIT NONE
8
9!==================================================================
10!     
11!     Purpose
12!     -------
13!     Calculates sedimentation of aerosols depending on their
14!     density and radius.
15!     
16!     Authors
17!     -------
18!     F. Forget (1999)
19!     Tracer generalisation by E. Millour (2009)
20!     
21!==================================================================
22
23c-----------------------------------------------------------------------
24c   declarations:
25c   -------------
26
27#include "dimensions.h"
28#include "dimphys.h"
29#include "comcstfi.h"
30#include "callkeys.h"
31
32c   arguments:
33c   ----------
34
35      INTEGER ngrid              ! number of horizontal grid points
36      INTEGER nlay               ! number of atmospheric layers
37      REAL ptimestep             ! physics time step (s)
38      REAL pplev(ngrid,nlay+1)   ! pressure at inter-layers (Pa)
39      REAL pt(ngrid,nlay)        ! temperature at mid-layer (K)
40      REAL zlev(ngrid,nlay+1)    ! altitude at layer boundaries
41
42
43c    Traceurs :
44      integer nq             ! number of tracers
45      real pq(ngrid,nlay,nq)  ! tracers (kg/kg)
46      real pdqfi(ngrid,nlay,nq)  ! tendency before sedimentation (kg/kg.s-1)
47      real pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1)
48      real pdqs_sed(ngrid,nq)    ! flux at surface (kg.m-2.s-1)
49     
50c   local:
51c   ------
52
53      INTEGER l,ig, iq
54
55      real zqi(ngrid,nlayermx,nq) ! to locally store tracers
56      real masse (ngrid,nlayermx) ! Layer mass (kg.m-2)
57      real epaisseur (ngrid,nlayermx) ! Layer thickness (m)
58      real wq(ngrid,nlayermx+1) ! displaced tracer mass (kg.m-2)
59c      real dens(ngrid,nlayermx) ! Mean density of the ice part. accounting for dust core
60      real rfall(ngrid,nlayermx)
61
62
63      LOGICAL firstcall
64      SAVE firstcall
65      DATA firstcall/.true./
66
67c    ** un petit test de coherence
68c       --------------------------
69
70      IF (firstcall) THEN
71        firstcall=.false.
72      ENDIF ! of IF (firstcall)
73     
74!=======================================================================
75!     Preliminary calculation of the layer characteristics
76!     (mass (kg.m-2), thickness (m), etc.
77
78      do  l=1,nlay
79        do ig=1, ngrid
80          masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
81          epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
82        end do
83      end do
84
85!======================================================================
86! Calculate the transport due to sedimentation for each tracer
87! [This has been rearranged by L. Kerber to allow the sedimentation
88!  of general tracers.]
89 
90      zqi(1:ngrid,1:nlay,1:nq) = 0.
91      do iq=1,nq
92       if((radius(iq).gt.1.e-9).and.(iq.ne.igcm_co2_ice)) then   
93!         (no sedim for gases; co2_ice sedim is done elsewhere)     
94
95! store locally updated tracers
96
97          do l=1,nlay
98            do ig=1, ngrid
99              zqi(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
100            enddo
101          enddo ! of do l=1,nlay
102       
103!======================================================================
104! Sedimentation
105!======================================================================
106! Water         
107          if (water.and.(iq.eq.igcm_h2o_ice)) then
108             call newsedim(ngrid,nlay,ngrid*nlay,ptimestep,
109     &            pplev,masse,epaisseur,pt,rfall,rho_q(iq),zqi,wq)
110
111! General Case
112          else
113             call newsedim(ngrid,nlay,1,ptimestep,
114     &            pplev,masse,epaisseur,pt,radius(iq),rho_q(iq),
115     &            zqi,wq)
116          endif
117
118!=======================================================================
119!     Calculate the tendencies
120!======================================================================
121
122          do ig=1,ngrid
123!     Ehouarn: with new way of tracking tracers by name, this is simply
124            pdqs_sed(ig,iq) = wq(ig,1)/ptimestep
125          end do
126
127          DO l = 1, nlay
128            DO ig=1,ngrid
129              pdqsed(ig,l,iq)=(zqi(ig,l,iq)-
130     &        (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
131            ENDDO
132          ENDDO
133       endif ! of no gases no co2_ice
134      enddo ! of do iq=1,nq
135      return
136      end
Note: See TracBrowser for help on using the repository browser.