Changeset 5433


Ignore:
Timestamp:
Dec 19, 2024, 9:57:53 PM (12 hours ago)
Author:
evignon
Message:

modification et correction de la prise en compte de l'advection verticale dans le 1D.
Il y a 3 aspects:

  • on tient compte de l'advection verticale des 3 phases de l'eau
  • on corrige l'application de la tendance qui etait calculee avant l'appel de la physique mais appliquee apres
  • on utilise desormais un schema amont par defaut (le schema centre, mis sous flag, conduisait a des oscillations). C'est un peu fort. Implementation a venir d'un schema van leer

ATTENTION: on perd la convergence sur les cas 1D avec vitesse verticale grande echelle (RICO, SANDU ...)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/dyn1d/scm.f90

    r5302 r5433  
    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
Note: See TracChangeset for help on using the changeset viewer.