Ignore:
Timestamp:
Dec 23, 2024, 12:04:08 PM (9 hours ago)
Author:
aborella
Message:

Merge with trunk

Location:
LMDZ6/branches/contrails
Files:
10 deleted
15 edited
12 copied

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/contrails

  • LMDZ6/branches/contrails/libf/misc/lmdz_cosp_wrappers.F90

    r5329 r5450  
    55END SUBROUTINE lmdz_cosp_wrapper_abort
    66
    7 #ifndef COSP
     7#ifndef CPP_COSP
    88
    99SUBROUTINE phys_cosp(itap, dtime, freq_cosp, &
    … …  
    4444#endif
    4545
    46 #ifndef COSP2
     46#ifndef CPP_COSP2
    4747
    4848subroutine phys_cosp2( itap,dtime,freq_cosp, ok_mensuelCOSP, ok_journeCOSP,   &
    … …  
    8585#endif
    8686
    87 #ifndef COSPV2
     87#ifndef CPP_COSPV2
    8888!
    8989subroutine lmdz_cosp_interface(itap, dtime, freq_cosp, ok_mensuelCOSP, ok_journeCOSP, &

  • LMDZ6/branches/contrails/libf/phylmd/cvltr_scav.f90

    r5292 r5450  
    33!
    44SUBROUTINE 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,                       &
    1112     zmfd1a,zmfphi2,zmfdam)
    1213  !
    … …  
    7273  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: dtrsat    ! tendance trsp+sat scav
    7374  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
    7476  !
    7577  ! Variables locales
    … …  
    685687     ENDDO
    686688  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
    687692
    688693  ! test de conservation du traceur

  • LMDZ6/branches/contrails/libf/phylmd/dyn1d/scm.f90

    r5302 r5450  
    11SUBROUTINE scm
    22
    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
     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
    77   USE phys_state_var_mod, ONLY : phys_state_var_init, phys_state_var_end, &
    88       clwcon, detr_therm, &
    … …  
    3333   USE indice_sol_mod
    3434   USE phyaqua_mod
    35 !  USE mod_1D_cases_read
    3635   USE mod_1D_cases_read_std
    37    !USE mod_1D_amma_read
    3836   USE print_control_mod, ONLY: lunout, prt_level
    3937   USE iniphysiq_mod, ONLY: iniphysiq
    … …  
    4947   USE dimensions_mod, ONLY: iim, jjm, llm, ndm
    5048   USE dimsoil_mod_h, ONLY: nsoilmx
    51   USE yomcst_mod_h
    52   USE tsoilnudge_mod_h
    53   USE fcg_gcssold_mod_h
    54   USE compar1d_mod_h
    55   USE date_cas_mod_h
     49   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
    5654implicit none
    5755
    … …  
    8684
    8785        integer :: kmax = llm
    88         integer llm700,nq1,nq2
     86        integer llm700,nq1,nq2,iflag_1d_vert_adv
    8987        INTEGER, PARAMETER :: nlev_max=1000, nqmx=1000
    9088        real timestep, frac
    … …  
    164162      real :: dt_cooling(llm),d_t_adv(llm),d_t_nudge(llm)
    165163      real :: d_u_nudge(llm),d_v_nudge(llm)
    166 !      real :: d_u_adv(llm),d_v_adv(llm)
    167164      real :: d_u_age(llm),d_v_age(llm)
    168165      real :: alpha
    … …  
    174171      REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_adv
    175172      REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_nudge
    176 !      REAL, ALLOCATABLE, DIMENSION(:):: d_th_adv
    177173
    178174!---------------------------------------------------------------------
    … …  
    182178      real :: fder(1),snsrf(1,nbsrf),qsurfsrf(1,nbsrf)
    183179      real :: tsoil(1,nsoilmx,nbsrf)
    184 !     real :: agesno(1,nbsrf)
    185180
    186181!---------------------------------------------------------------------
    … …  
    220215      integer jcode
    221216      INTEGER read_climoz
    222 !
    223217      integer :: it_end ! iteration number of the last call
    224 !Al1,plev,play,phi,phis,presnivs,
    225218      integer ecrit_slab_oc !1=ecrit,-1=lit,0=no file
    226219      data ecrit_slab_oc/-1/
    227 !
     220
    228221!     if flag_inhib_forcing = 0, tendencies of forcing are added
    229222!                           <> 0, tendencies of forcing are not added
    … …  
    282275       mth_ini_cas=1 ! pour le moment on compte depuis le debut de l'annee
    283276       call getin('time_ini',heure_ini_cas)
    284 
    285277      print*,'NATURE DE LA SURFACE ',nat_surf
    286 !
     278
    287279! Initialization of the logical switch for nudging
    288280
    … …  
    292284       jcode = jcode/10
    293285     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
    294291!-----------------------------------------------------------------------
    295292!  Definition of the run
    … …  
    398395      allocate(d_q_adv(llm,nqtot))
    399396      allocate(d_q_nudge(llm,nqtot))
    400 !      allocate(d_th_adv(llm))
    401397
    402398      q(:,:) = 0.
    … …  
    406402      d_q_nudge(:,:) = 0.
    407403
    408 !
    409404!   No ozone climatology need be read in this pre-initialization
    410405!          (phys_state_var_init is called again in physiq)
    … …  
    424419!!! Feedback forcing values for Gateaux differentiation (al1)
    425420!!!=====================================================================
    426 !!
     421     
    427422      qsol = qsolinp
    428423      qsurf = fq_sat(tsurf,psurf/100.)
    … …  
    433428      rho(1)=psurf/(rd*tsurf*(1.+(rv/rd-1.)*qsurf))
    434429
    435 !
    436430!! mpl et jyg le 22/08/2012 :
    437431!!  pour que les cas a flux de surface imposes marchent
    … …  
    487481
    488482      INCLUDE "1D_read_forc_cases.h"
    489    print*,'A d_t_adv ',d_t_adv(1:20)*86400
    490483
    491484      if (forcing_GCM2SCM) then
    … …  
    567560        agesno  = xagesno
    568561        tsoil(:,:,:)=tsurf
    569 !------ AMMA 2e run avec modele sol et rayonnement actif (MPL 23052012)
    570 !       tsoil(1,1,1)=299.18
    571 !       tsoil(1,2,1)=300.08
    572 !       tsoil(1,3,1)=301.88
    573 !       tsoil(1,4,1)=305.48
    574 !       tsoil(1,5,1)=308.00
    575 !       tsoil(1,6,1)=308.00
    576 !       tsoil(1,7,1)=308.00
    577 !       tsoil(1,8,1)=308.00
    578 !       tsoil(1,9,1)=308.00
    579 !       tsoil(1,10,1)=308.00
    580 !       tsoil(1,11,1)=308.00
    581562!-----------------------------------------------------------------------
    582563        call pbl_surface_init(fder, snsrf, qsurfsrf, tsoil)
    … …  
    635616          rvc_ancien = 0.
    636617        ENDIF
    637 !jyg<
    638 ! Etienne: comment those lines since now the TKE is inialized in 1D_read_forc_cases
    639 !!      pbl_tke(:,:,:)=1.e-8
    640 !        pbl_tke(:,:,:)=0.
    641 !        pbl_tke(:,2,:)=1.e-2
    642 !>jyg
    643618        rain_fall=0.
    644619        snow_fall=0.
    … …  
    754729
    755730      call phys_state_var_end
    756 !Al1
    757731      if (restart) then
    758732        print*,'call to restart dyn 1d'
    … …  
    762736       print*,'fnday,annee_ref,day_ref,day_ini',                            &
    763737     &     fnday,annee_ref,day_ref,day_ini
    764 !**      call ymds2ju(annee_ref,mois,day_ini,heure,day)
    765738       day = day_ini
    766739       day_end = day_ini + nday
    … …  
    780753! raz for safety
    781754       do l=1,llm
    782          d_q_vert_adv(l,1) = 0.
     755         d_q_vert_adv(l,:) = 0.
    783756       enddo
    784757      endif
    … …  
    819792!---------------------------------------------------------------------
    820793        phis(1)=zsurf*RG
    821 !        phi(1)=phis(1)+RD*temp(1)*(plev(1)-play(1))/(.5*(plev(1)+play(1)))
    822794
    823795        ! Calculate geopotential from the ground surface since phi and phis are added in physiq_mod
    … …  
    842814         z_adv=play
    843815      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
    846822      do l=2,llm-1
    847         ! vertical tendencies computed as d X / d t = -W  d X / d z
    848823        d_u_vert_adv(l)=-w_adv(l)*(u(l+1)-u(l-1))/(z_adv(l+1)-z_adv(l-1))
    849824        d_v_vert_adv(l)=-w_adv(l)*(v(l+1)-v(l-1))/(z_adv(l+1)-z_adv(l-1))
    850825        ! 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)*(teta(l+1)-teta(l-1))/(z_adv(l+1)-z_adv(l-1)) / (pzero/play(l))**rkappa
    852         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))
    853828      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
    858860   
    859861   ENDIF
    … …  
    940942!---------------------------------------------------------------------
    941943!  Nudging
    942 !  EV: rewrite the section to account for a time- and height-varying
    943 !  nudging
    944944!---------------------------------------------------------------------
    945945      d_t_nudge(:) = 0.
    … …  
    10001000
    10011001      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.
    10021028
    10031029!---------------------------------------------------------------------
    … …  
    10221048END IF
    10231049
    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 !!
    10611053
    10621054        teta=temp*(pzero/play)**rkappa
    … …  
    10711063      ENDIF
    10721064
    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
    10901067
    10911068!---------------------------------------------------------------------
    … …  
    11011078        daytime = daytime+1./day_step
    11021079        day = int(daytime+0.1/day_step)
    1103 !        time = max(daytime-day,0.0)
    1104 !Al1&jyg: correction de bug
    1105 !cc        time = real(mod(it,day_step))/day_step
    11061080        time = time_ini/24.+real(mod(it,day_step))/day_step
    11071081        it=it+1

  • LMDZ6/branches/contrails/libf/phylmd/iotd_ecrit.f90

    r5390 r5450  
    123123!Case of a 3D variable
    124124!---------------------
    125         if (llm==lmax) then
     125       if (llm==lmax) then
    126126           ndim=4
    127127           corner(1)=1
    … …  
    134134           edges(4)=1
    135135           dim_cc=dim_coord
    136 
    137136
    138137!Case of a 2D variable
    … …  
    159158      if (ntime==1) then
    160159          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)
    162161          !print*,'DEF ',nom,nid,varid
    163162          ierr = nf90_enddef(nid)

  • LMDZ6/branches/contrails/libf/phylmd/iotd_ini.f90

    r5390 r5450  
    124124      ierr=nf90_def_var(nid, "lev", nf90_float, dim_coord(3),nvarid)
    125125      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
    127128         ierr=nf90_put_att(nid,nvarid,"long_name","pseudo-alt")
    128129         ierr=nf90_put_att(nid,nvarid,'positive',"up")

  • LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp.f90

    r5431 r5450  
    268268  REAL, DIMENSION(klon,klev) :: ctot
    269269  REAL, DIMENSION(klon,klev) :: ctot_vol
    270   REAL, DIMENSION(klon) :: zqs, zdqs
     270  REAL, DIMENSION(klon) :: zqs, zdqs, zqsl, zdqsl, zqsi, zdqsi
    271271  REAL :: zdelta
    272272  REAL, DIMENSION(klon) :: zdqsdT_raw
    … …  
    278278  REAL, DIMENSION(klon) :: qcloud, qincloud_mpc
    279279  REAL, DIMENSION(klon) :: zrfl, zifl
    280   REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn
     280  REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn, zqnl
    281281  REAL, DIMENSION(klon) :: zoliql, zoliqi
    282282  REAL, DIMENSION(klon) :: zt
    283   REAL, DIMENSION(klon) :: zfice,zneb
     283  REAL, DIMENSION(klon) :: zfice,zneb, znebl
    284284  REAL, DIMENSION(klon) :: dzfice
    285285  REAL, DIMENSION(klon) :: zfice_turb, dzfice_turb
    … …  
    306306 
    307307  ! for condensation and ice supersaturation
    308   REAL, DIMENSION(klon) :: qvc, shear
     308  REAL, DIMENSION(klon) :: qvc, qvcl, shear
    309309  REAL :: delta_z
    310310  !--Added for ice supersaturation (ok_ice_supersat) and contrails (ok_plane_contrails)
    … …  
    377377rain(:) = 0.0
    378378snow(:) = 0.0
    379 zfice(:)=0.0
     379zfice(:)=1.0 ! initialized at 1 as by default we assume mpc to be at ice saturation
    380380dzfice(:)=0.0
    381381zfice_turb(:)=0.0
    … …  
    668668                  ! Calculation of saturation specific humidity and ice fraction
    669669                  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
    670683                  CALL calc_gammasat(klon,Tbef,qtot,pplay(:,k),gammasat,dgammasatdt)
    671684                  ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
    672685                  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
    681688                  !--AB Activates a condensation scheme that allows for
    682689                  !--ice supersaturation and contrails evolution from aviation
    … …  
    733740                  ENDIF ! .NOT. ok_ice_supersat
    734741
     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
    735757                  DO i=1,klon
    736758                      IF (keepgoing(i)) THEN
    … …  
    756778                        ENDIF
    757779                       
    758                         ! LEA_R : check formule
    759780                        IF ( ok_unadjusted_clouds ) THEN
    760781                          !--AB We relax the saturation adjustment assumption
    … …  
    796817
    797818
     819        ! Cloud phase final determination
     820        !--------------------------------
    798821        ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
    799822        IF (iflag_t_glace.GE.4) THEN
    … …  
    802825           temp_cltop(:,k)=ztemp_cltop(:)
    803826        ENDIF
    804 
    805         ! Partition function depending on temperature
     827        ! Partition function depending on temperature for all clouds (shallow convective and stratiform)
    806828        CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
    807829
    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
    809831        IF (iflag_icefrac .GE. 1) THEN
    810832           CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), omega(:,k), qice_ini, ziflcld, zqn, &
    … …  
    812834        ENDIF
    813835
    814         ! Water vapor update, Phase determination and subsequent latent heat exchange
     836        ! Water vapor update, subsequent latent heat exchange for each cloud type
     837        !------------------------------------------------------------------------
    815838        DO i=1, klon
    816839            ! Overwrite phase partitioning in boundary layer mixed phase clouds when the
    … …  
    855878                        IF (lognormale(i)) THEN 
    856879                           zcond(i)  = zqliq(i) + zqice(i)
    857                            zfice(i)=zfice_turb(i)
     880                           zfice(i)  = zfice_turb(i)
    858881                           rhcl(i,k) = zqvapcl(i) * rneb(i,k) + (zq(i) - zqn(i)) * (1.-rneb(i,k))
    859882                        ENDIF
    … …  
    865888
    866889            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           
    871892            ! temperature update due to phase change
    872893            zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
    … …  
    896917      zoliql(i) = zcond(i) * ( 1. - zfice(i) )
    897918      zoliqi(i) = zcond(i) * zfice(i)
    898       ! c_iso : initialisation of zoliq* also for isotopes
    899919    ENDDO
    900920
    … …  
    9841004            ENDIF
    9851005
    986             zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
    987             frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
     1006            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
    9881008
    9891009            ! 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))
    9911011
    9921012        ENDIF
    … …  
    10081028                ENDIF
    10091029
    1010               zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
    1011               frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
     1030              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
    10121032
    10131033            ENDIF

  • LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90

    r5406 r5450  
    230230REAL :: dqvc_mix_sub, dqvc_mix_issr
    231231REAL :: dqt_mix
    232 REAL :: a_mix, bovera, N_cld_mix, L_mix
     232REAL :: a_mix, bovera, Povera, N_cld_mix, L_mix
    233233REAL :: envfra_mix, cldfra_mix
    234234REAL :: L_shear, shear_fra
    … …  
    292292  !--formed elsewhere)
    293293  IF (keepgoing(i)) THEN
     294
     295    !--Initialisation
     296    issrfra(i)  = 0.
     297    qissr(i)    = 0.
    294298
    295299    !--If the temperature is higher than the threshold below which
    … …  
    364368      dqvc_mix(i) = 0.
    365369
    366       issrfra(i)  = 0.
    367       qissr(i)    = 0.
    368       subfra(i)   = 0.
    369       qsub(i)     = 0.
    370 
    371370      !--Initialisation of the cell properties
    372371      !--Dry density [kg/m3]
    … …  
    665664        !--Clouds within the mesh are assumed to be ellipses. The length of the
    666665        !--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.
    668668        !--clouds_perim is the total perimeter of the clouds within the mesh,
    669669        !--not considering interfaces with other meshes (only the interfaces with clear
    … …  
    684684        !-- b / a = bovera = MAX(0.1, cldfra)
    685685        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) ) )
    686689        !--The clouds perimeter is imposed using the formula from Morcrette 2012,
    687690        !--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
    706700
    707701        !--The time required for turbulent diffusion to homogenize a region of size

  • LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_ini.f90

    r5412 r5450  
    195195  !$OMP THREADPRIVATE(delta_hetfreez)
    196196 
    197   REAL, SAVE, PROTECTED :: coef_mixing_lscp=9.E-8            ! [-] tuning coefficient for the mixing process
     197  REAL, SAVE, PROTECTED :: coef_mixing_lscp=1.E-3            ! [-] tuning coefficient for the mixing process
    198198  !$OMP THREADPRIVATE(coef_mixing_lscp)
    199199 

  • LMDZ6/branches/contrails/libf/phylmd/lmdz_thermcell_ini.f90

    r5400 r5450  
    88integer, protected :: dvdq=1,dqimpl=-1,prt_level=0,lunout
    99real   , protected :: RG,RD,RCPD,RKAPPA,RLVTT,RLvCp,RETV
     10
     11
    1012
    1113!$OMP THREADPRIVATE(dvdq,dqimpl,prt_level,lunout)
    … …  
    4749!$OMP THREADPRIVATE(iflag_thermals_tenv)
    4850
     51integer, protected       :: thermals_subsid_advect_more_than_one=1
     52character*6, protected :: thermals_subsid_advect_scheme = 'upwind' ! or 'center'
     53
     54!$OMP THREADPRIVATE(thermals_subsid_advect_scheme,thermals_subsid_advect_more_than_one)
    4955
    5056CONTAINS
    … …  
    101107   CALL getin_p('thermals_flag_alim',thermals_flag_alim)
    102108   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)
    104111
    105112
    … …  
    134141write(lunout,*) 'thermcell_ini ,thermals_flag_alim       =',  thermals_flag_alim
    135142write(lunout,*) 'thermcell_ini ,iflag_thermals_tenv      =',  iflag_thermals_tenv
     143write(lunout,*) 'thermcell_ini ,thermals_subsid_advect_scheme=',thermals_subsid_advect_scheme
     144write(lunout,*) 'thermcell_ini ,thermals_subsid_advect_more_than_one=',thermals_subsid_advect_more_than_one
    136145
    137146 RETURN

  • LMDZ6/branches/contrails/libf/phylmd/phys_output_ctrlout_mod.F90

    r5383 r5450  
    20012001  TYPE(ctrl_out), SAVE, ALLOCATABLE :: o_dtr_sat(:)
    20022002  TYPE(ctrl_out), SAVE, ALLOCATABLE :: o_dtr_uscav(:)
     2003  TYPE(ctrl_out), SAVE, ALLOCATABLE :: o_dtr_wet_con(:)
    20032004  TYPE(ctrl_out), SAVE, ALLOCATABLE :: o_dtr_dry(:)
    20042005

  • LMDZ6/branches/contrails/libf/phylmd/phys_output_mod.F90

    r5394 r5450  
    172172    ALLOCATE(o_dtr_evapls(nqtot),o_dtr_ls(nqtot),o_dtr_trsp(nqtot))
    173173    ALLOCATE(o_dtr_sscav(nqtot),o_dtr_sat(nqtot),o_dtr_uscav(nqtot))
     174    ALLOCATE(o_dtr_wet_con(nqtot))
    174175    ALLOCATE(o_dtr_dry(nqtot),o_dtr_vdf(nqtot))
    175176IF (CPPKEY_STRATAER) THEN
    … …  
    540541            tnam = TRIM(dn)//'uscav';       o_dtr_uscav     (itr) = ctrl_out(flag, tnam, lnam, "-", [('',i=1,nfiles)])
    541542
     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)])
    542545            lnam = 'tracer tendency dry deposition'//TRIM(tracers(iq)%longName)
    543546            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  
    66  USE phytrac_mod, ONLY : d_tr_cl, d_tr_th, d_tr_cv, d_tr_lessi_impa, &
    77       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_dry
     8       d_tr_trsp, d_tr_sscav, d_tr_sat, d_tr_uscav, flux_tr_wet, flux_tr_dry
    99
    1010  ! Author: Abderrahmane IDELKADI (original include file)
    … …  
    189189         o_dtr_insc, o_dtr_bcscav, o_dtr_evapls, &
    190190         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, &
    192193         o_ustr_gwd_hines,o_vstr_gwd_hines,o_ustr_gwd_rando,o_vstr_gwd_rando, &
    193194         o_ustr_gwd_front,o_vstr_gwd_front, &
    … …  
    28562857             CALL histwrite_phy(o_dtr_uscav(itr),d_tr_uscav(:,:,itr))
    28572858            !--2D fields
     2859             CALL histwrite_phy(o_dtr_wet_con(itr), flux_tr_wet(:,itr))
    28582860             CALL histwrite_phy(o_dtr_dry(itr), flux_tr_dry(:,itr))
    28592861             zx_tmp_fi2d=0.

  • LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90

    r5431 r5450  
    8080    USE lmdz_call_blowing_snow, ONLY : call_blowing_snow_sublim_sedim
    8181    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
    8284    USE yamada_ini_mod, ONLY : yamada_ini
    8385    USE lmdz_atke_turbulence_ini, ONLY : atke_ini
    … …  
    12631265    CHARACTER(len=512) :: namelist_ecrad_file
    12641266
     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
    12651276    !======================================================================!
    12661277    ! Bifurcation vers un nouveau moniteur physique pour experimenter      !
    … …  
    18181829       CALL iniradia(klon,klev,paprs(1,1:klev+1))
    18191830
     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   
    18201836!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    18211837       CALL wake_ini(rg,rd,rv,prt_level)
    … …  
    38023818    !
    38033819    !===================================================================
     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    !===================================================================
    38043824    ! Computation of ratqs, the width (normalized) of the subrid scale
    38053825    ! water distribution

  • LMDZ6/branches/contrails/libf/phylmd/phytrac_mod.f90

    r5330 r5450  
    3535  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_sat
    3636  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_uscav
     37  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE   :: flux_tr_wet ! tracer wet deposit (surface)                    jyg
    3738  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: qPr,qDi ! concentration tra dans pluie,air descente insaturee
    3839  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: qPa,qMel
    … …  
    4748
    4849!$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)
    5051!$OMP THREADPRIVATE(d_tr_insc,d_tr_bcscav,d_tr_evapls,d_tr_ls,qPrls)
    5152!$OMP THREADPRIVATE(d_tr_cl,d_tr_dry,flux_tr_dry,d_tr_dec,d_tr_cv)
    … …  
    6869    ALLOCATE(d_tr_sscav(klon,klev,nbtr),d_tr_sat(klon,klev,nbtr))
    6970    ALLOCATE(d_tr_uscav(klon,klev,nbtr),qPr(klon,klev,nbtr),qDi(klon,klev,nbtr))
     71    ALLOCATE(flux_tr_wet(klon,nbtr))
    7072    ALLOCATE(qPa(klon,klev,nbtr),qMel(klon,klev,nbtr))
    7173    ALLOCATE(qTrdi(klon,klev,nbtr),dtrcvMA(klon,klev,nbtr))
    … …  
    408410          d_tr_dry(i,it)=0.
    409411          flux_tr_dry(i,it)=0.
     412          flux_tr_wet(i,it)=0.
    410413       ENDDO
    411414    ENDDO
    … …  
    697700                !--with the full array tr_seri even if only item it is processed
    698701
    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,                                    &
    706710                     zmfd1a,zmfphi2,zmfdam)
    707711
Note: See TracChangeset for help on using the changeset viewer.