Changeset 5450 for LMDZ6/branches/contrails/libf/phylmd
- Timestamp:
- Dec 23, 2024, 12:04:08 PM (7 hours ago)
- Location:
- LMDZ6/branches/contrails
- Files:
-
- 5 deleted
- 14 edited
- 7 copied
-
. (modified) (1 prop)
-
libf/phylmd/cvltr_scav.f90 (modified) (3 diffs)
-
libf/phylmd/dyn1d/scm.f90 (modified) (27 diffs)
-
libf/phylmd/iotd_ecrit.f90 (modified) (3 diffs)
-
libf/phylmd/iotd_ini.f90 (modified) (1 diff)
-
libf/phylmd/lmdz_lscp.f90 (modified) (15 diffs)
-
libf/phylmd/lmdz_lscp_condensation.f90 (modified) (5 diffs)
-
libf/phylmd/lmdz_lscp_ini.f90 (modified) (1 diff)
-
libf/phylmd/lmdz_surf_wind.f90 (copied) (copied from LMDZ6/trunk/libf/phylmd/lmdz_surf_wind.f90)
-
libf/phylmd/lmdz_surf_wind_ini.f90 (copied) (copied from LMDZ6/trunk/libf/phylmd/lmdz_surf_wind_ini.f90)
-
libf/phylmd/lmdz_thermcell_dq.F90 (deleted)
-
libf/phylmd/lmdz_thermcell_dq.f90 (copied) (copied from LMDZ6/trunk/libf/phylmd/lmdz_thermcell_dq.f90)
-
libf/phylmd/lmdz_thermcell_ini.f90 (modified) (4 diffs)
-
libf/phylmd/lmdz_thermcell_old.F90 (deleted)
-
libf/phylmd/lmdz_thermcell_old.f90 (copied) (copied from LMDZ6/trunk/libf/phylmd/lmdz_thermcell_old.f90)
-
libf/phylmd/lmdz_thermcell_plume.F90 (deleted)
-
libf/phylmd/lmdz_thermcell_plume.f90 (copied) (copied from LMDZ6/trunk/libf/phylmd/lmdz_thermcell_plume.f90)
-
libf/phylmd/lmdz_thermcell_plume_6A.F90 (deleted)
-
libf/phylmd/lmdz_thermcell_plume_6A.f90 (copied) (copied from LMDZ6/trunk/libf/phylmd/lmdz_thermcell_plume_6A.f90)
-
libf/phylmd/lmdz_wake_ini.F90 (deleted)
-
libf/phylmd/lmdz_wake_ini.f90 (copied) (copied from LMDZ6/trunk/libf/phylmd/lmdz_wake_ini.f90)
-
libf/phylmd/phys_output_ctrlout_mod.F90 (modified) (1 diff)
-
libf/phylmd/phys_output_mod.F90 (modified) (2 diffs)
-
libf/phylmd/phys_output_write_mod.F90 (modified) (3 diffs)
-
libf/phylmd/physiq_mod.F90 (modified) (4 diffs)
-
libf/phylmd/phytrac_mod.f90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/contrails
- Property svn:mergeinfo changed
/LMDZ6/trunk merged: 5433-5437,5440,5443,5445-5448
- Property svn:mergeinfo changed
-
LMDZ6/branches/contrails/libf/phylmd/cvltr_scav.f90
r5292 r5450 3 3 ! 4 4 SUBROUTINE cvltr_scav(pdtime, da, phi,phi2,d1a,dam, mpIN,epIN, & 5 sigd,sij,wght_cvfd,clw,elij,epmlmMm,eplaMm, & 6 pmflxrIN,pmflxsIN,ev,te,wdtrainA,wdtrainM, & 7 paprs,it,tr,upd,dnd,inb,icb, & 8 ccntrAA_3d,ccntrENV_3d,coefcoli_3d, & 9 dtrcv,trsptd,dtrSscav,dtrsat,dtrUscav,qDi,qPr, & 10 qPa,qMel,qTrdi,dtrcvMA,Mint, & 5 sigd,sij,wght_cvfd,clw,elij,epmlmMm,eplaMm, & 6 pmflxrIN,pmflxsIN,ev,te,wdtrainA,wdtrainM, & 7 paprs,it,tr,upd,dnd,inb,icb, & 8 ccntrAA_3d,ccntrENV_3d,coefcoli_3d, & 9 dtrcv,trsptd,dtrSscav,dtrsat,dtrUscav,flux_tr_wet, & 10 qDi,qPr, & 11 qPa,qMel,qTrdi,dtrcvMA,Mint, & 11 12 zmfd1a,zmfphi2,zmfdam) 12 13 ! … … 72 73 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrsat ! tendance trsp+sat scav 73 74 REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrUscav ! tendance du lessivage courant unsat 75 REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: flux_tr_wet ! wet deposit 74 76 ! 75 77 ! Variables locales … … 685 687 ENDDO 686 688 ENDDO 689 DO i=1, klon 690 flux_tr_wet(i,it) = (pmflxr(i,1)+pmflxs(i,1))*qPr(i,1,it)*pdtime ! wet deposit 691 ENDDO 687 692 688 693 ! test de conservation du traceur -
LMDZ6/branches/contrails/libf/phylmd/dyn1d/scm.f90
r5302 r5450 1 1 SUBROUTINE scm 2 2 3 USE flux_arp_mod_h4 5 6 3 USE flux_arp_mod_h 4 USE compbl_mod_h 5 USE clesphys_mod_h 6 USE ioipsl, only: ju2ymds, ymds2ju, ioconf_calendar,getin 7 7 USE phys_state_var_mod, ONLY : phys_state_var_init, phys_state_var_end, & 8 8 clwcon, detr_therm, & … … 33 33 USE indice_sol_mod 34 34 USE phyaqua_mod 35 ! USE mod_1D_cases_read36 35 USE mod_1D_cases_read_std 37 !USE mod_1D_amma_read38 36 USE print_control_mod, ONLY: lunout, prt_level 39 37 USE iniphysiq_mod, ONLY: iniphysiq … … 49 47 USE dimensions_mod, ONLY: iim, jjm, llm, ndm 50 48 USE dimsoil_mod_h, ONLY: nsoilmx 51 USE yomcst_mod_h52 USE tsoilnudge_mod_h53 USE fcg_gcssold_mod_h54 USE compar1d_mod_h55 USE date_cas_mod_h49 USE yomcst_mod_h 50 USE tsoilnudge_mod_h 51 USE fcg_gcssold_mod_h 52 USE compar1d_mod_h 53 USE date_cas_mod_h 56 54 implicit none 57 55 … … 86 84 87 85 integer :: kmax = llm 88 integer llm700,nq1,nq2 86 integer llm700,nq1,nq2,iflag_1d_vert_adv 89 87 INTEGER, PARAMETER :: nlev_max=1000, nqmx=1000 90 88 real timestep, frac … … 164 162 real :: dt_cooling(llm),d_t_adv(llm),d_t_nudge(llm) 165 163 real :: d_u_nudge(llm),d_v_nudge(llm) 166 ! real :: d_u_adv(llm),d_v_adv(llm)167 164 real :: d_u_age(llm),d_v_age(llm) 168 165 real :: alpha … … 174 171 REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_adv 175 172 REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_nudge 176 ! REAL, ALLOCATABLE, DIMENSION(:):: d_th_adv177 173 178 174 !--------------------------------------------------------------------- … … 182 178 real :: fder(1),snsrf(1,nbsrf),qsurfsrf(1,nbsrf) 183 179 real :: tsoil(1,nsoilmx,nbsrf) 184 ! real :: agesno(1,nbsrf)185 180 186 181 !--------------------------------------------------------------------- … … 220 215 integer jcode 221 216 INTEGER read_climoz 222 !223 217 integer :: it_end ! iteration number of the last call 224 !Al1,plev,play,phi,phis,presnivs,225 218 integer ecrit_slab_oc !1=ecrit,-1=lit,0=no file 226 219 data ecrit_slab_oc/-1/ 227 ! 220 228 221 ! if flag_inhib_forcing = 0, tendencies of forcing are added 229 222 ! <> 0, tendencies of forcing are not added … … 282 275 mth_ini_cas=1 ! pour le moment on compte depuis le debut de l'annee 283 276 call getin('time_ini',heure_ini_cas) 284 285 277 print*,'NATURE DE LA SURFACE ',nat_surf 286 ! 278 287 279 ! Initialization of the logical switch for nudging 288 280 … … 292 284 jcode = jcode/10 293 285 enddo 286 ! numerical scheme for vertical advection 287 288 iflag_1d_vert_adv=1 289 call getin('iflag_1d_vert_adv',iflag_1d_vert_adv) 290 294 291 !----------------------------------------------------------------------- 295 292 ! Definition of the run … … 398 395 allocate(d_q_adv(llm,nqtot)) 399 396 allocate(d_q_nudge(llm,nqtot)) 400 ! allocate(d_th_adv(llm))401 397 402 398 q(:,:) = 0. … … 406 402 d_q_nudge(:,:) = 0. 407 403 408 !409 404 ! No ozone climatology need be read in this pre-initialization 410 405 ! (phys_state_var_init is called again in physiq) … … 424 419 !!! Feedback forcing values for Gateaux differentiation (al1) 425 420 !!!===================================================================== 426 !! 421 427 422 qsol = qsolinp 428 423 qsurf = fq_sat(tsurf,psurf/100.) … … 433 428 rho(1)=psurf/(rd*tsurf*(1.+(rv/rd-1.)*qsurf)) 434 429 435 !436 430 !! mpl et jyg le 22/08/2012 : 437 431 !! pour que les cas a flux de surface imposes marchent … … 487 481 488 482 INCLUDE "1D_read_forc_cases.h" 489 print*,'A d_t_adv ',d_t_adv(1:20)*86400490 483 491 484 if (forcing_GCM2SCM) then … … 567 560 agesno = xagesno 568 561 tsoil(:,:,:)=tsurf 569 !------ AMMA 2e run avec modele sol et rayonnement actif (MPL 23052012)570 ! tsoil(1,1,1)=299.18571 ! tsoil(1,2,1)=300.08572 ! tsoil(1,3,1)=301.88573 ! tsoil(1,4,1)=305.48574 ! tsoil(1,5,1)=308.00575 ! tsoil(1,6,1)=308.00576 ! tsoil(1,7,1)=308.00577 ! tsoil(1,8,1)=308.00578 ! tsoil(1,9,1)=308.00579 ! tsoil(1,10,1)=308.00580 ! tsoil(1,11,1)=308.00581 562 !----------------------------------------------------------------------- 582 563 call pbl_surface_init(fder, snsrf, qsurfsrf, tsoil) … … 635 616 rvc_ancien = 0. 636 617 ENDIF 637 !jyg<638 ! Etienne: comment those lines since now the TKE is inialized in 1D_read_forc_cases639 !! pbl_tke(:,:,:)=1.e-8640 ! pbl_tke(:,:,:)=0.641 ! pbl_tke(:,2,:)=1.e-2642 !>jyg643 618 rain_fall=0. 644 619 snow_fall=0. … … 754 729 755 730 call phys_state_var_end 756 !Al1757 731 if (restart) then 758 732 print*,'call to restart dyn 1d' … … 762 736 print*,'fnday,annee_ref,day_ref,day_ini', & 763 737 & fnday,annee_ref,day_ref,day_ini 764 !** call ymds2ju(annee_ref,mois,day_ini,heure,day)765 738 day = day_ini 766 739 day_end = day_ini + nday … … 780 753 ! raz for safety 781 754 do l=1,llm 782 d_q_vert_adv(l, 1) = 0.755 d_q_vert_adv(l,:) = 0. 783 756 enddo 784 757 endif … … 819 792 !--------------------------------------------------------------------- 820 793 phis(1)=zsurf*RG 821 ! phi(1)=phis(1)+RD*temp(1)*(plev(1)-play(1))/(.5*(plev(1)+play(1)))822 794 823 795 ! Calculate geopotential from the ground surface since phi and phis are added in physiq_mod … … 842 814 z_adv=play 843 815 ENDIF 844 845 teta=temp*(pzero/play)**rkappa 816 817 ! vertical tendencies computed as d X / d t = -W d X / d z 818 819 IF (iflag_1d_vert_adv .EQ. 0) THEN 820 821 ! old centered numerical scheme 846 822 do l=2,llm-1 847 ! vertical tendencies computed as d X / d t = -W d X / d z848 823 d_u_vert_adv(l)=-w_adv(l)*(u(l+1)-u(l-1))/(z_adv(l+1)-z_adv(l-1)) 849 824 d_v_vert_adv(l)=-w_adv(l)*(v(l+1)-v(l-1))/(z_adv(l+1)-z_adv(l-1)) 850 825 ! d theta / dt = -W d theta / d z, transformed into d temp / d t dividing by (pzero/play(l))**rkappa 851 d_t_vert_adv(l)=-w_adv(l)*(te ta(l+1)-teta(l-1))/(z_adv(l+1)-z_adv(l-1)) / (pzero/play(l))**rkappa852 d_q_vert_adv(l, 1)=-w_adv(l)*(q(l+1,1)-q(l-1,1))/(z_adv(l+1)-z_adv(l-1))826 d_t_vert_adv(l)=-w_adv(l)*(temp(l+1)-temp(l-1))/(z_adv(l+1)-z_adv(l-1)) 827 d_q_vert_adv(l,:)=-w_adv(l)*(q(l+1,:)-q(l-1,:))/(z_adv(l+1)-z_adv(l-1)) 853 828 enddo 854 d_u_adv(:)=d_u_adv(:)+d_u_vert_adv(:) 855 d_v_adv(:)=d_v_adv(:)+d_v_vert_adv(:) 856 d_t_adv(:)=d_t_adv(:)+d_t_vert_adv(:) 857 d_q_adv(:,1)=d_q_adv(:,1)+d_q_vert_adv(:,1) 829 830 831 ELSE IF (iflag_1d_vert_adv .EQ. 1) THEN 832 ! upstream numerical scheme 833 834 do l=2,llm-1 835 IF ( ( ( forc_w .EQ. 1 ) .AND. ( w_adv(l) .GT. 0. ) ) .OR. & 836 ( ( forc_w .NE. 1 ) .AND. ( w_adv(l) .LT. 0. ) ) ) THEN 837 d_u_vert_adv(l)=-w_adv(l)*(u(l)-u(l-1))/(z_adv(l)-z_adv(l-1)) 838 d_v_vert_adv(l)=-w_adv(l)*(v(l)-v(l-1))/(z_adv(l)-z_adv(l-1)) 839 d_t_vert_adv(l)=-w_adv(l)*(temp(l)-temp(l-1))/(z_adv(l)-z_adv(l-1)) 840 d_q_vert_adv(l,:)=-w_adv(l)*(q(l,:)-q(l-1,:))/(z_adv(l)-z_adv(l-1)) 841 ELSE 842 d_u_vert_adv(l)=-w_adv(l)*(u(l+1)-u(l))/(z_adv(l+1)-z_adv(l)) 843 d_v_vert_adv(l)=-w_adv(l)*(v(l+1)-v(l))/(z_adv(l+1)-z_adv(l)) 844 d_t_vert_adv(l)=-w_adv(l)*(temp(l+1)-temp(l))/(z_adv(l+1)-z_adv(l)) 845 d_q_vert_adv(l,:)=-w_adv(l)*(q(l+1,:)-q(l,:))/(z_adv(l+1)-z_adv(l)) 846 ENDIF 847 enddo 848 849 850 ENDIF ! numerical scheme for vertical advection 851 852 853 IF (flag_inhib_forcing == 0) then ! if tendency of forcings should be added 854 u(1:mxcalc)=u(1:mxcalc) + timestep * d_u_vert_adv(1:mxcalc) 855 v(1:mxcalc)=v(1:mxcalc) + timestep * d_v_vert_adv(1:mxcalc) 856 q(1:mxcalc,:)=q(1:mxcalc,:) + timestep * d_q_vert_adv(1:mxcalc,:) 857 temp(1:mxcalc)=temp(1:mxcalc) + timestep * d_t_vert_adv(1:mxcalc) 858 teta=temp*(pzero/play)**rkappa 859 ENDIF 858 860 859 861 ENDIF … … 940 942 !--------------------------------------------------------------------- 941 943 ! Nudging 942 ! EV: rewrite the section to account for a time- and height-varying943 ! nudging944 944 !--------------------------------------------------------------------- 945 945 d_t_nudge(:) = 0. … … 1000 1000 1001 1001 ENDDO 1002 1003 !----------------------------------------------------------- 1004 ! horizontal forcings (advection) and nudging 1005 !----------------------------------------------------------- 1006 1007 IF (flag_inhib_forcing == 0) then ! if tendency of forcings should be added 1008 1009 u(1:mxcalc)=u(1:mxcalc) + timestep*( & 1010 & du_phys(1:mxcalc) & 1011 & +d_u_age(1:mxcalc)+d_u_adv(1:mxcalc) & 1012 & +d_u_nudge(1:mxcalc) ) 1013 v(1:mxcalc)=v(1:mxcalc) + timestep*( & 1014 & dv_phys(1:mxcalc) & 1015 & +d_v_age(1:mxcalc)+d_v_adv(1:mxcalc) & 1016 & +d_v_nudge(1:mxcalc) ) 1017 q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*( & 1018 & dq(1:mxcalc,:) & 1019 & +d_q_adv(1:mxcalc,:) & 1020 & +d_q_nudge(1:mxcalc,:) ) 1021 1022 1023 temp(1:mxcalc)=temp(1:mxcalc)+timestep*( & 1024 & dt_phys(1:mxcalc) & 1025 & +d_t_adv(1:mxcalc) & 1026 & +d_t_nudge(1:mxcalc) & 1027 & +dt_cooling(1:mxcalc)) ! Taux de chauffage ou refroid. 1002 1028 1003 1029 !--------------------------------------------------------------------- … … 1022 1048 END IF 1023 1049 1024 IF (flag_inhib_forcing == 0) then ! if tendency of forcings should be added 1025 1026 u(1:mxcalc)=u(1:mxcalc) + timestep*( & 1027 & du_phys(1:mxcalc) & 1028 & +d_u_age(1:mxcalc)+d_u_adv(1:mxcalc) & 1029 & +d_u_nudge(1:mxcalc) ) 1030 v(1:mxcalc)=v(1:mxcalc) + timestep*( & 1031 & dv_phys(1:mxcalc) & 1032 & +d_v_age(1:mxcalc)+d_v_adv(1:mxcalc) & 1033 & +d_v_nudge(1:mxcalc) ) 1034 q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*( & 1035 & dq(1:mxcalc,:) & 1036 & +d_q_adv(1:mxcalc,:) & 1037 & +d_q_nudge(1:mxcalc,:) ) 1038 1039 if (prt_level.ge.3) then 1040 print *, & 1041 & 'physiq-> temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1) ', & 1042 & temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1) 1043 print* ,'dv_phys=',dv_phys 1044 print* ,'d_v_age=',d_v_age 1045 print* ,'d_v_adv=',d_v_adv 1046 print* ,'d_v_nudge=',d_v_nudge 1047 print*, v 1048 print*, vg 1049 endif 1050 1051 temp(1:mxcalc)=temp(1:mxcalc)+timestep*( & 1052 & dt_phys(1:mxcalc) & 1053 & +d_t_adv(1:mxcalc) & 1054 & +d_t_nudge(1:mxcalc) & 1055 & +dt_cooling(1:mxcalc)) ! Taux de chauffage ou refroid. 1056 1057 1058 !======================================================================= 1059 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !! 1060 !======================================================================= 1050 1051 1052 ! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !! 1061 1053 1062 1054 teta=temp*(pzero/play)**rkappa … … 1071 1063 ENDIF 1072 1064 1073 !--------------------------------------------------------------------- 1074 ! Add large-scale tendencies (advection, etc) : 1075 !--------------------------------------------------------------------- 1076 1077 !cc nrlmd 1078 !cc tmpvar=teta 1079 !cc call advect_vert(llm,omega,timestep,tmpvar,plev) 1080 !cc 1081 !cc teta(1:mxcalc)=tmpvar(1:mxcalc) 1082 !cc tmpvar(:)=q(:,1) 1083 !cc call advect_vert(llm,omega,timestep,tmpvar,plev) 1084 !cc q(1:mxcalc,1)=tmpvar(1:mxcalc) 1085 !cc tmpvar(:)=q(:,2) 1086 !cc call advect_vert(llm,omega,timestep,tmpvar,plev) 1087 !cc q(1:mxcalc,2)=tmpvar(1:mxcalc) 1088 1089 END IF ! end if tendency of tendency should be added 1065 1066 END IF ! end if in case tendency should be added 1090 1067 1091 1068 !--------------------------------------------------------------------- … … 1101 1078 daytime = daytime+1./day_step 1102 1079 day = int(daytime+0.1/day_step) 1103 ! time = max(daytime-day,0.0)1104 !Al1&jyg: correction de bug1105 !cc time = real(mod(it,day_step))/day_step1106 1080 time = time_ini/24.+real(mod(it,day_step))/day_step 1107 1081 it=it+1 -
LMDZ6/branches/contrails/libf/phylmd/iotd_ecrit.f90
r5390 r5450 123 123 !Case of a 3D variable 124 124 !--------------------- 125 125 if (llm==lmax) then 126 126 ndim=4 127 127 corner(1)=1 … … 134 134 edges(4)=1 135 135 dim_cc=dim_coord 136 137 136 138 137 !Case of a 2D variable … … 159 158 if (ntime==1) then 160 159 ierr = nf90_redef (nid) 161 ierr = nf90_def_var(nid,nom,nf90_float,dim_cc ,varid)160 ierr = nf90_def_var(nid,nom,nf90_float,dim_cc(1:ndim),varid) 162 161 !print*,'DEF ',nom,nid,varid 163 162 ierr = nf90_enddef(nid) -
LMDZ6/branches/contrails/libf/phylmd/iotd_ini.f90
r5390 r5450 124 124 ierr=nf90_def_var(nid, "lev", nf90_float, dim_coord(3),nvarid) 125 125 ierr=nf90_put_att(nid,nvarid,"long_name","vert level") 126 if ( coordv(2)>coordv(1) ) then 126 127 if ( coordv(min(2,llm))>=coordv(1) ) then 127 128 ierr=nf90_put_att(nid,nvarid,"long_name","pseudo-alt") 128 129 ierr=nf90_put_att(nid,nvarid,'positive',"up") -
LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp.f90
r5431 r5450 268 268 REAL, DIMENSION(klon,klev) :: ctot 269 269 REAL, DIMENSION(klon,klev) :: ctot_vol 270 REAL, DIMENSION(klon) :: zqs, zdqs 270 REAL, DIMENSION(klon) :: zqs, zdqs, zqsl, zdqsl, zqsi, zdqsi 271 271 REAL :: zdelta 272 272 REAL, DIMENSION(klon) :: zdqsdT_raw … … 278 278 REAL, DIMENSION(klon) :: qcloud, qincloud_mpc 279 279 REAL, DIMENSION(klon) :: zrfl, zifl 280 REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn 280 REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn, zqnl 281 281 REAL, DIMENSION(klon) :: zoliql, zoliqi 282 282 REAL, DIMENSION(klon) :: zt 283 REAL, DIMENSION(klon) :: zfice,zneb 283 REAL, DIMENSION(klon) :: zfice,zneb, znebl 284 284 REAL, DIMENSION(klon) :: dzfice 285 285 REAL, DIMENSION(klon) :: zfice_turb, dzfice_turb … … 306 306 307 307 ! for condensation and ice supersaturation 308 REAL, DIMENSION(klon) :: qvc, shear308 REAL, DIMENSION(klon) :: qvc, qvcl, shear 309 309 REAL :: delta_z 310 310 !--Added for ice supersaturation (ok_ice_supersat) and contrails (ok_plane_contrails) … … 377 377 rain(:) = 0.0 378 378 snow(:) = 0.0 379 zfice(:)= 0.0379 zfice(:)=1.0 ! initialized at 1 as by default we assume mpc to be at ice saturation 380 380 dzfice(:)=0.0 381 381 zfice_turb(:)=0.0 … … 668 668 ! Calculation of saturation specific humidity and ice fraction 669 669 CALL calc_qsat_ecmwf(klon,Tbef,qtot,pplay(:,k),RTT,0,.false.,zqs,zdqs) 670 671 IF (iflag_icefrac .GE. 3) THEN 672 ! consider a ice weighted qs to ensure that liquid clouds at T<0 have a consistent cloud fraction 673 ! and cloud condensed water content. idea following Dietlitcher et al. 2018, GMD 674 ! we make this option works only for the physically-based and tke-depdenent param (iflag_icefrac>=1) 675 DO i=1,klon 676 CALL calc_qsat_ecmwf(klon,Tbef,qtot,pplay(:,k),RTT,1,.false.,zqsl,zdqsl) 677 CALL calc_qsat_ecmwf(klon,Tbef,qtot,pplay(:,k),RTT,2,.false.,zqsi,zdqsi) 678 zqs(i)=zfice(i)*zqsi(i)+(1.-zfice(i))*zqsl(i) 679 zdqs(i)=zfice(i)*zdqsi(i)+zqsi(i)*dzfice(i)+(1.-zfice(i))*zdqsl(i)-zqsl(i)*dzfice(i) 680 ENDDO 681 ENDIF 682 670 683 CALL calc_gammasat(klon,Tbef,qtot,pplay(:,k),gammasat,dgammasatdt) 671 684 ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs 672 685 zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:) 673 ! cloud phase determination 674 IF (iflag_t_glace.GE.4) THEN 675 ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top 676 CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop) 677 ENDIF 678 679 CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop,ztemp_cltop,zfice,dzfice) 680 686 687 ! Cloud condensation based on subgrid pdf 681 688 !--AB Activates a condensation scheme that allows for 682 689 !--ice supersaturation and contrails evolution from aviation … … 733 740 ENDIF ! .NOT. ok_ice_supersat 734 741 742 ! cloud phase determination 743 IF (iflag_icefrac .GE. 2) THEN 744 ! phase partitioning depending on temperature. activates here in the iteration process if iflag_icefrac > 2 745 CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), omega(:,k), qice_ini, ziflcld, zqn, & 746 rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice, dzfice, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k)) 747 ELSE 748 ! phase partitioning depending on temperature and eventually distance to cloud top 749 IF (iflag_t_glace.GE.4) THEN 750 ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top 751 CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop) 752 ENDIF 753 CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop,ztemp_cltop,zfice,dzfice) 754 ENDIF 755 756 735 757 DO i=1,klon 736 758 IF (keepgoing(i)) THEN … … 756 778 ENDIF 757 779 758 ! LEA_R : check formule759 780 IF ( ok_unadjusted_clouds ) THEN 760 781 !--AB We relax the saturation adjustment assumption … … 796 817 797 818 819 ! Cloud phase final determination 820 !-------------------------------- 798 821 ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top 799 822 IF (iflag_t_glace.GE.4) THEN … … 802 825 temp_cltop(:,k)=ztemp_cltop(:) 803 826 ENDIF 804 805 ! Partition function depending on temperature 827 ! Partition function depending on temperature for all clouds (shallow convective and stratiform) 806 828 CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice) 807 829 808 ! Partition function depending on tke for non shallow-convective clouds 830 ! Partition function depending on tke for non shallow-convective clouds, erase previous estimation 809 831 IF (iflag_icefrac .GE. 1) THEN 810 832 CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), omega(:,k), qice_ini, ziflcld, zqn, & … … 812 834 ENDIF 813 835 814 ! Water vapor update, Phase determination and subsequent latent heat exchange 836 ! Water vapor update, subsequent latent heat exchange for each cloud type 837 !------------------------------------------------------------------------ 815 838 DO i=1, klon 816 839 ! Overwrite phase partitioning in boundary layer mixed phase clouds when the … … 855 878 IF (lognormale(i)) THEN 856 879 zcond(i) = zqliq(i) + zqice(i) 857 zfice(i) =zfice_turb(i)880 zfice(i) = zfice_turb(i) 858 881 rhcl(i,k) = zqvapcl(i) * rneb(i,k) + (zq(i) - zqn(i)) * (1.-rneb(i,k)) 859 882 ENDIF … … 865 888 866 889 ENDIF 867 868 ! c_iso : routine that computes in-cloud supersaturation 869 ! c_iso condensation of isotopes (zcond, zsursat, zfice, zq in input) 870 890 891 871 892 ! temperature update due to phase change 872 893 zt(i) = zt(i) + (1.-zfice(i))*zcond(i) & … … 896 917 zoliql(i) = zcond(i) * ( 1. - zfice(i) ) 897 918 zoliqi(i) = zcond(i) * zfice(i) 898 ! c_iso : initialisation of zoliq* also for isotopes899 919 ENDDO 900 920 … … 984 1004 ENDIF 985 1005 986 zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/ zneb(i))987 frac_nucl(i,k)= 1.- zneb(i)*zfrac_lessi1006 zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb)) 1007 frac_nucl(i,k)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi 988 1008 989 1009 ! Nucleation with a factor of -1 instead of -0.5 990 zfrac_lessi = 1. - EXP(-zprec_cond(i)/ zneb(i))1010 zfrac_lessi = 1. - EXP(-zprec_cond(i)/MAX(rneb(i,k),seuil_neb)) 991 1011 992 1012 ENDIF … … 1008 1028 ENDIF 1009 1029 1010 zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/ zneb(i))1011 frac_impa(i,kk)= 1.- zneb(i)*zfrac_lessi1030 zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb)) 1031 frac_impa(i,kk)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi 1012 1032 1013 1033 ENDIF -
LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90
r5406 r5450 230 230 REAL :: dqvc_mix_sub, dqvc_mix_issr 231 231 REAL :: dqt_mix 232 REAL :: a_mix, bovera, N_cld_mix, L_mix232 REAL :: a_mix, bovera, Povera, N_cld_mix, L_mix 233 233 REAL :: envfra_mix, cldfra_mix 234 234 REAL :: L_shear, shear_fra … … 292 292 !--formed elsewhere) 293 293 IF (keepgoing(i)) THEN 294 295 !--Initialisation 296 issrfra(i) = 0. 297 qissr(i) = 0. 294 298 295 299 !--If the temperature is higher than the threshold below which … … 364 368 dqvc_mix(i) = 0. 365 369 366 issrfra(i) = 0.367 qissr(i) = 0.368 subfra(i) = 0.369 qsub(i) = 0.370 371 370 !--Initialisation of the cell properties 372 371 !--Dry density [kg/m3] … … 665 664 !--Clouds within the mesh are assumed to be ellipses. The length of the 666 665 !--semi-major axis is a and the length of the semi-minor axis is b. 667 !--N_cld_mix is the number of clouds within the mesh, and 666 !--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer. 667 !--In particular, it is 0 if cldfra = 1. 668 668 !--clouds_perim is the total perimeter of the clouds within the mesh, 669 669 !--not considering interfaces with other meshes (only the interfaces with clear … … 684 684 !-- b / a = bovera = MAX(0.1, cldfra) 685 685 bovera = MAX(0.1, cldfra(i)) 686 !--P / a is a function of b / a only, that we can calculate 687 !-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) ) 688 Povera = RPI * ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) ) 686 689 !--The clouds perimeter is imposed using the formula from Morcrette 2012, 687 690 !--based on observations. 688 !-- clouds_perim_normalized = alpha * cldfra * ( 1. - cldfra ) = clouds_perim / ( norm_constant * cell_area ) 689 !-- 690 !--With all this, we have 691 !-- a = SQRT( cell_area * cldfra / ( b / a * N_cld_mix * RPI ) ) 692 !-- P = RPI * a * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) ) 693 !--and therefore, using the perimeter 694 !-- alpha * cldfra * ( 1. - cldfra ) * norm_constant * cell_area 695 !-- = N_cld_mix * RPI & 696 !-- * SQRT( cell_area * cldfra / ( b / a * N_cld_mix * RPI ) ) & 697 !-- * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) ) 698 !--and finally 699 N_cld_mix = coef_mixing_lscp * cldfra(i) * ( 1. - cldfra(i) )**2 & 700 * cell_area(i) * bovera / RPI & 701 / ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )**2 702 !--where coef_mix_lscp = ( alpha * norm_constant )**2 703 !--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer 704 !--In particular, it is 0 if cldfra = 1 705 a_mix = SQRT( cell_area(i) * cldfra(i) / bovera / N_cld_mix / RPI ) 691 !-- clouds_perim / cell_area = N_cld_mix * ( P / a * a ) / cell_area = coef_mix_lscp * cldfra * ( 1. - cldfra ) 692 !--With cldfra = a * ( b / a * a ) * RPI * N_cld_mix / cell_area, we have: 693 !-- cldfra = a * b / a * RPI / (P / a) * coef_mix_lscp * cldfra * ( 1. - cldfra ) 694 !-- a = (P / a) / ( coef_mix_lscp * RPI * ( 1. - cldfra ) * (b / a) ) 695 a_mix = Povera / coef_mixing_lscp / RPI / ( 1. - cldfra(i) ) / bovera 696 !--and finally, 697 !-- N_cld_mix = coef_mix_lscp * cldfra * ( 1. - cldfra ) * cell_area / ( P / a * a ) 698 N_cld_mix = coef_mixing_lscp * cldfra(i) * ( 1. - cldfra(i) ) * cell_area(i) & 699 / Povera / a_mix 706 700 707 701 !--The time required for turbulent diffusion to homogenize a region of size -
LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_ini.f90
r5412 r5450 195 195 !$OMP THREADPRIVATE(delta_hetfreez) 196 196 197 REAL, SAVE, PROTECTED :: coef_mixing_lscp= 9.E-8! [-] tuning coefficient for the mixing process197 REAL, SAVE, PROTECTED :: coef_mixing_lscp=1.E-3 ! [-] tuning coefficient for the mixing process 198 198 !$OMP THREADPRIVATE(coef_mixing_lscp) 199 199 -
LMDZ6/branches/contrails/libf/phylmd/lmdz_thermcell_ini.f90
r5400 r5450 8 8 integer, protected :: dvdq=1,dqimpl=-1,prt_level=0,lunout 9 9 real , protected :: RG,RD,RCPD,RKAPPA,RLVTT,RLvCp,RETV 10 11 10 12 11 13 !$OMP THREADPRIVATE(dvdq,dqimpl,prt_level,lunout) … … 47 49 !$OMP THREADPRIVATE(iflag_thermals_tenv) 48 50 51 integer, protected :: thermals_subsid_advect_more_than_one=1 52 character*6, protected :: thermals_subsid_advect_scheme = 'upwind' ! or 'center' 53 54 !$OMP THREADPRIVATE(thermals_subsid_advect_scheme,thermals_subsid_advect_more_than_one) 49 55 50 56 CONTAINS … … 101 107 CALL getin_p('thermals_flag_alim',thermals_flag_alim) 102 108 CALL getin_p('iflag_thermals_tenv',iflag_thermals_tenv) 103 109 CALL getin_p('thermals_subsid_advect_scheme',thermals_subsid_advect_scheme) 110 CALL getin_p('thermals_subsid_advect_more_than_one',thermals_subsid_advect_more_than_one) 104 111 105 112 … … 134 141 write(lunout,*) 'thermcell_ini ,thermals_flag_alim =', thermals_flag_alim 135 142 write(lunout,*) 'thermcell_ini ,iflag_thermals_tenv =', iflag_thermals_tenv 143 write(lunout,*) 'thermcell_ini ,thermals_subsid_advect_scheme=',thermals_subsid_advect_scheme 144 write(lunout,*) 'thermcell_ini ,thermals_subsid_advect_more_than_one=',thermals_subsid_advect_more_than_one 136 145 137 146 RETURN -
LMDZ6/branches/contrails/libf/phylmd/phys_output_ctrlout_mod.F90
r5383 r5450 2001 2001 TYPE(ctrl_out), SAVE, ALLOCATABLE :: o_dtr_sat(:) 2002 2002 TYPE(ctrl_out), SAVE, ALLOCATABLE :: o_dtr_uscav(:) 2003 TYPE(ctrl_out), SAVE, ALLOCATABLE :: o_dtr_wet_con(:) 2003 2004 TYPE(ctrl_out), SAVE, ALLOCATABLE :: o_dtr_dry(:) 2004 2005 -
LMDZ6/branches/contrails/libf/phylmd/phys_output_mod.F90
r5394 r5450 172 172 ALLOCATE(o_dtr_evapls(nqtot),o_dtr_ls(nqtot),o_dtr_trsp(nqtot)) 173 173 ALLOCATE(o_dtr_sscav(nqtot),o_dtr_sat(nqtot),o_dtr_uscav(nqtot)) 174 ALLOCATE(o_dtr_wet_con(nqtot)) 174 175 ALLOCATE(o_dtr_dry(nqtot),o_dtr_vdf(nqtot)) 175 176 IF (CPPKEY_STRATAER) THEN … … 540 541 tnam = TRIM(dn)//'uscav'; o_dtr_uscav (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)]) 541 542 543 lnam = 'tracer convective wet deposition'//TRIM(tracers(iq)%longName) 544 tnam = TRIM(dn)//'wet_con'; o_dtr_wet_con (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)]) 542 545 lnam = 'tracer tendency dry deposition'//TRIM(tracers(iq)%longName) 543 546 tnam = 'cum'//TRIM(dn)//'dry'; o_dtr_dry (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)]) -
LMDZ6/branches/contrails/libf/phylmd/phys_output_write_mod.F90
r5384 r5450 6 6 USE phytrac_mod, ONLY : d_tr_cl, d_tr_th, d_tr_cv, d_tr_lessi_impa, & 7 7 d_tr_lessi_nucl, d_tr_insc, d_tr_bcscav, d_tr_evapls, d_tr_ls, & 8 d_tr_trsp, d_tr_sscav, d_tr_sat, d_tr_uscav, flux_tr_ dry8 d_tr_trsp, d_tr_sscav, d_tr_sat, d_tr_uscav, flux_tr_wet, flux_tr_dry 9 9 10 10 ! Author: Abderrahmane IDELKADI (original include file) … … 189 189 o_dtr_insc, o_dtr_bcscav, o_dtr_evapls, & 190 190 o_dtr_ls, o_dtr_trsp, o_dtr_sscav, o_dtr_dry, & 191 o_dtr_sat, o_dtr_uscav, o_trac_cum, o_du_gwd_rando, o_dv_gwd_rando, & 191 o_dtr_sat, o_dtr_uscav, o_dtr_wet_con, & 192 o_trac_cum, o_du_gwd_rando, o_dv_gwd_rando, & 192 193 o_ustr_gwd_hines,o_vstr_gwd_hines,o_ustr_gwd_rando,o_vstr_gwd_rando, & 193 194 o_ustr_gwd_front,o_vstr_gwd_front, & … … 2856 2857 CALL histwrite_phy(o_dtr_uscav(itr),d_tr_uscav(:,:,itr)) 2857 2858 !--2D fields 2859 CALL histwrite_phy(o_dtr_wet_con(itr), flux_tr_wet(:,itr)) 2858 2860 CALL histwrite_phy(o_dtr_dry(itr), flux_tr_dry(:,itr)) 2859 2861 zx_tmp_fi2d=0. -
LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90
r5431 r5450 80 80 USE lmdz_call_blowing_snow, ONLY : call_blowing_snow_sublim_sedim 81 81 USE lmdz_wake_ini, ONLY : wake_ini 82 USE lmdz_surf_wind_ini, ONLY : surf_wind_ini, iflag_surf_wind 83 USE lmdz_surf_wind, ONLY : surf_wind 82 84 USE yamada_ini_mod, ONLY : yamada_ini 83 85 USE lmdz_atke_turbulence_ini, ONLY : atke_ini … … 1263 1265 CHARACTER(len=512) :: namelist_ecrad_file 1264 1266 1267 1268 ! Subgrid scale wind : 1269 ! Need to be allocatable/save because the number of bin is not known (provided by surf_wind_ini) 1270 integer, save :: nsrfwnd=1 1271 real, dimension(:,:), allocatable, save :: surf_wind_value, surf_wind_proba ! module and probability of sugrdi wind wind sample 1272 !$OMP THREADPRIVATE(nsrfwnd,surf_wind_value, surf_wind_proba) 1273 1274 1275 1265 1276 !======================================================================! 1266 1277 ! Bifurcation vers un nouveau moniteur physique pour experimenter ! … … 1818 1829 CALL iniradia(klon,klev,paprs(1,1:klev+1)) 1819 1830 1831 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1832 CALL surf_wind_ini(klon,lunout) 1833 CALL getin_p('nsrfwnd',nsrfwnd) 1834 allocate(surf_wind_value(klon,nsrfwnd),surf_wind_proba(klon,nsrfwnd)) 1835 1820 1836 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1821 1837 CALL wake_ini(rg,rd,rv,prt_level) … … 3802 3818 ! 3803 3819 !=================================================================== 3820 ! Computation of subrgid scale near-surface wind distribution 3821 call surf_wind(klon,nsrfwnd,u10m,v10m,wake_s,wake_Cstar,ustar,wstar,surf_wind_value,surf_wind_proba) 3822 3823 !=================================================================== 3804 3824 ! Computation of ratqs, the width (normalized) of the subrid scale 3805 3825 ! water distribution -
LMDZ6/branches/contrails/libf/phylmd/phytrac_mod.f90
r5330 r5450 35 35 REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_sat 36 36 REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_uscav 37 REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: flux_tr_wet ! tracer wet deposit (surface) jyg 37 38 REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: qPr,qDi ! concentration tra dans pluie,air descente insaturee 38 39 REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: qPa,qMel … … 47 48 48 49 !$OMP THREADPRIVATE(qPa,qMel,qTrdi,dtrcvMA,d_tr_th,d_tr_lessi_impa,d_tr_lessi_nucl) 49 !$OMP THREADPRIVATE(d_tr_trsp,d_tr_sscav,d_tr_sat,d_tr_uscav, qPr,qDi)50 !$OMP THREADPRIVATE(d_tr_trsp,d_tr_sscav,d_tr_sat,d_tr_uscav,flux_tr_wet,qPr,qDi) 50 51 !$OMP THREADPRIVATE(d_tr_insc,d_tr_bcscav,d_tr_evapls,d_tr_ls,qPrls) 51 52 !$OMP THREADPRIVATE(d_tr_cl,d_tr_dry,flux_tr_dry,d_tr_dec,d_tr_cv) … … 68 69 ALLOCATE(d_tr_sscav(klon,klev,nbtr),d_tr_sat(klon,klev,nbtr)) 69 70 ALLOCATE(d_tr_uscav(klon,klev,nbtr),qPr(klon,klev,nbtr),qDi(klon,klev,nbtr)) 71 ALLOCATE(flux_tr_wet(klon,nbtr)) 70 72 ALLOCATE(qPa(klon,klev,nbtr),qMel(klon,klev,nbtr)) 71 73 ALLOCATE(qTrdi(klon,klev,nbtr),dtrcvMA(klon,klev,nbtr)) … … 408 410 d_tr_dry(i,it)=0. 409 411 flux_tr_dry(i,it)=0. 412 flux_tr_wet(i,it)=0. 410 413 ENDDO 411 414 ENDDO … … 697 700 !--with the full array tr_seri even if only item it is processed 698 701 699 CALL cvltr_scav(pdtphys, da, phi,phi2,d1a,dam, mp,ep, & 700 sigd,sij,wght_cvfd,clw,elij,epmlmMm,eplaMm, & 701 pmflxr,pmflxs,evap,t_seri,wdtrainA,wdtrainM, & 702 paprs,it,tr_seri,upwd,dnwd,itop_con,ibas_con, & 703 ccntrAA_3d,ccntrENV_3d,coefcoli_3d, & 704 d_tr_cv,d_tr_trsp,d_tr_sscav,d_tr_sat,d_tr_uscav,qDi,qPr,& 705 qPa,qMel,qTrdi,dtrcvMA,Mint, & 702 CALL cvltr_scav(pdtphys, da, phi,phi2,d1a,dam, mp,ep, & 703 sigd,sij,wght_cvfd,clw,elij,epmlmMm,eplaMm, & 704 pmflxr,pmflxs,evap,t_seri,wdtrainA,wdtrainM, & 705 paprs,it,tr_seri,upwd,dnwd,itop_con,ibas_con, & 706 ccntrAA_3d,ccntrENV_3d,coefcoli_3d, & 707 d_tr_cv,d_tr_trsp,d_tr_sscav,d_tr_sat,d_tr_uscav,flux_tr_wet, & 708 qDi,qPr, & 709 qPa,qMel,qTrdi,dtrcvMA,Mint, & 706 710 zmfd1a,zmfphi2,zmfdam) 707 711
Note: See TracChangeset
for help on using the changeset viewer.
