source: trunk/LMDZ.PLUTO/libf/phypluto/mp2m_calmufi.F90 @ 3955

Last change on this file since 3955 was 3683, checked in by debatzbr, 10 months ago

Solves random variations of the SW heating at the model top
BBT

File size: 8.2 KB
Line 
1MODULE mp2m_calmufi
2    use tracer_h
3    use callkeys_mod, only : call_haze_prod_pCH4, haze_rho
4   
5    ! Microphysical model MP2M
6    use mp2m_intgcm
7    use mp2m_diagnostics
8   
9
10    implicit none
11   
12    !============================================================================
13    !
14    !     Purpose
15    !     -------
16    !     Interface subroutine to YAMMS model for LMD PCM.
17    !
18    !     The subroutine computes the microphysics processes for a single vertical column.
19    !       - All input vectors are assumed to be defined from GROUND to TOP of the atmosphere.
20    !       - All output vectors are defined from GROUND to TOP of the atmosphere.
21    !       - Only tendencies are returned.
22    !
23    !     @important
24    !     The method assumes global initialization of YAMMS model (and extras) has been already
25    !     done elsewhere.
26    !
27    !     @warning
28    !     Microphysical tracers from physics must be in X/kg_of_air and convert into X/m2 for microphysics.
29    !
30    !     @warning
31    !     We suppose a given order of tracers (1. mu_m0as, 2. mu_m3as, 3. mu_m0af, 4. mu_m3af)!
32    !
33    !     Authors
34    !     -------
35    !     B. de Batz de Trenquelléon (11/2024)
36    !
37    !============================================================================
38
39    CONTAINS
40   
41    SUBROUTINE calmufi(dt, plev, zlev, play, zlay, g3d, temp, pq, zdqfi, zdqmufi_prod, zdqmufi)
42        !! Interface subroutine to YAMMS model for LMD PCM.
43        !!
44        !! The subroutine computes the microphysics processes for a single vertical column.
45        !!
46        REAL(kind=8), INTENT(IN) :: dt  !! Physics timestep (s).
47       
48        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: plev            ! Pressure levels (Pa).
49        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlev            ! Altitude levels (m).
50        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: play            ! Pressure layers (Pa).
51        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlay            ! Altitude at the center of each layer (m).
52        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: g3d             ! Latitude-Altitude depending gravitational acceleration (m.s-2).
53        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: temp            ! Temperature at the center of each layer (K).
54       
55        REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: pq           ! Tracers (X.kg-1).
56        REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: zdqfi        ! Tendency from former processes for tracers (X.kg-1.s-1).
57        REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: zdqmufi_prod ! Aerosols production tendency (kg/kg_of_air/s).
58        REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdqmufi      ! Microphysical tendency for tracers (X.m-2 --> X.kg-1.s-1).
59       
60        ! Local tracers:
61        !~~~~~~~~~~~~~~~
62        REAL(kind=8), DIMENSION(:,:,:), ALLOCATABLE :: zq    ! Local tracers updated from former processes (X.kg-1).
63       
64        REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m0as      ! 0th order moment of the spherical mode (m-2).
65        REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m3as      ! 3rd order moment of the spherical mode (m3.m-2).
66        REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m0af      ! 0th order moment of the fractal mode (m-2).
67        REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m3af      ! 3rd order moment of the fractal mode (m3.m-2).
68
69        REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m3as_prod ! Production of 3rd order moment of the spherical mode (m3.m-2).
70       
71        REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: dm0as   ! Tendency of the 0th order moment of the spherical mode distribution (m-2).
72        REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: dm3as   ! Tendency of the 3rd order moment of the spherical mode distribution (m3.m-2).
73        REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: dm0af   ! Tendency of the 0th order moment of the fractal mode distribution (m-2).
74        REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: dm3af   ! Tendency of the 3rd order moment of the fractal mode distribution (m3.m-2).
75       
76        ! Local variables:
77        !~~~~~~~~~~~~~~~~~
78        REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: int2ext ! Conversion intensive to extensive (kg.m-2).
79        REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: tmp
80        TYPE(error) :: err
81       
82        INTEGER           :: ilon,iq
83        INTEGER           :: nq,nlon,nlay
84        CHARACTER(len=10) :: tname
85       
86        ! Read size of arrays:
87        !~~~~~~~~~~~~~~~~~~~~~
88        nq    = size(pq,DIM=3)
89        nlon  = size(play,DIM=1)
90        nlay  = size(play,DIM=2)
91       
92        ! Allocate arrays:
93        !~~~~~~~~~~~~~~~~~
94        ALLOCATE(zq(nlon,nlay,nq))
95       
96        ALLOCATE(m0as(nlay))
97        ALLOCATE(m3as(nlay))
98        ALLOCATE(m0af(nlay))
99        ALLOCATE(m3af(nlay))
100
101        ALLOCATE(m3as_prod(nlay))
102       
103        ALLOCATE(dm0as(nlay))
104        ALLOCATE(dm3as(nlay))
105        ALLOCATE(dm0af(nlay))
106        ALLOCATE(dm3af(nlay))
107
108        ALLOCATE(int2ext(nlon,nlay)) 
109       
110        !------------------
111        ! 1. Initialization
112        !------------------
113
114        ! Initialization of zdqmufi here since intent=out and no action performed on every tracers
115        zdqmufi(:,:,:) = 0.D0
116       
117        ! Initialize tracers updated with former processes from physics
118        zq(:,:,:) = pq(:,:,:) + zdqfi(:,:,:)*dt
119       
120        ! Loop on horizontal grid points
121        DO ilon = 1, nlon
122
123            ! Convert tracers to extensive [X.kg-1 --> X.m-2]
124            int2ext(ilon,:) = (plev(ilon,1:nlay)-plev(ilon,2:nlay+1)) / g3d(ilon,1:nlay)
125           
126            m0as(:) = zq(ilon,:,micro_indx(1)) * int2ext(ilon,:)
127            m3as(:) = zq(ilon,:,micro_indx(2)) * int2ext(ilon,:)
128           
129            m0af(:) = zq(ilon,:,micro_indx(3)) * int2ext(ilon,:)
130            m3af(:) = zq(ilon,:,micro_indx(4)) * int2ext(ilon,:)
131
132            ! Production of haze in the atmosphere by photolysis of CH4
133            if (call_haze_prod_pCH4) then
134                do iq = 1, nq
135                    tname = noms(iq)
136                    if (tname(1:4).eq."haze") then
137                        m3as_prod(:) = zdqmufi_prod(ilon,:,iq) * (int2ext(ilon,:) / haze_rho) * dt
138                    endif
139                enddo
140            else
141                m3as_prod(:) = 0.D0
142            endif
143           
144            ! Hackin the pressure level
145            tmp = plev(ilon,:)
146            if (tmp(nlay+1) == 0.0) then
147            tmp(nlay+1) = 2*tmp(nlay) - tmp(nlay-1)
148            endif
149       
150            ! Initialize YAMMS atmospheric column
151            err = mm_column_init(tmp,zlev(ilon,:),play(ilon,:),zlay(ilon,:),temp(ilon,:)) ; IF (err /= 0) call abort_program(err)
152            ! Initialize YAMMS aerosols moments column
153            err = mm_aerosols_init(m0as,m3as,m0af,m3af) ; IF (err /= 0) call abort_program(err)
154           
155            ! Initializes tendencies
156            dm0as(:) = 0._mm_wp ; dm3as(:) = 0._mm_wp ; dm0af(:) = 0._mm_wp ; dm3af(:) = 0._mm_wp
157           
158        !----------------------------
159        ! 2. Call microphysical model
160        !----------------------------
161           
162            ! Call microphysics
163            IF (.NOT.mm_muphys(m3as_prod,dm0as,dm3as,dm0af,dm3af)) THEN
164                call abort_program(error("mm_muphys aborted -> initialization not done !",-1))
165            ENDIF
166           
167            ! Save diagnostics
168            call mm_diagnostics(dt,mp2m_aer_s_prec(ilon),mp2m_aer_f_prec(ilon),  &
169                                mp2m_aer_s_w(ilon,:),mp2m_aer_f_w(ilon,:),       &
170                                mp2m_aer_s_flux(ilon,:),mp2m_aer_f_flux(ilon,:), &
171                                mp2m_rc_sph(ilon,:),mp2m_rc_fra(ilon,:))
172       
173            ! Convert tracers back to intensives
174            zdqmufi(ilon,:,micro_indx(1)) = dm0as(:) / int2ext(ilon,:)
175            zdqmufi(ilon,:,micro_indx(2)) = dm3as(:) / int2ext(ilon,:)
176            zdqmufi(ilon,:,micro_indx(3)) = dm0af(:) / int2ext(ilon,:)
177            zdqmufi(ilon,:,micro_indx(4)) = dm3af(:) / int2ext(ilon,:)
178           
179        END DO ! End loop on ilon
180
181        ! YAMMS gives a tendency which is integrated for all the timestep but in the GCM
182        ! we want to have routines spitting tendencies in s-1 -> let's divide !
183        zdqmufi(:,:,:) = zdqmufi(:,:,:) / dt
184       
185    END SUBROUTINE calmufi
186
187END MODULE mp2m_calmufi
Note: See TracBrowser for help on using the repository browser.