source: trunk/LMDZ.GENERIC/libf/phystd/condensation_generic_mod.F90 @ 2708

Last change on this file since 2708 was 2708, checked in by aslmd, 2 years ago

error messages when setting a traceur.def with wrong names , or with vap/ice tracer missing for the same specie , or if the specie is not in table_tracers_condensable

File size: 8.5 KB
Line 
1module condensation_generic_mod
2    implicit none
3   
4contains
5   
6    subroutine condensation_generic(ngrid,nlayer,nq,ptimestep, pplev, pplay,    &
7                pt, pq, pdt, pdq, pdtlsc, pdqvaplsc, pdqliqlsc)
8        use ioipsl_getin_p_mod, only: getin_p !-> to get the metallicity
9        use generic_cloud_common_h, only : RLVTT, cpp, &
10        Psat_generic,Lcpdqsat_generic,specie_parameters,specie_parameters_table
11        USE tracer_h
12        IMPLICIT none
13
14!=======================================================================
15!     
16!     Purpose
17!     -------
18!     Calculates large-scale condensation of generic tracer "tname".
19!     By convention, tname ends with the suffix "_vap", as it represents the
20!     gas phase of the generic tracer
21!     
22!     Authors
23!     -------
24!     Adapted from largescale.F90 by Lucas Teinturier (2022)
25!     largescale.F90 adapted from the LMDTERRE code by R. Wordsworth (2009)
26!     Original author Z. X. Li (1993)
27!     
28!=========================================================================
29
30        INTEGER, intent(in) :: ngrid,nlayer,nq
31
32!       Arguments
33        REAL, intent(in) :: ptimestep             ! intervalle du temps (s)
34        REAL, intent(in) :: pplev(ngrid,nlayer+1) ! pression a inter-couche
35        REAL, intent(in) :: pplay(ngrid,nlayer)   ! pression au milieu de couche
36        REAL, intent(in) :: pt(ngrid,nlayer)      ! temperature (K)
37        REAL, intent(in) :: pq(ngrid,nlayer,nq)   ! tracer mixing ratio (kg/kg)
38        REAL, intent(in) :: pdt(ngrid,nlayer)     ! physical temperature tendency (K/s)
39        REAL, intent(in) :: pdq(ngrid,nlayer,nq)  ! physical tracer tendency (K/s)
40        ! CHARACTER(*), intent(in) :: tname_vap     ! name of the tracer we consider. BY convention, it ends with _vap !!!
41        REAL, intent(out) :: pdtlsc(ngrid,nlayer)  ! incrementation de la temperature (K)
42        REAL, intent(out) :: pdqvaplsc(ngrid,nlayer,nq) ! incrementation de la vapeur du traceur
43        REAL, intent(out) :: pdqliqlsc(ngrid,nlayer,nq) ! incrementation du traceur liquide
44
45!       Options :
46        real, save :: metallicity !metallicity of planet
47!$OMP THREADPRIVATE(metallicity)
48
49!       Local variables
50        INTEGER i, k , nn, iq
51        INTEGER,PARAMETER :: nitermax=5000
52        DOUBLE PRECISION,PARAMETER :: alpha=.1,qthreshold=1.d-8
53        ! JL13: if "careful, T<Tmin in psat water" appears often, you may want to stabilise the model by
54        !                   decreasing alpha and increasing nitermax accordingly
55        DOUBLE PRECISION zq(ngrid)
56        DOUBLE PRECISION zcond(ngrid),zcond_iter
57        DOUBLE PRECISION zqs(ngrid)
58        real zt(ngrid),local_p,psat_tmp,dlnpsat_tmp,Lcp,zqs_temp,zdqs
59        integer igcm_generic_vap, igcm_generic_ice! index of the vap and ice of generic_tracer
60        ! CHARACTER(len=*) :: tname_ice
61        ! evaporation calculations
62        REAL dqevap(ngrid,nlayer),dtevap(ngrid,nlayer)     
63        REAL qevap(ngrid,nlayer,nq)
64        REAL tevap(ngrid,nlayer)
65
66        DOUBLE PRECISION zx_q(ngrid)
67        LOGICAL,SAVE :: firstcall=.true.
68!$OMP THREADPRIVATE(firstcall)
69        IF (firstcall) THEN
70                write(*,*) "value for metallicity? "
71                metallicity=0.0 ! default value
72                call getin_p("metallicity",metallicity)
73                write(*,*) " metallicity = ",metallicity
74                firstcall = .false.
75        ENDIF
76!       Initialisation of outputs and local variables
77        pdtlsc(1:ngrid,1:nlayer)  = 0.0
78        pdqvaplsc(1:ngrid,1:nlayer,1:nq)  = 0.0
79        pdqliqlsc(1:ngrid,1:nlayer,1:nq) = 0.0
80        dqevap(1:ngrid,1:nlayer)=0.0
81        dtevap(1:ngrid,1:nlayer)=0.0
82        qevap(1:ngrid,1:nlayer,1:nq)=0.0
83        tevap(1:ngrid,1:nlayer)=0.0
84        ! Let's loop on tracers
85        do iq=1,nq
86                if((is_condensable(iq)==1) .and. (index(noms(iq),"vap") .ne. 0)) then
87                        write(*,*) "There is a specie which is condensable, for generic condensation : ", noms(iq)
88
89!                       Let's get the index of our tracers (we look for igcm _generic_vap and igcm_generic_ice)
90                        ! tname_ice = trim(noms(iq)(1:len(tname_ice)-3))//"ice"
91                        ! print*,trim(adjustl(trim(noms(iq))(9:len(trim(noms(iq)))-4))) !testing here, should go away
92                        ! stop
93
94                        igcm_generic_vap=iq
95                        igcm_generic_ice = -1
96
97                        ! We look for the corresponding ice traceur (before or after in the list of traceurs, maybe could be generalised to the whole list)
98                        if (iq .ne. nq) then
99                                if ((noms(iq)(1:len(trim(noms(iq)))-4) .eq. noms(iq+1)(1:len(trim(noms(iq+1)))-4)) .and. (index(noms(iq+1),"ice") .ne. 0)) then
100                                        igcm_generic_ice = iq+1
101                                end if
102                        end if
103                        if ((iq .gt. 1)) then
104                                if ((noms(iq)(1:len(trim(noms(iq)))-4) .eq. noms(iq-1)(1:len(trim(noms(iq-1)))-4)) .and. (index(noms(iq-1),"ice") .ne. 0)) then
105                                        igcm_generic_ice = iq-1
106                                end if
107                        end if
108                        if (igcm_generic_ice .eq. -1) then
109                                write(*,*) "ERROR : You set a vap traceur but you forgot to set the corresponding ice traceur"
110                        endif
111
112                        ! Need to call specie_parameters of the considered specie
113                        !call specie_parameters(noms(iq)(9:len(trim(noms(iq)))-4))
114                        write(*,*) 'looking specie parameters for : ',noms(iq)(1:len(trim(noms(iq)))-4)
115                        call specie_parameters_table(noms(iq)(1:len(trim(noms(iq)))-4))
116                       
117                        Lcp=RLVTT/cpp ! need to be init here
118
119                        !  Vertical loop (from top to bottom)
120                        DO k = nlayer, 1, -1
121                                zt(1:ngrid)=pt(1:ngrid,k)+pdt(1:ngrid,k)*ptimestep
122
123                                ! Computes Psat and the partial condensation
124                                DO i = 1, ngrid
125                                        local_p=pplay(i,k)
126                                        if(zt(i).le.15.) then
127                                                print*,'in lsc',i,k,zt(i)
128                                        !           zt(i)=15.   ! check too low temperatures
129                                        endif
130                                        zx_q(i) = pq(i,k,igcm_generic_vap)+pdq(i,k,igcm_generic_vap)*ptimestep
131                                        ! iterative process to stabilize the scheme when large water amounts JL12
132                                        zcond(i) = 0.0d0
133                                        Do nn=1,nitermax 
134                                                call Psat_generic(zt(i),local_p,metallicity,psat_tmp,zqs_temp)
135                                                zqs(i)=zqs_temp
136                                                call Lcpdqsat_generic(zt(i),local_p,psat_tmp,zqs_temp,zdqs,dlnpsat_tmp)
137                                                zcond_iter = alpha*(zx_q(i)-zqs(i))/(1.d0+zdqs)
138                                                !zcond can be negative here
139                                                zx_q(i) = zx_q(i) - zcond_iter
140                                                zcond(i) = zcond(i) + zcond_iter
141                                                zt(i) = zt(i) + zcond_iter*Lcp
142                                                if (ABS(zcond_iter/alpha/zqs(i)).lt.qthreshold) exit
143                                                if (nn.eq.nitermax) print*,'itermax in largescale'
144                                        End do ! niter
145                                        zcond(i)=MAX(zcond(i),-(pq(i,k,igcm_generic_ice)+pdq(i,k,igcm_generic_ice)*ptimestep))
146                                        zcond(i) = zcond(i)/ptimestep
147                                ENDDO ! i=1,ngrid
148
149                                !Tendances de t et q
150                                pdqvaplsc(1:ngrid,k,igcm_generic_vap)  = - zcond(1:ngrid)
151                                pdqliqlsc(1:ngrid,k,igcm_generic_ice) = - pdqvaplsc(1:ngrid,k,igcm_generic_vap)
152                                pdtlsc(1:ngrid,k)  = pdtlsc(1:ngrid,k) + pdqliqlsc(1:ngrid,k,igcm_generic_ice)*Lcp
153
154                        Enddo ! k= nlayer, 1, -1
155                endif !(is_condensable(iq)==1) .and. (index(noms(iq),"vap") .ne. 0))
156        enddo ! iq=1,nq
157
158    end subroutine condensation_generic
159end module condensation_generic_mod
Note: See TracBrowser for help on using the repository browser.