| 1 | |
|---|
| 2 | |
|---|
| 3 | SUBROUTINE calmufi(dt, plev, zlev, play, zlay, g3d, temp, pq, zdqfi, zdq) |
|---|
| 4 | |
|---|
| 5 | ! Old compilers do not handle modern Fortran in microphysics |
|---|
| 6 | #ifndef OLD_COMPILO |
|---|
| 7 | |
|---|
| 8 | !! Interface subroutine to YAMMS model for Titan LMDZ GCM. |
|---|
| 9 | !! |
|---|
| 10 | !! The subroutine computes the microphysics processes for a single vertical column. |
|---|
| 11 | !! |
|---|
| 12 | !! - All input vectors are assumed to be defined from GROUND to TOP of the atmosphere. |
|---|
| 13 | !! - All output vectors are defined from GROUND to TOP of the atmosphere. |
|---|
| 14 | !! - Only tendencies are returned. |
|---|
| 15 | !! |
|---|
| 16 | !! @important |
|---|
| 17 | !! The method assumes global initialization of YAMMS model (and extras) has been already |
|---|
| 18 | !! done elsewhere. |
|---|
| 19 | !! |
|---|
| 20 | !! Authors : J.Burgalat, J.Vatant d'Ollone - 2017 |
|---|
| 21 | !! Modified : B. de Batz de Trenquelléon (2023) |
|---|
| 22 | !! |
|---|
| 23 | USE MMP_GCM |
|---|
| 24 | USE tracer_h |
|---|
| 25 | USE callkeys_mod, only : callclouds |
|---|
| 26 | USE muphy_diag |
|---|
| 27 | IMPLICIT NONE |
|---|
| 28 | |
|---|
| 29 | REAL(kind=8), INTENT(IN) :: dt !! Physics timestep (s). |
|---|
| 30 | |
|---|
| 31 | REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: plev !! Pressure levels (Pa). |
|---|
| 32 | REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlev !! Altitude levels (m). |
|---|
| 33 | REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: play !! Pressure layers (Pa). |
|---|
| 34 | REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlay !! Altitude at the center of each layer (m). |
|---|
| 35 | REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: g3d !! Latitude-Altitude depending gravitational acceleration (m.s-2). |
|---|
| 36 | REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: temp !! Temperature at the center of each layer (K). |
|---|
| 37 | |
|---|
| 38 | REAL(kind=8), DIMENSION(:,:,:), INTENT(IN) :: pq !! Tracers (\(X.kg^{-1}}\)). |
|---|
| 39 | REAL(kind=8), DIMENSION(:,:,:), INTENT(IN) :: zdqfi !! Tendency from former processes for tracers (\(X.kg^{-1}}\)). |
|---|
| 40 | REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdq !! Microphysical tendency for tracers (\(X.kg^{-1}}\)). |
|---|
| 41 | |
|---|
| 42 | REAL(kind=8), DIMENSION(:,:,:), ALLOCATABLE :: zq !! Local tracers updated from former processes (\(X.kg^{-1}}\)). |
|---|
| 43 | |
|---|
| 44 | REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m0as !! 0th order moment of the spherical mode (\(m^{-2}\)). |
|---|
| 45 | REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m3as !! 3rd order moment of the spherical mode (\(m^{3}.m^{-2}\)). |
|---|
| 46 | REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m0af !! 0th order moment of the fractal mode (\(m^{-2}\)). |
|---|
| 47 | REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m3af !! 3rd order moment of the fractal mode (\(m^{3}.m^{-2}\)). |
|---|
| 48 | |
|---|
| 49 | REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m0n !! 0th order moment of the CCN distribution (\(m^{-2}\)). |
|---|
| 50 | REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m3n !! 3rd order moment of the CCN distribution (\(m^{3}.m^{-2}\)). |
|---|
| 51 | REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: m3i !! 3rd order moments of the ice components (\(m^{3}.m^{-2}\)). |
|---|
| 52 | REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: gazs !! Condensible species gazs molar fraction (\(mol.mol^{-1}\)). |
|---|
| 53 | |
|---|
| 54 | REAL(kind=8), DIMENSION(:), ALLOCATABLE :: dm0as !! Tendency of the 0th order moment of the spherical mode distribution (\(m^{-2}\)). |
|---|
| 55 | REAL(kind=8), DIMENSION(:), ALLOCATABLE :: dm3as !! Tendency of the 3rd order moment of the spherical mode distribution (\(m^{3}.m^{-2}\)). |
|---|
| 56 | REAL(kind=8), DIMENSION(:), ALLOCATABLE :: dm0af !! Tendency of the 0th order moment of the fractal mode distribution (\(m^{-2}\)). |
|---|
| 57 | REAL(kind=8), DIMENSION(:), ALLOCATABLE :: dm3af !! Tendency of the 3rd order moment of the fractal mode distribution (\(m^{3}.m^{-2}\)). |
|---|
| 58 | REAL(kind=8), DIMENSION(:), ALLOCATABLE :: dm0n !! Tendency of the 0th order moment of the _CCN_ distribution (\(m^{-2}\)). |
|---|
| 59 | REAL(kind=8), DIMENSION(:), ALLOCATABLE :: dm3n !! Tendency of the 3rd order moment of the _CCN_ distribution (\(m^{3}.m^{-2}\)). |
|---|
| 60 | REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: dm3i !! Tendencies of the 3rd order moments of each ice components (\(m^{3}.m^{-2}\)). |
|---|
| 61 | REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: dgazs !! Tendencies of each condensible gaz species !(\(mol.mol^{-1}\)). |
|---|
| 62 | |
|---|
| 63 | REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: int2ext !! (\(m^{-2}\)). |
|---|
| 64 | REAL(kind=8), DIMENSION(:), ALLOCATABLE :: tmp |
|---|
| 65 | TYPE(error) :: err |
|---|
| 66 | |
|---|
| 67 | INTEGER :: ilon, i,nices |
|---|
| 68 | INTEGER :: nq,nlon,nlay |
|---|
| 69 | |
|---|
| 70 | ! Read size of arrays |
|---|
| 71 | nq = size(pq,DIM=3) |
|---|
| 72 | nlon = size(play,DIM=1) |
|---|
| 73 | nlay = size(play,DIM=2) |
|---|
| 74 | nices = size(ices_indx) |
|---|
| 75 | ! Conversion intensive to extensive |
|---|
| 76 | ALLOCATE( int2ext(nlon,nlay) ) |
|---|
| 77 | |
|---|
| 78 | ! Allocate arrays |
|---|
| 79 | ALLOCATE( zq(nlon,nlay,nq) ) |
|---|
| 80 | |
|---|
| 81 | ALLOCATE( m0as(nlay) ) |
|---|
| 82 | ALLOCATE( m3as(nlay) ) |
|---|
| 83 | ALLOCATE( m0af(nlay) ) |
|---|
| 84 | ALLOCATE( m3af(nlay) ) |
|---|
| 85 | ALLOCATE( m0n(nlay) ) |
|---|
| 86 | ALLOCATE( m3n(nlay) ) |
|---|
| 87 | ALLOCATE( m3i(nlay,nices) ) |
|---|
| 88 | ALLOCATE( gazs(nlay,nices) ) |
|---|
| 89 | |
|---|
| 90 | ALLOCATE( dm0as(nlay) ) |
|---|
| 91 | ALLOCATE( dm3as(nlay) ) |
|---|
| 92 | ALLOCATE( dm0af(nlay) ) |
|---|
| 93 | ALLOCATE( dm3af(nlay) ) |
|---|
| 94 | ALLOCATE( dm0n(nlay) ) |
|---|
| 95 | ALLOCATE( dm3n(nlay) ) |
|---|
| 96 | ALLOCATE( dm3i(nlay,nices) ) |
|---|
| 97 | ALLOCATE( dgazs(nlay,nices) ) |
|---|
| 98 | |
|---|
| 99 | ! Initialization of zdq here since intent=out and no action performed on every tracers |
|---|
| 100 | zdq(:,:,:) = 0.D0 |
|---|
| 101 | |
|---|
| 102 | ! Initialize tracers updated with former processes from physics |
|---|
| 103 | zq(:,:,:) = pq(:,:,:) + zdqfi(:,:,:)*dt |
|---|
| 104 | |
|---|
| 105 | ! Loop on horizontal grid points |
|---|
| 106 | DO ilon = 1, nlon |
|---|
| 107 | ! Convert tracers to extensive ( except for gazs where we work with molar mass ratio ) |
|---|
| 108 | ! We suppose a given order of tracers ! |
|---|
| 109 | int2ext(ilon,:) = ( plev(ilon,1:nlay)-plev(ilon,2:nlay+1) ) / g3d(ilon,1:nlay) |
|---|
| 110 | |
|---|
| 111 | ! Check because of the threshold of small tracers values in the dynamics |
|---|
| 112 | ! WARNING : With this patch it enables to handles the small values required |
|---|
| 113 | ! by YAMMS, but still it might still leads to some unphysical values of |
|---|
| 114 | ! radii inside YAMMS, harmless for now, but who might be a problem the day |
|---|
| 115 | ! you'll want to compute optics from the radii. |
|---|
| 116 | WHERE (pq(ilon,:,1) > 2.D-200 .AND. pq(ilon,:,2) > 2.D-200) ! Test on both moments to avoid divergences if one hit the threshold but not the other |
|---|
| 117 | m0as(:) = zq(ilon,:,1) * int2ext(ilon,:) ! It can still be a pb if both m0 and m3 has been set to epsilon at the beginning of dynamics |
|---|
| 118 | m3as(:) = zq(ilon,:,2) * int2ext(ilon,:) ! then mixed, even though both are above the threshold, their ratio can be nonsense |
|---|
| 119 | ELSEWHERE |
|---|
| 120 | m0as(:) = 0.D0 |
|---|
| 121 | m3as(:) = 0.D0 |
|---|
| 122 | ENDWHERE |
|---|
| 123 | WHERE (pq(ilon,:,3) > 2.D-200 .AND. pq(ilon,:,4) > 2.D-200) |
|---|
| 124 | m0af(:) = zq(ilon,:,3) * int2ext(ilon,:) |
|---|
| 125 | m3af(:) = zq(ilon,:,4) * int2ext(ilon,:) |
|---|
| 126 | ELSEWHERE |
|---|
| 127 | m0af(:) = 0.D0 |
|---|
| 128 | m3af(:) = 0.D0 |
|---|
| 129 | ENDWHERE |
|---|
| 130 | |
|---|
| 131 | if (callclouds) then ! if call clouds |
|---|
| 132 | m0n(:) = zq(ilon,:,5) * int2ext(ilon,:) |
|---|
| 133 | m3n(:) = zq(ilon,:,6) * int2ext(ilon,:) |
|---|
| 134 | do i=1,nices |
|---|
| 135 | m3i(:,i) = zq(ilon,:,ices_indx(i)) * int2ext(ilon,:) |
|---|
| 136 | gazs(:,i) = zq(ilon,:,gazs_indx(i)) / rat_mmol(gazs_indx(i)) ! For gazs we work on the full tracer array !! |
|---|
| 137 | ! We use the molar mass ratio from GCM in case there is discrepancy with the mm one |
|---|
| 138 | enddo |
|---|
| 139 | endif |
|---|
| 140 | |
|---|
| 141 | ! Hackin the pressure level |
|---|
| 142 | tmp = plev(ilon,:) |
|---|
| 143 | if (tmp(nlay+1) == 0.0) then |
|---|
| 144 | tmp(nlay+1) = 2*tmp(nlay) - tmp(nlay-1) |
|---|
| 145 | endif |
|---|
| 146 | |
|---|
| 147 | ! Initialize YAMMS atmospheric column |
|---|
| 148 | err = mm_column_init(tmp,zlev(ilon,:),play(ilon,:),zlay(ilon,:),temp(ilon,:)) ; IF (err /= 0) call abort_program(err) |
|---|
| 149 | ! Initialize YAMMS aerosols moments column |
|---|
| 150 | err = mm_aerosols_init(m0as,m3as,m0af,m3af) ; IF (err /= 0) call abort_program(err) |
|---|
| 151 | IF (callclouds) THEN ! call clouds |
|---|
| 152 | err = mm_clouds_init(m0n,m3n,m3i,gazs) ; IF (err /= 0) call abort_program(err) |
|---|
| 153 | ENDIF |
|---|
| 154 | |
|---|
| 155 | ! Check on size (???) |
|---|
| 156 | |
|---|
| 157 | ! initializes tendencies: |
|---|
| 158 | dm0as(:) = 0._mm_wp ; dm3as(:) = 0._mm_wp ; dm0af(:) = 0._mm_wp ; dm3af(:) = 0._mm_wp |
|---|
| 159 | dm0n(:) = 0._mm_wp ; dm3n(:) = 0._mm_wp ; dm3i(:,:) = 0._mm_wp ; dgazs(:,:) = 0._mm_wp |
|---|
| 160 | |
|---|
| 161 | ! call microphysics |
|---|
| 162 | IF (callclouds) THEN ! call clouds |
|---|
| 163 | IF(.NOT.mm_muphys(dm0as,dm3as,dm0af,dm3af,dm0n,dm3n,dm3i,dgazs)) & |
|---|
| 164 | call abort_program(error("mm_muphys aborted -> initialization not done !",-1)) |
|---|
| 165 | ELSE |
|---|
| 166 | IF (.NOT.mm_muphys(dm0as,dm3as,dm0af,dm3af)) & |
|---|
| 167 | call abort_program(error("mm_muphys aborted -> initialization not done !",-1)) |
|---|
| 168 | ENDIF |
|---|
| 169 | ! save diags (if no clouds, relevant arrays will be set to 0 !) |
|---|
| 170 | call mm_diagnostics(dt,mmd_aer_prec(ilon),mmd_aer_s_w(ilon,:),mmd_aer_f_w(ilon,:),mmd_aer_s_flux(ilon,:),mmd_aer_f_flux(ilon,:), & |
|---|
| 171 | mmd_ccn_prec(ilon),mmd_ccn_w(ilon,:),mmd_ccn_flux(ilon,:),mmd_ice_prec(ilon,:), & |
|---|
| 172 | mmd_ice_fluxes(ilon,:,:),mmd_gazs_sat(ilon,:,:)) |
|---|
| 173 | call mm_get_radii(mmd_rc_sph(ilon,:),mmd_rc_fra(ilon,:),mmd_rc_cld(ilon,:)) |
|---|
| 174 | |
|---|
| 175 | ! Convert tracers back to intensives ( except for gazs where we work with molar mass ratio ) |
|---|
| 176 | ! We suppose a given order of tracers ! |
|---|
| 177 | |
|---|
| 178 | zdq(ilon,:,1) = dm0as(:) / int2ext(ilon,:) |
|---|
| 179 | zdq(ilon,:,2) = dm3as(:) / int2ext(ilon,:) |
|---|
| 180 | zdq(ilon,:,3) = dm0af(:) / int2ext(ilon,:) |
|---|
| 181 | zdq(ilon,:,4) = dm3af(:) / int2ext(ilon,:) |
|---|
| 182 | |
|---|
| 183 | if (callclouds) then ! if call clouds |
|---|
| 184 | zdq(ilon,:,5) = dm0n(:) / int2ext(ilon,:) |
|---|
| 185 | zdq(ilon,:,6) = dm3n(:) / int2ext(ilon,:) |
|---|
| 186 | do i=1,nices |
|---|
| 187 | zdq(ilon,:,ices_indx(i)) = dm3i(:,i) / int2ext(ilon,:) |
|---|
| 188 | zdq(ilon,:,gazs_indx(i)) = dgazs(:,i) * rat_mmol(gazs_indx(i)) ! For gazs we work on the full tracer array !! |
|---|
| 189 | ! We use the molar mass ratio from GCM in case there is discrepancy with the mm one |
|---|
| 190 | enddo |
|---|
| 191 | endif |
|---|
| 192 | |
|---|
| 193 | END DO ! loop on ilon |
|---|
| 194 | |
|---|
| 195 | ! YAMMS gives a tendency which is integrated for all the timestep but in the GCM |
|---|
| 196 | ! we want to have routines spitting tendencies in s-1 -> let's divide ! |
|---|
| 197 | zdq(:,:,:) = zdq(:,:,:) / dt |
|---|
| 198 | |
|---|
| 199 | #endif |
|---|
| 200 | |
|---|
| 201 | END SUBROUTINE calmufi |
|---|
