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

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

Print tendencies in diagfi.nc and small correction

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