Changeset 2242
- Timestamp:
- Feb 17, 2020, 5:44:12 PM (5 years ago)
- Location:
- trunk
- Files:
-
- 8 edited
-
LMDZ.COMMON/libf/dyn3dpar/addfi_p.F (modified) (2 diffs)
-
LMDZ.TITAN/README (modified) (1 diff)
-
LMDZ.TITAN/libf/muphytitan/mm_globals.f90 (modified) (1 diff)
-
LMDZ.TITAN/libf/phytitan/calmufi.F90 (modified) (4 diffs)
-
LMDZ.TITAN/libf/phytitan/datafile_mod.F90 (modified) (2 diffs)
-
LMDZ.TITAN/libf/phytitan/optci.F90 (modified) (7 diffs)
-
LMDZ.TITAN/libf/phytitan/optcv.F90 (modified) (9 diffs)
-
LMDZ.TITAN/libf/phytitan/physiq_mod.F90 (modified) (16 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/dyn3dpar/addfi_p.F
r1422 r2242 78 78 REAL,PARAMETER :: qtestw = 1.0e-15 79 79 REAL,PARAMETER :: qtestt = 1.0e-40 80 REAL,PARAMETER :: qtestt2 = 1.0D-200 80 81 81 82 REAL SSUM … … 192 193 c$OMP END DO 193 194 ENDDO 195 else if (planet_type=="titan") then 196 ! Titan : needs to be able to deal with very low values of tracers 197 ! values for microphysics -> threshold at 1D-200 - JVO 20 198 DO iq = 1, nqtot 199 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 200 DO k = 1,llm 201 DO j = ijb,ije 202 pq(j,k,iq)= pq(j,k,iq) + pdqfi(j,k,iq) * pdt 203 pq(j,k,iq)= AMAX1( pq(j,k,iq), qtestt2 ) 204 ENDDO 205 ENDDO 206 c$OMP END DO 207 ENDDO 194 208 else 195 209 ! general case, treat all tracers equally) -
trunk/LMDZ.TITAN/README
r2138 r2242 1563 1563 but would require to be able to have writediag in 5D) 1564 1564 EDIT (22/05/19) : To have correct calculations, output is now exp(-tau), you need to postproc it to have tau. 1565 1566 == 17/02/2020 == JVO (r2242) 1567 - Update the interface with YAMMS so it now correctly handles the small values of the moments, 1568 requiring dynamics to have a threshold quite low (set to 1D-200) and a sanity check in calmufi 1569 corresponding to this value. Thus it removes 'most' of the unphysical radius obtained in 1570 YAMMS. There are still some, but at least there is no more problem of model stability for the moments 1571 and the code can run. 1572 **** Still, take care the day you want to calculate opacities from the radii and not the moments. 1573 **** Although, note that there are some negative values, in the output files, but theses negatives are 1574 harmless, as they are only present in output files, dynamics reseting to epsilon after. 1575 - TBD ? Add a sanity check within YAMMS filtering high unphysical radii ? 1576 - Given theses pbs with the radii, also update the optics so that it computes the opacities in 1577 a simpler way, directly for M3, through look-up tables, M3 being a good proxy. -
trunk/LMDZ.TITAN/libf/muphytitan/mm_globals.f90
r1897 r2242 1066 1066 ! il faudrait peut etre revoir la gestion des zeros ici et la 1067 1067 ! remplacer par une valeur seuil des moments. 1068 ! 1069 !-> JVO 19 : Done. Zero threshold set at espilon from dynamics on the 1070 ! input moments in calmufi (safer than here). Might still be some unphysical 1071 ! values after the dynamics near the threshold. Could be a could idea to add 1072 ! a sanity check filtering too high values of radii. 1073 ! 1074 ! TBD : Add a sanity check for high radii ???? 1068 1075 WHERE(mm_m3aer_s > 0._mm_wp .AND. mm_m0aer_s > 0._mm_wp) 1069 1076 mm_rcs = mm_get_rcs(mm_m0aer_s,mm_m3aer_s) -
trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90
r2109 r2242 92 92 93 93 ! Initialization of zdq here since intent=out and no action performed on every tracers 94 zdq(:,:,:) = 0. 094 zdq(:,:,:) = 0.D0 95 95 96 96 ! Initialize tracers updated with former processes from physics … … 102 102 ! We suppose a given order of tracers ! 103 103 int2ext(ilon,:) = ( plev(ilon,1:nlay)-plev(ilon,2:nlay+1) ) / g3d(ilon,1:nlay) 104 m0as(:) = zq(ilon,:,1) * int2ext(ilon,:) 105 m3as(:) = zq(ilon,:,2) * int2ext(ilon,:) 106 m0af(:) = zq(ilon,:,3) * int2ext(ilon,:) 107 m3af(:) = zq(ilon,:,4) * int2ext(ilon,:) 104 105 ! Check because of the threshold of small tracers values in the dynamics 106 ! WARNING : With this patch it enables to handles the small values required 107 ! by YAMMS, but still it might still leads to some unphysical values of 108 ! radii inside YAMMS, harmless for now, but who might be a problem the day 109 ! you'll want to compute optics from the radii. 110 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 111 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 112 m3as(:) = zq(ilon,:,2) * int2ext(ilon,:) ! then mixed, even though both are above the threshold, their ratio can be nonsense 113 ELSEWHERE 114 m0as(:)=0.D0 115 m3as(:)=0.D0 116 ENDWHERE 117 WHERE (pq(ilon,:,3) > 2.D-200 .AND. pq(ilon,:,4) > 2.D-200) 118 m0af(:) = zq(ilon,:,3) * int2ext(ilon,:) 119 m3af(:) = zq(ilon,:,4) * int2ext(ilon,:) 120 ELSEWHERE 121 m0af(:)=0.D0 122 m3af(:)=0.D0 123 ENDWHERE 108 124 109 125 if (callclouds) then ! if call clouds … … 131 147 !dm0n(:) = 0._mm_wp ; dm3n(:) = 0._mm_wp ; dm3i(:,:) = 0._mm_wp ; dgazs(:,:) = 0._mm_wp 132 148 133 dm0as(:) = 0. 0 ; dm3as(:) = 0.0 ; dm0af(:) = 0.0 ; dm3af(:) = 0.0134 dm0n(:) = 0. 0 ; dm3n(:) = 0.0 ; dm3i(:,:) = 0.0 ; dgazs(:,:) = 0.0149 dm0as(:) = 0.D0 ; dm3as(:) = 0.D0 ; dm0af(:) = 0.D0 ; dm3af(:) = 0.D0 150 dm0n(:) = 0.D0 ; dm3n(:) = 0.D0 ; dm3i(:,:) = 0.D0 ; dgazs(:,:) = 0.D0 135 151 ! call microphysics 136 152 … … 165 181 enddo 166 182 endif 167 168 ! Sanity check ( way safer to be done here rather than within YAMMS ) 169 WHERE (zq+zdq < 0.0) ; zdq = -zq ; END WHERE 170 183 171 184 END DO ! loop on ilon 172 185 -
trunk/LMDZ.TITAN/libf/phytitan/datafile_mod.F90
r1822 r2242 15 15 16 16 ! Default directories for correlated-k data 17 ! Set in inifis_mod17 ! Read in inifis_mod, set in physiq_mod 18 18 character(LEN=100),save :: corrkdir = 'datagcm/corrk_data/50X50' 19 19 character(LEN=100),save :: banddir = 'datagcm/corrk_data/50X50/50x50' … … 25 25 !$OMP THREADPRIVATE(config_mufi) 26 26 27 ! Default file for coupled microphysics optical properties 28 ! Set in physiq_mod 29 character(LEN=100),save :: haze_opt_file ='datagcm/HAZE_OPT_23x23.DAT' 30 !$OMP THREADPRIVATE(haze_opt_file) 31 27 32 end module datafile_mod 28 33 !----------------------------------------------------------------------- -
trunk/LMDZ.TITAN/libf/phytitan/optci.F90
r2133 r2242 6 6 tgasref,pfgasref,wnoi,scalep,indi 7 7 use gases_h 8 use datafile_mod, only: haze_opt_file 8 9 use comcstfi_mod, only: r 9 10 use callkeys_mod, only: continuum,graybody,corrk_recombin, & 10 11 callclouds,callmufi,seashaze,uncoupl_optic_haze 11 12 use tracer_h, only : nmicro,nice 12 use MMP_OPTICS13 13 14 14 implicit none … … 83 83 ! variables for k in units m^-1 84 84 real*8 dz(L_LEVELS) 85 !real*8 rho !! see test below86 85 87 86 integer igas, jgas, ilay … … 89 88 integer interm 90 89 91 real*8 m0as,m3as,m0af,m3af 92 real*8 ext_s,sca_s,ssa_s,asf_s 93 real*8 ext_f,sca_f,ssa_f,asf_f 90 ! Variables for haze optics 91 character(len=200) file_path 92 logical file_ok 93 integer dumch 94 real*8 dumwvl 95 96 real*8 m3as,m3af 97 real*8 dtauaer_s,dtauaer_f 98 real*8,save :: rhoaer_s(L_NSPECTI),ssa_s(L_NSPECTI),asf_s(L_NSPECTI) 99 real*8,save :: rhoaer_f(L_NSPECTI),ssa_f(L_NSPECTI),asf_f(L_NSPECTI) 100 !$OMP THREADPRIVATE(rhoaer_s,rhoaer_f,ssa_s,ssa_f,asf_s,asf_f) 101 94 102 logical,save :: firstcall=.true. 95 !$OMP THREADPRIVATE(firstcall) 103 !$OMP THREADPRIVATE(firstcall) 104 96 105 97 106 !! AS: to save time in computing continuum (see bilinearbig) … … 101 110 ENDIF 102 111 103 ! Some initialisation be acause there's a pb with disr_haze at the limits (nw=1)112 ! Some initialisation because there's a pb with disr_haze at the limits (nw=1) 104 113 ! I should check this - For now we set vars to zero : better than nans - JVO 2017 105 106 dhaze_t(:,:) = 0. 107 ssa_t(:,:) = 0. 108 asf_t(:,:) = 0. 114 DHAZE_T(:,:) = 0.0 115 SSA_T(:,:) = 0.0 116 ASF_T(:,:) = 0.0 117 118 ! Load tabulated haze optical properties if needed. 119 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 120 IF (firstcall .AND. callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN 121 OPEN(12,file=TRIM(haze_opt_file),form='formatted') ! The file has been inquired in physiq_mod firstcall 122 READ(12,*) ! dummy header 123 DO NW=1,L_NSPECTI 124 READ(12,*) dumch, dumwvl, rhoaer_f(nw), ssa_f(nw), asf_f(nw), rhoaer_s(nw), ssa_s(nw), asf_s(nw) 125 ENDDO 126 CLOSE(12) 127 ENDIF 109 128 110 129 !======================================================================= … … 118 137 do K=2,L_LEVELS 119 138 120 ilay = k / 2 ! int. arithmetic => gives the gcm layer index139 ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed) 121 140 122 141 DPR(k) = PLEV(K)-PLEV(K-1) … … 137 156 do K=2,L_LEVELS 138 157 139 ilay = k / 2 ! int. arithmetic => gives the gcm layer index158 ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed) 140 159 141 ! Optical coupling of YAMMS is plugged but inactivated for now142 ! as long as the microphysics only isn't fully debugged -- JVO 01/18143 160 IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN 144 m0as = pqmo(ilay,1) 145 m3as = pqmo(ilay,2) 146 m0af = pqmo(ilay,3) 147 m3af = pqmo(ilay,4) 148 149 IF (.NOT.mmp_sph_optics_ir(m0as,m3as,nw,ext_s,sca_s,ssa_s,asf_s)) & 150 CALL abort_gcm("optcv", "Fatal error in mmp_sph_optics_ir", 12) 151 IF (.NOT.mmp_fra_optics_ir(m0af,m3af,nw,ext_f,sca_f,ssa_f,asf_f)) & 152 CALL abort_gcm("optcv", "Fatal error in mmp_fra_optics_ir", 12) 153 dhaze_T(k,nw) = ext_s+ext_f 154 SSA_T(k,nw) = (sca_s+sca_f)/dhaze_T(k,nw) 155 ASF_T(k,nw) = (asf_s*sca_s + asf_f*sca_f) /(sca_s+sca_f) 161 m3as = pqmo(ilay,2) / 2.0 162 m3af = pqmo(ilay,4) / 2.0 163 164 IF ( ilay .lt. 18 ) THEN 165 m3as = pqmo(18,2) / 2.0 *dz(k)/dz(18) 166 m3af = pqmo(18,4) / 2.0 *dz(k)/dz(18) 167 ENDIF 168 169 dtauaer_s = m3as*rhoaer_s(nw) 170 dtauaer_f = m3af*rhoaer_f(nw) 171 DHAZE_T(k,nw) = dtauaer_s + dtauaer_f 172 173 IF ( dtauaer_s + dtauaer_f .GT. 1.D-30 ) THEN 174 SSA_T(k,nw) = ( dtauaer_s*ssa_s(nw) + dtauaer_f*ssa_f(nw) ) / ( dtauaer_s+dtauaer_f ) 175 ASF_T(k,nw) = ( dtauaer_s*ssa_s(nw)*asf_s(nw) + dtauaer_f*ssa_f(nw)*asf_f(nw) ) & 176 / ( ssa_s(nw)*dtauaer_s + ssa_f(nw)*dtauaer_f ) 177 ELSE 178 DHAZE_T(k,nw) = 0.D0 179 SSA_T(k,nw) = 1.0 180 ASF_T(k,nw) = 1.0 181 ENDIF 182 156 183 IF (callclouds.and.firstcall) & 157 184 WRITE(*,*) 'WARNING: In optci, optical properties & … … 159 186 ELSE 160 187 ! Call fixed vertical haze profile of extinction - same for all columns 161 call disr_haze(dz(k),plev(k),wnoi(nw), dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))162 if (seashaze) dhaze_T(k,nw) = dhaze_T(k,nw)*seashazefact(k)188 call disr_haze(dz(k),plev(k),wnoi(nw),DHAZE_T(k,nw),SSA_T(k,nw),ASF_T(k,nw)) 189 if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k) 163 190 ENDIF 164 191 -
trunk/LMDZ.TITAN/libf/phytitan/optcv.F90
r2133 r2242 6 6 tgasref,pfgasref,wnov,scalep,indv 7 7 use gases_h 8 use datafile_mod, only: haze_opt_file 8 9 use comcstfi_mod, only: r 9 10 use callkeys_mod, only: continuum,graybody,callgasvis,corrk_recombin, & … … 99 100 integer interm 100 101 101 real*8 m0as,m3as,m0af,m3af 102 real*8 ext_s,sca_s,ssa_s,asf_s 103 real*8 ext_f,sca_f,ssa_f,asf_f 102 ! Variables for haze optics 103 character(len=200) file_path 104 logical file_ok 105 integer dumch 106 real*8 dumwvl 107 108 real*8 m3as,m3af 109 real*8 dtauaer_s,dtauaer_f 110 real*8,save :: rhoaer_s(L_NSPECTV),ssa_s(L_NSPECTV),asf_s(L_NSPECTV) 111 real*8,save :: rhoaer_f(L_NSPECTV),ssa_f(L_NSPECTV),asf_f(L_NSPECTV) 112 !$OMP THREADPRIVATE(rhoaer_s,rhoaer_f,ssa_s,ssa_f,asf_s,asf_f) 113 104 114 logical,save :: firstcall=.true. 105 115 !$OMP THREADPRIVATE(firstcall) 106 116 107 117 … … 114 124 ! Some initialisation beacause there's a pb with disr_haze at the limits (nw=1) 115 125 ! I should check this - For now we set vars to zero : better than nans - JVO 2017 116 117 dhaze_t(:,:) = 0. 118 ssa_t(:,:) = 0. 119 asf_t(:,:) = 0. 120 126 DHAZE_T(:,:) = 0.0 127 SSA_T(:,:) = 0.0 128 ASF_T(:,:) = 0.0 129 130 ! Load tabulated haze optical properties if needed. 131 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 132 IF (firstcall .AND. callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN 133 OPEN(12,file=TRIM(haze_opt_file),form='formatted') ! The file has been inquired in physiq_mod firstcall 134 READ(12,*) ! dummy header 135 DO NW=1,L_NSPECTI 136 READ(12,*) ! there's IR 1st 137 ENDDO 138 DO NW=1,L_NSPECTV 139 READ(12,*) dumch, dumwvl, rhoaer_f(nw), ssa_f(nw), asf_f(nw), rhoaer_s(nw), ssa_s(nw), asf_s(nw) 140 ENDDO 141 CLOSE(12) 142 ENDIF 121 143 122 144 !======================================================================= … … 132 154 do K=2,L_LEVELS 133 155 134 ilay = k / 2 ! int. arithmetic => gives the gcm layer index156 ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed) 135 157 136 158 DPR(k) = PLEV(K)-PLEV(K-1) … … 150 172 ! Rayleigh scattering 151 173 do NW=1,L_NSPECTV 152 TRAY(1:4,NW) = 1 d-30174 TRAY(1:4,NW) = 1.d-30 153 175 do K=5,L_LEVELS 154 176 TRAY(K,NW) = TAURAY(NW) * DPR(K) … … 159 181 do K=2,L_LEVELS 160 182 161 ilay = k / 2 ! int. arithmetic => gives the gcm layer index183 ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed) 162 184 163 185 do NW=1,L_NSPECTV 164 186 165 ! Optical coupling of YAMMS is plugged but inactivated (if false) for now166 ! as long as the microphysics only isn't fully debugged -- JVO 01/18167 187 IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN 168 m0as = pqmo(ilay,1) 169 m3as = pqmo(ilay,2) 170 m0af = pqmo(ilay,3) 171 m3af = pqmo(ilay,4) 172 173 IF (.NOT.mmp_sph_optics_vis(m0as,m3as,nw,ext_s,sca_s,ssa_s,asf_s)) & 174 CALL abort_gcm("optcv", "Fatal error in mmp_sph_optics_vis", 12) 175 IF (.NOT.mmp_fra_optics_vis(m0af,m3af,nw,ext_f,sca_f,ssa_f,asf_f)) & 176 CALL abort_gcm("optcv", "Fatal error in mmp_fra_optics_vis", 12) 177 dhaze_T(k,nw) = ext_s+ext_f 178 SSA_T(k,nw) = (sca_s+sca_f)/dhaze_T(k,nw) 179 ASF_T(k,nw) = (asf_s*sca_s + asf_f*sca_f) /(sca_s+sca_f) 188 m3as = pqmo(ilay,2) / 2.0 189 m3af = pqmo(ilay,4) / 2.0 190 191 IF ( ilay .lt. 18 ) THEN 192 m3as = pqmo(18,2) / 2.0 * dz(k) / dz(18) 193 m3af = pqmo(18,4) / 2.0 * dz(k) / dz(18) 194 ENDIF 195 196 dtauaer_s = m3as*rhoaer_s(nw) 197 dtauaer_f = m3af*rhoaer_f(nw) 198 DHAZE_T(k,nw) = dtauaer_s + dtauaer_f 199 200 IF ( dtauaer_s + dtauaer_f .GT. 1.D-30 ) THEN 201 SSA_T(k,nw) = ( dtauaer_s*ssa_s(nw) + dtauaer_f*ssa_f(nw) ) / ( dtauaer_s+dtauaer_f ) 202 ASF_T(k,nw) = ( dtauaer_s*ssa_s(nw)*asf_s(nw) + dtauaer_f*ssa_f(nw)*asf_f(nw) ) & 203 / ( ssa_s(nw)*dtauaer_s + ssa_f(nw)*dtauaer_f ) 204 ELSE 205 DHAZE_T(k,nw) = 0.D0 206 SSA_T(k,nw) = 1.0 207 ASF_T(k,nw) = 1.0 208 ENDIF 209 180 210 IF (callclouds.and.firstcall) & 181 211 WRITE(*,*) 'WARNING: In optcv, optical properties & … … 183 213 ELSE 184 214 ! Call fixed vertical haze profile of extinction - same for all columns 185 call disr_haze(dz(k),plev(k),wnov(nw), dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))186 if (seashaze) dhaze_T(k,nw) = dhaze_T(k,nw)*seashazefact(k)215 call disr_haze(dz(k),plev(k),wnov(nw),DHAZE_T(k,nw),SSA_T(k,nw),ASF_T(k,nw)) 216 if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k) 187 217 ENDIF 188 218 189 219 !JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR 190 191 192 193 194 if (k<5) dhaze_T(K,:) = 0.0220 ! but visible does not handle very well diffusion in first layer. 221 ! The tauaero and tauray are thus set to 0 (a small value for rayleigh because the code crashes otherwise) 222 ! in the 4 first semilayers in optcv, but not optci. 223 ! This solves random variations of the sw heating at the model top. 224 if (k<5) DHAZE_T(K,:) = 0.0 195 225 196 226 DRAYAER = TRAY(K,NW) … … 316 346 ! Haze scattering 317 347 !JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR 318 319 320 348 ! but not in the visible 349 ! The dhaze_s is thus set to 0 in the 4 first semilayers in optcv, but not optci. 350 ! This solves random variations of the sw heating at the model top. 321 351 DO NW=1,L_NSPECTV 322 352 DHAZES_T(1:4,NW) = 0.d0 323 353 DO K=5,L_LEVELS 324 354 DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW) ! effect of scattering albedo on haze … … 386 416 END DO ! end full gauss loop 387 417 388 389 418 if(firstcall) firstcall = .false. 390 419 -
trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90
r2138 r2242 21 21 use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen 22 22 use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat 23 use datafile_mod, only: datadir, corrkdir, banddir 23 use datafile_mod, only: datadir, corrkdir, banddir, haze_opt_file 24 24 use geometry_mod, only: latitude, longitude, cell_area 25 25 USE comgeomfi_h, only: totarea, totarea_planet … … 402 402 real,dimension(ngrid,nlayer,nq) :: dtpq ! (temporary in 01/18) 403 403 404 logical file_ok 404 405 405 406 !----------------------------------------------------------------------------- … … 417 418 REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: g3d !! Latitude-Altitude depending gravitational acceleration (m.s-2). 418 419 REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: temp !! Temperature at the center of each layer (K). 419 REAL(kind=8), DIMENSION(:,:,:), INTENT(IN) :: pq !! Tracers (\( kg.kg^{-1}}\)).420 REAL(kind=8), DIMENSION(:,:,:), INTENT(IN) :: zdqfi !! Tendency from former processes for tracers (\( kg.kg^{-1}}\)).421 REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdq !! Microphysical tendency for tracers (\( kg.kg^{-1}}\)).420 REAL(kind=8), DIMENSION(:,:,:), INTENT(IN) :: pq !! Tracers (\(X.kg^{-1}}\)). 421 REAL(kind=8), DIMENSION(:,:,:), INTENT(IN) :: zdqfi !! Tendency from former processes for tracers (\(X.kg^{-1}}\)). 422 REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdq !! Microphysical tendency for tracers (\(X.kg^{-1}}\)). 422 423 END SUBROUTINE calmufi 423 424 END INTERFACE … … 484 485 ! Variables set to 0 485 486 ! ~~~~~~~~~~~~~~~~~~ 486 dtrad(:,:) = 0. 0487 fluxrad(:) = 0. 0488 zdtsw(:,:) = 0. 0489 zdtlw(:,:) = 0. 0487 dtrad(:,:) = 0.D0 488 fluxrad(:) = 0.D0 489 zdtsw(:,:) = 0.D0 490 zdtlw(:,:) = 0.D0 490 491 491 492 ! Initialize setup for correlated-k radiative transfer 492 ! JVO 17 : Was in callcorrk firstcall, but we need spectral intervals for microphysics 493 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 493 ! JVO 17 : Was in callcorrk firstcall, but we need spectral intervals for microphysics. 494 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 494 495 495 496 if (corrk) then … … 512 513 int_dtauv(:,:,:) = 0.D0 513 514 515 IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN 516 haze_opt_file=trim(datadir)//'/HAZE_OPTIC_'//trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))//'.DAT' 517 inquire(file=trim(haze_opt_file),exist=file_ok) 518 if(.not.file_ok) then 519 write(*,*) 'The file ',TRIM(haze_opt_file),' with the haze optical properties' 520 write(*,*) 'was not found by optci.F90 ! Check in ', TRIM(datadir) 521 write(*,*) 'that you have the one corresponding to the given spectral resolution !!' 522 write(*,*) 'Meanwhile I abort ...' 523 call abort 524 endif 525 ENDIF 526 514 527 endif 515 528 … … 536 549 537 550 IF ( callmufi ) THEN 538 ! WARNING JB (27/03/2018): inimufi AND mmp_initialize_opticsshould be wrapped in an OMP SINGLE statement.551 ! WARNING JB (27/03/2018): inimufi should be wrapped in an OMP SINGLE statement. 539 552 call inimufi(ptimestep) 540 ! Optical coupling of YAMMS is plugged but inactivated for now541 ! as long as the microphysics only isn't fully debugged -- JVO 01/18542 ! NOTE JB (03/18): mmp_optic_file is initialized in inimufi => mmp_initialize_optics must be called AFTER inimufi543 IF (.NOT.uncoupl_optic_haze) call mmp_initialize_optics(mmp_optic_file)544 553 545 554 ! initialize microphysics diagnostics arrays. … … 547 556 548 557 ENDIF 549 550 558 551 559 #ifdef CPP_XIOS … … 582 590 ! Initialize albedo calculation. 583 591 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 584 albedo(:,:)=0. 0585 albedo_bareground(:)=0. 0592 albedo(:,:)=0.D0 593 albedo_bareground(:)=0.D0 586 594 call surfini(ngrid,nq,qsurf,albedo,albedo_bareground) 587 595 … … 660 668 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 661 669 662 pdt(:,:) = 0. 0663 zdtsurf(:) = 0. 0664 pdq(:,:,:) = 0. 0665 dqsurf(:,:) = 0. 0666 pdu(:,:) = 0. 0667 pdv(:,:) = 0. 0668 pdpsrf(:) = 0. 0669 zflubid(:) = 0. 0670 taux(:) = 0. 0671 tauy(:) = 0. 0670 pdt(:,:) = 0.D0 671 zdtsurf(:) = 0.D0 672 pdq(:,:,:) = 0.D0 673 dqsurf(:,:) = 0.D0 674 pdu(:,:) = 0.D0 675 pdv(:,:) = 0.D0 676 pdpsrf(:) = 0.D0 677 zflubid(:) = 0.D0 678 taux(:) = 0.D0 679 tauy(:) = 0.D0 672 680 673 681 zday=pday+ptime ! Compute time, in sols (and fraction thereof). … … 906 914 fluxtop_lw(:) = emis(:)*sigma*tsurf(:)**4 907 915 908 dtrad(:,:)=0. 0 ! no atmospheric radiative heating916 dtrad(:,:)=0.D0 ! no atmospheric radiative heating 909 917 910 918 endif ! end of corrk … … 1044 1052 1045 1053 zdh(:,:) = pdt(:,:)/zpopsk(:,:) 1046 zduadj(:,:)=0. 01047 zdvadj(:,:)=0. 01048 zdhadj(:,:)=0. 01054 zduadj(:,:)=0.D0 1055 zdvadj(:,:)=0.D0 1056 zdhadj(:,:)=0.D0 1049 1057 1050 1058 … … 1088 1096 1089 1097 if (callmufi) then 1090 zzlev(:,nlayer+1)=zzlay(:,nlayer)+(zzlay(:,nlayer)-zzlev(:,nlayer)) ! JVO 19 : We assume zzlev isn't reused later on 1098 1099 zzlev(:,nlayer+1)=zzlay(:,nlayer)+(zzlay(:,nlayer)-zzlev(:,nlayer)) ! JVO 19 : We assume zzlev isn't reused later on (could be done cleaner) 1100 1091 1101 #ifdef USE_QTEST 1092 dtpq(:,:,:) = 0. 0 ! we want tpq to go only through mufi1102 dtpq(:,:,:) = 0.D0 ! we want tpq to go only through mufi 1093 1103 call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,tpq,dtpq,zdqmufi) 1094 1104 tpq(:,:,:) = tpq(:,:,:) + zdqmufi(:,:,:)*ptimestep ! only manipulation of tpq->*ptimestep here 1095 1105 #else 1106 1096 1107 call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,pq,pdq,zdqmufi) 1108 1097 1109 pdq(:,:,:) = pdq(:,:,:) + zdqmufi(:,:,:) 1110 1111 ! Sanity check ( way safer to be done here rather than within YAMMS ) 1112 ! Important : the sanity check intentionally include the former processes tendency ! 1113 ! NB : Despite this sanity check there might be still some unphysical values going through : 1114 ! - Negatives, but harmless as it will be only for the output files 1115 ! just remove them in post-proc. 1116 ! - Weird unphysical ratio of m0 and m3, ok for now, but take care of them if 1117 ! you want to compute optics from radii. 1118 WHERE ( (pq(:,:,1)+pdq(:,:,1)*ptimestep < 0.D0) .OR. (pq(:,:,2)+pdq(:,:,2)*ptimestep < 0.D0) ) 1119 pdq(:,:,1) = (epsilon(1.0)-1.D0)*pq(:,:,1)/ptimestep 1120 pdq(:,:,2) = (epsilon(1.0)-1.D0)*pq(:,:,2)/ptimestep 1121 ENDWHERE 1122 WHERE ( (pq(:,:,3)+pdq(:,:,3)*ptimestep < 0.D0) .OR. (pq(:,:,4)+pdq(:,:,4)*ptimestep < 0.D0) ) 1123 pdq(:,:,3) = (epsilon(1.0)-1.D0)*pq(:,:,3)/ptimestep 1124 pdq(:,:,4) = (epsilon(1.0)-1.D0)*pq(:,:,4)/ptimestep 1125 ENDWHERE 1098 1126 #endif 1099 1127 … … 1103 1131 call calmufi(ptimestep,zplevbar,zzlevbar,zplaybar,zzlaybar, & 1104 1132 gzlat,ztfibar,zqfibar,zdqfibar,zdqmufibar) 1133 ! TODO : Add a sanity check here ! 1105 1134 endif 1106 1135 … … 1215 1244 ! v. Updates and positivity check 1216 1245 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1217 zdqchi(:,:,:) = 0.D0 ! -> for the microphysical tracers we need to update to 0 each time ;-)1246 zdqchi(:,:,:) = 0.D0 ! -> dycchi is saved but for the nmicro tracers we must update to 0 at each step 1218 1247 1219 1248 do iq=1,nkim … … 1315 1344 1316 1345 if(firstcall)then 1317 zdtdyn(:,:)=0. 01318 zdudyn(:,:)=0. 01346 zdtdyn(:,:)=0.D0 1347 zdudyn(:,:)=0.D0 1319 1348 endif 1320 1349 … … 1674 1703 CALL send_xios_field("dqcond_C3H6",dyccond(:,:,17+nmicro)) 1675 1704 CALL send_xios_field("dqcond_C4H4",dyccond(:,:,21+nmicro)) 1676 CALL send_xios_field("dqcond_CH3CCH",dyccond(:,:,2 4+nmicro))1677 CALL send_xios_field("dqcond_C3H8",dyccond(:,:,2 5+nmicro))1678 CALL send_xios_field("dqcond_C4H2",dyccond(:,:,2 6+nmicro))1679 CALL send_xios_field("dqcond_C4H6",dyccond(:,:,2 7+nmicro))1680 CALL send_xios_field("dqcond_C4H10",dyccond(:,:,2 8+nmicro))1681 CALL send_xios_field("dqcond_AC6H6",dyccond(:,:,2 9+nmicro))1682 CALL send_xios_field("dqcond_HCN",dyccond(:,:,3 6+nmicro))1683 CALL send_xios_field("dqcond_CH3CN",dyccond(:,:, 40+nmicro))1684 CALL send_xios_field("dqcond_HC3N",dyccond(:,:,4 2+nmicro))1685 CALL send_xios_field("dqcond_NCCN",dyccond(:,:,4 3+nmicro))1686 CALL send_xios_field("dqcond_C4N2",dyccond(:,:,4 4+nmicro))1705 CALL send_xios_field("dqcond_CH3CCH",dyccond(:,:,23+nmicro)) 1706 CALL send_xios_field("dqcond_C3H8",dyccond(:,:,24+nmicro)) 1707 CALL send_xios_field("dqcond_C4H2",dyccond(:,:,25+nmicro)) 1708 CALL send_xios_field("dqcond_C4H6",dyccond(:,:,26+nmicro)) 1709 CALL send_xios_field("dqcond_C4H10",dyccond(:,:,27+nmicro)) 1710 CALL send_xios_field("dqcond_AC6H6",dyccond(:,:,28+nmicro)) 1711 CALL send_xios_field("dqcond_HCN",dyccond(:,:,35+nmicro)) 1712 CALL send_xios_field("dqcond_CH3CN",dyccond(:,:,39+nmicro)) 1713 CALL send_xios_field("dqcond_HC3N",dyccond(:,:,41+nmicro)) 1714 CALL send_xios_field("dqcond_NCCN",dyccond(:,:,42+nmicro)) 1715 CALL send_xios_field("dqcond_C4N2",dyccond(:,:,43+nmicro)) 1687 1716 1688 1717 ! Pseudo-evaporation flux (mol/mol/s)
Note: See TracChangeset
for help on using the changeset viewer.
