Ignore:
Timestamp:
Jan 17, 2015, 2:00:00 PM (10 years ago)
Author:
jyg
Message:

1/ Update of the loop on sub-surfaces in pbl_surface_mod.F90 : all initializations
are removed from the loop and put before the loop.
2/ The possibility of nudging RH, T, u, and v is added in lmdz1d.F90 (and in
1DUTILS.h).

Location:
LMDZ5/trunk/libf/phylmd
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/1DUTILS.h

    r2151 r2181  
    133133        ENDIF
    134134
     135!Config  Key  = iflag_nudge
     136!Config  Desc = atmospheric nudging ttype (decimal code)
     137!Config  Def  = 0
     138!Config  Help = 0 ==> no nudging
     139!  If digit number n of iflag_nudge is set, then nudging of type n is on
     140!  If digit number n of iflag_nudge is not set, then nudging of type n is off
     141!   (digits are numbered from the right)
     142       iflag_nudge = 0
     143       CALL getin('iflag_nudge',iflag_nudge)
     144
    135145!Config  Key  = ok_flux_surf
    136146!Config  Desc = forcage ou non par les flux de surface
    … …  
    39803990!=====================================================================
    39813991
     3992!  Subroutines for nudging
     3993
     3994      Subroutine Nudge_RHT_init (paprs,pplay,t,q,t_targ,rh_targ)
     3995! ========================================================
     3996  USE dimphy
     3997
     3998  implicit none
     3999
     4000! ========================================================
     4001      REAL paprs(klon,klevp1)
     4002      REAL pplay(klon,klev)
     4003!
     4004!      Variables d'etat
     4005      REAL t(klon,klev)
     4006      REAL q(klon,klev)
     4007!
     4008!   Profiles cible
     4009      REAL t_targ(klon,klev)
     4010      REAL rh_targ(klon,klev)
     4011!
     4012   INTEGER k,i
     4013   REAL zx_qs
     4014
     4015! Declaration des constantes et des fonctions thermodynamiques
     4016!
     4017include "YOMCST.h"
     4018include "YOETHF.h"
     4019!
     4020!  ----------------------------------------
     4021!  Statement functions
     4022include "FCTTRE.h"
     4023!  ----------------------------------------
     4024!
     4025        DO k = 1,klev
     4026         DO i = 1,klon
     4027           t_targ(i,k) = t(i,k)
     4028           IF (t(i,k).LT.RTT) THEN
     4029              zx_qs = qsats(t(i,k))/(pplay(i,k))
     4030           ELSE
     4031              zx_qs = qsatl(t(i,k))/(pplay(i,k))
     4032           ENDIF
     4033           rh_targ(i,k) = q(i,k)/zx_qs
     4034         ENDDO
     4035        ENDDO
     4036      print *, 't_targ',t_targ
     4037      print *, 'rh_targ',rh_targ
     4038!
     4039!
     4040      RETURN
     4041      END
     4042
     4043      Subroutine Nudge_UV_init (paprs,pplay,u,v,u_targ,v_targ)
     4044! ========================================================
     4045  USE dimphy
     4046
     4047  implicit none
     4048
     4049! ========================================================
     4050      REAL paprs(klon,klevp1)
     4051      REAL pplay(klon,klev)
     4052!
     4053!      Variables d'etat
     4054      REAL u(klon,klev)
     4055      REAL v(klon,klev)
     4056!
     4057!   Profiles cible
     4058      REAL u_targ(klon,klev)
     4059      REAL v_targ(klon,klev)
     4060!
     4061   INTEGER k,i
     4062!
     4063        DO k = 1,klev
     4064         DO i = 1,klon
     4065           u_targ(i,k) = u(i,k)
     4066           v_targ(i,k) = v(i,k)
     4067         ENDDO
     4068        ENDDO
     4069      print *, 'u_targ',u_targ
     4070      print *, 'v_targ',v_targ
     4071!
     4072!
     4073      RETURN
     4074      END
     4075
     4076      Subroutine Nudge_RHT (dtime,paprs,pplay,t_targ,rh_targ,t,q,          &
     4077     &                      d_t,d_q)
     4078! ========================================================
     4079  USE dimphy
     4080
     4081  implicit none
     4082
     4083! ========================================================
     4084      REAL dtime
     4085      REAL paprs(klon,klevp1)
     4086      REAL pplay(klon,klev)
     4087!
     4088!      Variables d'etat
     4089      REAL t(klon,klev)
     4090      REAL q(klon,klev)
     4091!
     4092! Tendances
     4093      REAL d_t(klon,klev)
     4094      REAL d_q(klon,klev)
     4095!
     4096!   Profiles cible
     4097      REAL t_targ(klon,klev)
     4098      REAL rh_targ(klon,klev)
     4099!
     4100!   Temps de relaxation
     4101      REAL tau
     4102!c      DATA tau /3600./
     4103!!      DATA tau /5400./
     4104      DATA tau /1800./
     4105!
     4106   INTEGER k,i
     4107   REAL zx_qs, rh, tnew, d_rh
     4108
     4109! Declaration des constantes et des fonctions thermodynamiques
     4110!
     4111include "YOMCST.h"
     4112include "YOETHF.h"
     4113!
     4114!  ----------------------------------------
     4115!  Statement functions
     4116include "FCTTRE.h"
     4117!  ----------------------------------------
     4118!
     4119        print *,'dtime, tau ',dtime,tau
     4120        print *, 't_targ',t_targ
     4121        print *, 'rh_targ',rh_targ
     4122        print *,'temp ',t
     4123        print *,'hum ',q
     4124        DO k = 1,klev
     4125         DO i = 1,klon
     4126!!           IF (paprs(i,1)-pplay(i,k) .GT. 10000.) THEN
     4127            IF (t(i,k).LT.RTT) THEN
     4128               zx_qs = qsats(t(i,k))/(pplay(i,k))
     4129            ELSE
     4130               zx_qs = qsatl(t(i,k))/(pplay(i,k))
     4131            ENDIF
     4132            rh = q(i,k)/zx_qs
     4133!
     4134            d_t(i,k) = d_t(i,k) + 1./tau*(t_targ(i,k)-t(i,k))
     4135            d_rh = 1./tau*(rh_targ(i,k)-rh)
     4136!
     4137            tnew = t(i,k)+d_t(i,k)
     4138            IF (tnew.LT.RTT) THEN
     4139               zx_qs = qsats(tnew)/(pplay(i,k))
     4140            ELSE
     4141               zx_qs = qsatl(tnew)/(pplay(i,k))
     4142            ENDIF
     4143            d_q(i,k) = d_q(i,k) + d_rh*zx_qs
     4144!
     4145            print *,' k,d_t,rh,d_rh,d_q ',    &
     4146                      k,d_t(i,k),rh,d_rh,d_q(i,k)
     4147!!           ENDIF
     4148!
     4149         ENDDO
     4150        ENDDO
     4151!
     4152      RETURN
     4153      END
     4154
     4155      Subroutine Nudge_UV (dtime,paprs,pplay,u_targ,v_targ,u,v,          &
     4156     &                      d_u,d_v)
     4157! ========================================================
     4158  USE dimphy
     4159
     4160  implicit none
     4161
     4162! ========================================================
     4163      REAL dtime
     4164      REAL paprs(klon,klevp1)
     4165      REAL pplay(klon,klev)
     4166!
     4167!      Variables d'etat
     4168      REAL u(klon,klev)
     4169      REAL v(klon,klev)
     4170!
     4171! Tendances
     4172      REAL d_u(klon,klev)
     4173      REAL d_v(klon,klev)
     4174!
     4175!   Profiles cible
     4176      REAL u_targ(klon,klev)
     4177      REAL v_targ(klon,klev)
     4178!
     4179!   Temps de relaxation
     4180      REAL tau
     4181!c      DATA tau /3600./
     4182      DATA tau /5400./
     4183!
     4184   INTEGER k,i
     4185
     4186!
     4187        print *,'dtime, tau ',dtime,tau
     4188        print *, 'u_targ',u_targ
     4189        print *, 'v_targ',v_targ
     4190        print *,'zonal velocity ',u
     4191        print *,'meridional velocity ',v
     4192        DO k = 1,klev
     4193         DO i = 1,klon
     4194           IF (paprs(i,1)-pplay(i,k) .GT. 10000.) THEN
     4195!
     4196            d_u(i,k) = d_u(i,k) + 1./tau*(u_targ(i,k)-u(i,k))
     4197            d_v(i,k) = d_v(i,k) + 1./tau*(v_targ(i,k)-v(i,k))
     4198!
     4199            print *,' k,u,d_u,v,d_v ',    &
     4200                      k,u(i,k),d_u(i,k),v(i,k),d_v(i,k)
     4201           ENDIF
     4202!
     4203         ENDDO
     4204        ENDDO
     4205!
     4206      RETURN
     4207      END
     4208

  • LMDZ5/trunk/libf/phylmd/change_srf_frac_mod.F90

    r1907 r2181  
    4949    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: u10m
    5050    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: v10m
    51     REAL, DIMENSION(klon,klev+1,nbsrf), INTENT(INOUT) :: pbl_tke
     51!jyg<
     52!!    REAL, DIMENSION(klon,klev+1,nbsrf), INTENT(INOUT) :: pbl_tke
     53    REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(INOUT) :: pbl_tke
     54!>jyg
    5255
    5356! Loccal variables

  • LMDZ5/trunk/libf/phylmd/compar1d.h

    r2019 r2181  
    33!
    44      integer :: forcing_type
     5      integer :: iflag_nudge
    56      real :: nat_surf
    67      real :: tsurf
    … …  
    3233     & wtsurf,wqsurf,restart_runoff,xagesno,qsolinp,zpicinp,            &
    3334     & forcing_type,                                                    &
     35     & iflag_nudge,                                                     &
    3436     & restart,ok_old_disvert
    3537

  • LMDZ5/trunk/libf/phylmd/lmdz1d.F90

    r2126 r2181  
    138138!
    139139!---------------------------------------------------------------------
     140!  Declarations related to nudging
     141!---------------------------------------------------------------------
     142     integer :: nudge_max
     143     parameter (nudge_max=9)
     144     integer :: inudge_RHT=1
     145     integer :: inudge_UV=2
     146     logical :: nudge(nudge_max)
     147     real :: t_targ(llm)
     148     real :: rh_targ(llm)
     149     real :: u_targ(llm)
     150     real :: v_targ(llm)
     151!
     152!---------------------------------------------------------------------
    140153!  Declarations related to vertical discretization:
    141154!---------------------------------------------------------------------
    … …  
    156169      real :: du_phys(llm),dv_phys(llm),dt_phys(llm)
    157170      real :: dt_dyn(llm)
    158       real :: dt_cooling(llm),d_th_adv(llm)
     171      real :: dt_cooling(llm),d_th_adv(llm),d_t_nudge(llm)
     172      real :: d_u_nudge(llm),d_v_nudge(llm)
    159173      real :: alpha
    160174      real :: ttt
    … …  
    164178      REAL, ALLOCATABLE, DIMENSION(:,:):: dq_dyn
    165179      REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_adv
    166 !     REAL, ALLOCATABLE, DIMENSION(:):: d_th_adv
     180      REAL, ALLOCATABLE, DIMENSION(:,:):: d_q_nudge
     181!      REAL, ALLOCATABLE, DIMENSION(:):: d_th_adv
    167182
    168183!---------------------------------------------------------------------
    … …  
    211226!---------------------------------------------------------------------
    212227      integer :: k,l,i,it=1,mxcalc
     228      integer jcode
    213229      integer jjmp1
    214230      parameter (jjmp1=jjm+1-1/jjm)
    … …  
    330346        if (forcing_toga.or.forcing_sandu.or.forcing_astex .or. forcing_dice)                 &
    331347     &    type_ts_forcing = 1
    332 
     348!
     349! Initialization of the logical switch for nudging
     350     jcode = iflag_nudge
     351     do i = 1,nudge_max
     352       nudge(i) = mod(jcode,10) .ge. 1
     353       jcode = jcode/10
     354     enddo
    333355!---------------------------------------------------------------------
    334356!  Definition of the run
    … …  
    444466      allocate(dq_dyn(llm,nqtot))
    445467      allocate(d_q_adv(llm,nqtot))
    446 !     allocate(d_th_adv(llm))
     468      allocate(d_q_nudge(llm,nqtot))
     469!      allocate(d_th_adv(llm))
    447470
    448471!
    … …  
    751774         open(97,file='div_slab.dat',STATUS='OLD')
    752775       endif
     776!
     777!---------------------------------------------------------------------
     778!    Initialize target profile for RHT nudging if needed
     779!---------------------------------------------------------------------
     780      if (nudge(inudge_RHT)) then
     781        call nudge_RHT_init(plev,play,temp,q(:,1),t_targ,rh_targ)
     782      endif
     783      if (nudge(inudge_UV)) then
     784        call nudge_UV_init(plev,play,u,v,u_targ,v_targ)
     785      endif
     786!
    753787!=====================================================================
    754788! START OF THE TEMPORAL LOOP :
    … …  
    876910!!     : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
    877911!
     912!!!!!!!!!!!!!!!!!!!!!!!!
     913!  Nudging
     914!!!!!!!!!!!!!!!!!!!!!!!!
     915      d_t_nudge(:) = 0.
     916      d_q_nudge(:,:) = 0.
     917      d_u_nudge(:) = 0.
     918      d_v_nudge(:) = 0.
     919      if (nudge(inudge_RHT)) then
     920        call nudge_RHT(timestep,plev,play,t_targ,rh_targ,temp,q(:,1),     &
     921    &                  d_t_nudge,d_q_nudge(:,1))
     922      endif
     923      if (nudge(inudge_UV)) then
     924        call nudge_UV(timestep,plev,play,u_targ,v_targ,u,v,     &
     925    &                  d_u_nudge,d_v_nudge)
     926      endif
     927!
    878928!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    879929!         call  writefield_phy('dv_age' ,dv_age,llm)
    … …  
    893943        u(1:mxcalc)=u(1:mxcalc) + timestep*(                                &
    894944     &              du_phys(1:mxcalc)                                       &
    895      &             +du_age(1:mxcalc) )           
     945     &             +du_age(1:mxcalc)                                        &
     946     &             +d_u_nudge(1:mxcalc) )           
    896947        v(1:mxcalc)=v(1:mxcalc) + timestep*(                                 &
    897948     &              dv_phys(1:mxcalc)                                       &
    898      &             +dv_age(1:mxcalc) )
     949     &             +dv_age(1:mxcalc)                                        &
     950     &             +d_v_nudge(1:mxcalc) )
    899951        q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*(                              &
    900952     &                dq(1:mxcalc,:)                                        &
    901      &               +d_q_adv(1:mxcalc,:) )
     953     &               +d_q_adv(1:mxcalc,:)                                   &
     954     &               +d_q_nudge(1:mxcalc,:) )
    902955
    903956        if (prt_level.ge.1) then
    … …  
    913966     &              dt_phys(1:mxcalc)                                       &
    914967     &             +d_th_adv(1:mxcalc)                                      &
     968     &             +d_t_nudge(1:mxcalc)                                      &
    915969     &             +dt_cooling(1:mxcalc))  ! Taux de chauffage ou refroid.
    916970

  • LMDZ5/trunk/libf/phylmd/pbl_surface_mod.F90

    r2177 r2181  
    1212  USE dimphy
    1313  USE mod_phys_lmdz_para,  ONLY : mpi_size
     14  USE mod_grid_phy_lmdz,   ONLY : klon_glo
    1415  USE ioipsl
    1516  USE surface_data,        ONLY : type_ocean, ok_veget
    … …  
    210211       wfbils,    wfbilo,    flux_t,   flux_u, flux_v,&
    211212       dflux_t,   dflux_q,   zxsnow,                  &
    212        zxfluxt,   zxfluxq,   q2m,      flux_q, tke,   &
     213!jyg<
     214!!       zxfluxt,   zxfluxq,   q2m,      flux_q, tke,   &
     215       zxfluxt,   zxfluxq,   q2m,      flux_q, tke_x,   &
     216!>jyg
    213217!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
    214218!!        tke_x,     tke_w                              &
    … …  
    267271! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
    268272!                    (orientation positive vers le bas)
    269 ! tke---input/output-R- tke (kg/m**2/s)
     273! tke_x---input/output-R- tke in the (x) region (kg/m**2/s)
    270274! wake_dltke-input/output-R- tke difference between (w) and (x) (kg/m**2/s)
    271275! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
    … …  
    347351    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: alb1    ! albedo in visible SW interval
    348352    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: alb2    ! albedo in near infra-red SW interval
     353!jyg Pourquoi ustar et wstar sont-elles INOUT ?
    349354    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: ustar   ! u* (m/s)
    350355    REAL, DIMENSION(klon, nbsrf+1), INTENT(INOUT)   :: wstar   ! w* (m/s)
    351356    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: u10m    ! u speed at 10m
    352357    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: v10m    ! v speed at 10m
    353     REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke
     358!jyg<
     359!!    REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke
     360    REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke_x
     361!>jyg
    354362
    355363!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
    356     REAL, DIMENSION(klon, klev+1, nbsrf), INTENT(INOUT) :: wake_dltke ! TKE_w - TKE_x
     364    REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: wake_dltke ! TKE_w - TKE_x
    357365!!!
    358366
    … …  
    604612!!!
    605613!!!jyg le 08/02/2012
     614    REAL, DIMENSION(klon, nbsrf)       :: windsp
     615!
    606616    REAL, DIMENSION(klon, nbsrf)       :: t2m_x
    607617    REAL, DIMENSION(klon, nbsrf)       :: q2m_x
    … …  
    797807     
    798808       ! Initialize ok_flux_surf (for 1D model)
    799        if (klon>1) ok_flux_surf=.FALSE.
     809       if (klon_glo>1) ok_flux_surf=.FALSE.
    800810       
    801811       ! Initilize debug IO
    … …  
    841851!****************************************************************************************
    842852! 2) Initialization to zero
    843 !    Done for all local variables that will be compressed later
    844 !    and argument with INTENT(OUT)
    845 !****************************************************************************************
    846     cdragh = 0.0  ; cdragm = 0.0     ; dflux_t = 0.0   ; dflux_q = 0.0
    847     ypct = 0.0    ; yts = 0.0        ; ysnow = 0.0
    848     zv1 = 0.0     ; yqsurf = 0.0     ; yalb1 = 0.0     ; yalb2 = 0.0   
    849     yrain_f = 0.0 ; ysnow_f = 0.0    ; yfder = 0.0     ; ysolsw = 0.0   
    850     ysollw = 0.0  ; yrugos = 0.0     ; yu1 = 0.0   
    851     yv1 = 0.0     ; ypaprs = 0.0     ; ypplay = 0.0
    852     ydelp = 0.0   ; yu = 0.0         ; yv = 0.0        ; yt = 0.0         
    853     yq = 0.0      ; y_dflux_t = 0.0  ; y_dflux_q = 0.0
    854     yrugoro = 0.0 ; ywindsp = 0.0   
    855     d_ts = 0.0    ; yfluxlat=0.0     ; flux_t = 0.0    ; flux_q = 0.0     
    856     flux_u = 0.0  ; flux_v = 0.0     ; d_t = 0.0       ; d_q = 0.0     
    857     d_t_diss= 0.0 ;d_u = 0.0     ; d_v = 0.0        ; yqsol = 0.0   
    858     ytherm = 0.0  ; ytke=0.
    859     ! Martin
    860     ysnowhgt = 0.0; yqsnow = 0.0     ; yrunoff = 0.0   ; ytoice =0.0
    861     yalb3_new = 0.0  ; ysissnow = 0.0  ; ysollwd = 0.0
    862     ypphi = 0.0   ; ycldt = 0.0      ; yrmu0 = 0.0
    863     ! Martin
    864 
    865 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
    866     ytke_x=0.     ; ytke_w=0.        ; ywake_dltke=0.
    867     y_d_t_x=0.    ; y_d_t_w=0.       ; y_d_q_x=0.      ; y_d_q_w=0.
    868     d_t_w=0.      ; d_q_w=0.         
    869     d_t_x=0.      ; d_q_x=0.
    870     d_t_diss_x = 0. ; d_t_diss_w = 0.
    871 !!    d_wake_dlt=0.    ; d_wake_dlq=0.
    872     d_u_x=0.      ; d_u_w=0.         ; d_v_x=0.        ; d_v_w=0.
    873     flux_t_x=0.   ; flux_t_w=0.      ; flux_q_x=0.     ; flux_q_w=0.
    874     yfluxlat_x=0. ; yfluxlat_w=0.
    875     ywake_s=0.    ; ywake_cstar=0.   ;ywake_dens=0.
    876 !!!
    877 !!! nrlmd le 13/06/2011
    878     tau_eq=0.     ; delta_coef=0.
    879     y_delta_flux_t1=0.
    880     ydtsurf_th=0.
    881     yts_x=0.      ; yts_w=0.
    882     y_delta_tsurf=0.
    883     cdragh_x=0.   ; cdragh_w=0.      ; cdragm_x=0.     ;cdragm_w=0.
    884     kh=0.         ; kh_x=0.          ; kh_w=0.
    885 !!!
    886     tke(:,:,is_ave)=0.
     853!****************************************************************************************
     854!
     855! 2a) Initialization of all argument variables with INTENT(OUT)
     856!****************************************************************************************
     857 lwdown_m(:)=0.
     858 cdragh(:)=0. ; cdragm(:)=0.
     859 zu1(:)=0. ; zv1(:)=0.
     860 alb1_m(:)=0. ; alb2_m(:)=0. ; alb3_lic(:)=0.
     861 zxsens(:)=0. ; zxevap(:)=0. ; zxtsol(:)=0.
     862 d_t_w(:,:)=0. ; d_q_w(:,:)=0. ; d_t_x(:,:)=0. ; d_q_x(:,:)=0.
     863 zxfluxlat(:)=0.
     864 zt2m(:)=0. ; zq2m(:)=0. ; qsat2m(:)=0. ; rh2m(:)=0.
     865 d_t(:,:)=0. ; d_t_diss(:,:)=0. ; d_q(:,:)=0. ; d_u(:,:)=0. ; d_v(:,:)=0.
     866 zcoefh(:,:,:)=0. ; zcoefm(:,:,:)=0.
     867 zxsens_x(:)=0. ; zxsens_w(:)=0. ; zxfluxlat_x(:)=0. ; zxfluxlat_w(:)=0.
     868 cdragh_x(:)=0. ; cdragh_w(:)=0. ; cdragm_x(:)=0. ; cdragm_w(:)=0.
     869 kh(:)=0. ; kh_x(:)=0. ; kh_w(:)=0.
     870 slab_wfbils(:)=0.
     871 qsol_d(:)=0.
     872 s_pblh(:)=0. ; s_pblh_x(:)=0. ; s_pblh_w(:)=0.
     873 s_plcl(:)=0. ; s_plcl_x(:)=0. ; s_plcl_w(:)=0.
     874 s_capCL(:)=0. ; s_oliqCL(:)=0. ; s_cteiCL(:)=0. ; s_pblT(:)=0.
     875 s_therm(:)=0.
     876 s_trmb1(:)=0. ; s_trmb2(:)=0. ; s_trmb3(:)=0.
     877 zxrugs(:)=0. ; zustar(:)=0.
     878 zu10m(:)=0. ; zv10m(:)=0.
     879 fder_print(:)=0.
     880 zxqsurf(:)=0.
     881 zxfluxu(:,:)=0. ; zxfluxv(:,:)=0.
     882 rugos_d(:,:)=0. ; agesno_d(:,:)=0.
     883 solsw(:,:)=0. ; sollw(:,:)=0.
     884 d_ts(:,:)=0.
     885 evap_d(:,:)=0.
     886 fluxlat(:,:)=0.
     887 wfbils(:,:)=0. ; wfbilo(:,:)=0.
     888 flux_t(:,:,:)=0. ; flux_q(:,:,:)=0. ; flux_u(:,:,:)=0. ; flux_v(:,:,:)=0.
     889 dflux_t(:)=0. ; dflux_q(:)=0.
     890 zxsnow(:)=0.
     891 zxfluxt(:,:)=0. ; zxfluxq(:,:)=0.
     892 qsnow(:)=0. ; snowhgt(:)=0. ; to_ice(:)=0. ; sissnow(:)=0.
     893 runoff(:)=0.
    887894    IF (iflag_pbl<20.or.iflag_pbl>=30) THEN
    888895       zcoefh(:,:,:) = 0.0
    … …  
    894901      zcoefh(:,:,is_ave)=0.
    895902    ENDIF
    896     ytsoil = 999999.
    897 
     903!!
     904!  The components "is_ave" of tke_x and wake_deltke are "OUT" variables
     905!jyg<
     906!!    tke(:,:,is_ave)=0.
     907    tke_x(:,:,is_ave)=0.
     908    wake_dltke(:,:,is_ave)=0.
     909!>jyg
    898910!!! jyg le 23/02/2013
    899     pblh(:,:)      = 999999.     ! pblh,plcl,cteiCL are meaningfull only over sub-surfaces
    900     plcl(:,:)      = 999999.     ! actually present in the grid cell.
    901     cteiCL(:,:)    = 999999.
    902     pblh_x(:,:)      = 999999. 
    903     plcl_x(:,:)      = 999999. 
    904     cteiCL_x(:,:)    = 999999. 
    905     pblh_w(:,:)      = 999999.   
    906     plcl_w(:,:)      = 999999.   
    907     cteiCL_w(:,:)    = 999999.   
    908 !
    909911    t2m(:,:)       = 999999.     ! t2m and q2m are meaningfull only over sub-surfaces
    910912    q2m(:,:)       = 999999.     ! actually present in the grid cell.
    911913!!!
    912     rh2m(:)        = 0.
    913     qsat2m(:)      = 0.
     914    rh2m(:) = 0. ; qsat2m(:) = 0.
    914915!!!
    915916!!! jyg le 10/02/2012
    916     rh2m_x(:)        = 0.
    917     qsat2m_x(:)      = 0.
    918     rh2m_w(:)        = 0.
    919     qsat2m_w(:)      = 0.
     917    rh2m_x(:) = 0. ; qsat2m_x(:) = 0. ; rh2m_w(:) = 0. ; qsat2m_w(:) = 0.
     918!!!
     919
     920! 2b) Initialization of all local variables that will be compressed later
     921!****************************************************************************************
     922!!    cdragh = 0.0  ; cdragm = 0.0     ; dflux_t = 0.0   ; dflux_q = 0.0
     923    ypct = 0.0    ; yts = 0.0        ; ysnow = 0.0
     924!!    zv1 = 0.0     ; yqsurf = 0.0     ; yalb1 = 0.0     ; yalb2 = 0.0   
     925    yqsurf = 0.0  ; yalb1 = 0.0      ; yalb2 = 0.0   
     926    yrain_f = 0.0 ; ysnow_f = 0.0    ; yfder = 0.0     ; ysolsw = 0.0   
     927    ysollw = 0.0  ; yrugos = 0.0     ; yu1 = 0.0   
     928    yv1 = 0.0     ; ypaprs = 0.0     ; ypplay = 0.0
     929    ydelp = 0.0   ; yu = 0.0         ; yv = 0.0        ; yt = 0.0         
     930    yq = 0.0      ; y_dflux_t = 0.0  ; y_dflux_q = 0.0
     931    yrugoro = 0.0 ; ywindsp = 0.0   
     932!!    d_ts = 0.0    ; yfluxlat=0.0     ; flux_t = 0.0    ; flux_q = 0.0     
     933    yfluxlat=0.0
     934!!    flux_u = 0.0  ; flux_v = 0.0     ; d_t = 0.0       ; d_q = 0.0     
     935!!    d_t_diss= 0.0 ;d_u = 0.0     ; d_v = 0.0
     936    yqsol = 0.0   
     937    ytherm = 0.0  ; ytke=0.
     938    ! Martin
     939    ysnowhgt = 0.0; yqsnow = 0.0     ; yrunoff = 0.0   ; ytoice =0.0
     940    yalb3_new = 0.0  ; ysissnow = 0.0  ; ysollwd = 0.0
     941    ypphi = 0.0   ; ycldt = 0.0      ; yrmu0 = 0.0
     942    ! Martin
     943
     944!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
     945    ytke_x=0.     ; ytke_w=0.        ; ywake_dltke=0.
     946    y_d_t_x=0.    ; y_d_t_w=0.       ; y_d_q_x=0.      ; y_d_q_w=0.
     947!!    d_t_w=0.      ; d_q_w=0.         
     948!!    d_t_x=0.      ; d_q_x=0.
     949!!    d_wake_dlt=0.    ; d_wake_dlq=0.
     950    yfluxlat_x=0. ; yfluxlat_w=0.
     951    ywake_s=0.    ; ywake_cstar=0.   ;ywake_dens=0.
     952!!!
     953!!! nrlmd le 13/06/2011
     954    tau_eq=0.     ; delta_coef=0.
     955    y_delta_flux_t1=0.
     956    ydtsurf_th=0.
     957    yts_x=0.      ; yts_w=0.
     958    y_delta_tsurf=0.
     959!!!
     960    ytsoil = 999999.
     961
     962
     963! 2c) Initialization of all local variables computed within the subsurface loop and used later on
     964!****************************************************************************************
     965    d_t_diss_x(:,:) = 0. ;        d_t_diss_w(:,:) = 0.
     966    d_u_x(:,:)=0. ;               d_u_w(:,:)=0.
     967    d_v_x(:,:)=0. ;               d_v_w(:,:)=0.
     968    flux_t_x(:,:,:)=0. ;          flux_t_w(:,:,:)=0.
     969    flux_q_x(:,:,:)=0. ;          flux_q_w(:,:,:)=0.
     970!
     971!jyg<
     972    flux_u_x(:,:,:)=0. ;          flux_u_w(:,:,:)=0.
     973    flux_v_x(:,:,:)=0. ;          flux_v_w(:,:,:)=0.
     974    fluxlat_x(:,:)=0. ;           fluxlat_w(:,:)=0.
     975!>jyg
     976!
     977!jyg<
     978! pblh,plcl,capCL,cteiCL ... are meaningfull only over sub-surfaces
     979! actually present in the grid cell  ==> value set to 999999.
     980!                           
     981!jyg<
     982       ustar(:,:)   = 999999.
     983       wstar(:,:)   = 999999.
     984       windsp(:,:)  = SQRT(u10m(:,:)**2 + v10m(:,:)**2 )
     985       u10m(:,:)    = 999999.
     986       v10m(:,:)    = 999999.
     987!>jyg
     988!
     989       pblh(:,:)   = 999999.        ! Hauteur de couche limite
     990       plcl(:,:)   = 999999.        ! Niveau de condensation de la CLA
     991       capCL(:,:)  = 999999.        ! CAPE de couche limite
     992       oliqCL(:,:) = 999999.        ! eau_liqu integree de couche limite
     993       cteiCL(:,:) = 999999.        ! cloud top instab. crit. couche limite
     994       pblt(:,:)   = 999999.        ! T a la Hauteur de couche limite
     995       therm(:,:)  = 999999.
     996       trmb1(:,:)  = 999999.        ! deep_cape
     997       trmb2(:,:)  = 999999.        ! inhibition
     998       trmb3(:,:)  = 999999.        ! Point Omega
     999!
     1000       t2m_x(:,:)    = 999999.
     1001       q2m_x(:,:)    = 999999.
     1002       ustar_x(:,:)   = 999999.
     1003       wstar_x(:,:)   = 999999.
     1004       u10m_x(:,:)   = 999999.
     1005       v10m_x(:,:)   = 999999.
     1006!                           
     1007       pblh_x(:,:)   = 999999.      ! Hauteur de couche limite
     1008       plcl_x(:,:)   = 999999.      ! Niveau de condensation de la CLA
     1009       capCL_x(:,:)  = 999999.      ! CAPE de couche limite
     1010       oliqCL_x(:,:) = 999999.      ! eau_liqu integree de couche limite
     1011       cteiCL_x(:,:) = 999999.      ! cloud top instab. crit. couche limite
     1012       pblt_x(:,:)   = 999999.      ! T a la Hauteur de couche limite
     1013       therm_x(:,:)  = 999999.     
     1014       trmb1_x(:,:)  = 999999.      ! deep_cape
     1015       trmb2_x(:,:)  = 999999.      ! inhibition
     1016       trmb3_x(:,:)  = 999999.      ! Point Omega
     1017!
     1018       t2m_w(:,:)    = 999999.
     1019       q2m_w(:,:)    = 999999.
     1020       ustar_w(:,:)   = 999999.
     1021       wstar_w(:,:)   = 999999.
     1022       u10m_w(:,:)   = 999999.
     1023       v10m_w(:,:)   = 999999.
     1024                           
     1025       pblh_w(:,:)   = 999999.      ! Hauteur de couche limite
     1026       plcl_w(:,:)   = 999999.      ! Niveau de condensation de la CLA
     1027       capCL_w(:,:)  = 999999.      ! CAPE de couche limite
     1028       oliqCL_w(:,:) = 999999.      ! eau_liqu integree de couche limite
     1029       cteiCL_w(:,:) = 999999.      ! cloud top instab. crit. couche limite
     1030       pblt_w(:,:)   = 999999.      ! T a la Hauteur de couche limite
     1031       therm_w(:,:)  = 999999.     
     1032       trmb1_w(:,:)  = 999999.      ! deep_cape
     1033       trmb2_w(:,:)  = 999999.      ! inhibition
     1034       trmb3_w(:,:)  = 999999.      ! Point Omega
     1035!!!     
     1036!
    9201037!!!
    9211038!****************************************************************************************
    … …  
    10701187          yv1(j)     = v(i,1)
    10711188          ypaprs(j,klev+1) = paprs(i,klev+1)
    1072           ywindsp(j) = SQRT(u10m(i,nsrf)**2 + v10m(i,nsrf)**2 )
     1189!jyg<
     1190!!          ywindsp(j) = SQRT(u10m(i,nsrf)**2 + v10m(i,nsrf)**2 )
     1191          ywindsp(j) = windsp(i,nsrf)
     1192!>jyg
    10731193          ! Martin
    10741194          yzsig(j)   = zsig(i)
    … …  
    10931213          DO j = 1, knon
    10941214             i = ni(j)
    1095              ytke(j,k)   = tke(i,k,nsrf)
     1215!jyg<
     1216!!             ytke(j,k)   = tke(i,k,nsrf)
     1217             ytke(j,k)   = tke_x(i,k,nsrf)
     1218!>jyg
    10961219             yu(j,k) = u(i,k)
    10971220             yv(j,k) = v(i,k)
    … …  
    11211244          DO j = 1, knon
    11221245             i = ni(j)
    1123              ytke_x(j,k) = tke(i,k,nsrf)-wake_s(i)*wake_dltke(i,k,nsrf)
    1124              ytke_w(j,k) = tke(i,k,nsrf)+(1.-wake_s(i))*wake_dltke(i,k,nsrf)
     1246!jyg<
     1247!!             ytke_x(j,k) = tke(i,k,nsrf)-wake_s(i)*wake_dltke(i,k,nsrf)
     1248!!             ytke_w(j,k) = tke(i,k,nsrf)+(1.-wake_s(i))*wake_dltke(i,k,nsrf)
     1249!!             ywake_dltke(j,k) = wake_dltke(i,k,nsrf)
     1250!!             ytke(j,k)     = tke(i,k,nsrf)
     1251!
     1252             ytke_x(j,k)      = tke_x(i,k,nsrf)
     1253             ytke(j,k)        = tke_x(i,k,nsrf)+wake_s(i)*wake_dltke(i,k,nsrf)
     1254             ytke_w(j,k)      = tke_x(i,k,nsrf)+wake_dltke(i,k,nsrf)
    11251255             ywake_dltke(j,k) = wake_dltke(i,k,nsrf)
    1126              ytke(j,k)     = tke(i,k,nsrf)
     1256!>jyg
    11271257          ENDDO
    11281258        ENDDO
    … …  
    16401770          !     y_flux_u1, y_flux_v1)
    16411771
    1642           alb3_lic(:)=0.
     1772!jyg<
     1773!!          alb3_lic(:)=0.
     1774!>jyg
    16431775          DO j = 1, knon
    16441776             i = ni(j)
    … …  
    20532185!      print*,'Dans pbl OK1'
    20542186
    2055        evap(:,nsrf) = - flux_q(:,1,nsrf)
    2056        
    2057        alb1(:, nsrf) = 0.
    2058        alb2(:, nsrf) = 0.
    2059        snow(:, nsrf) = 0.
    2060        qsurf(:, nsrf) = 0.
    2061        rugos(:, nsrf) = 0.
    2062        fluxlat(:,nsrf) = 0.
     2187!jyg<
     2188!!       evap(:,nsrf) = - flux_q(:,1,nsrf)
     2189!>jyg
    20632190       DO j = 1, knon
    20642191          i = ni(j)
     2192          evap(i,nsrf) = - flux_q(i,1,nsrf)                  !jyg
    20652193          d_ts(i,nsrf) = y_d_ts(j)
    20662194          alb1(i,nsrf) = yalb1_new(j) 
    … …  
    20832211!!!
    20842212!!! nrlmd le 02/05/2011
    2085         fluxlat_x(:,nsrf) = 0.
    2086         fluxlat_w(:,nsrf) = 0.
    20872213        DO j = 1, knon
    20882214          i = ni(j)
    … …  
    21282254           DO j = 1, knon
    21292255              i = ni(j)
    2130               tke(i,k,nsrf)    = ytke(j,k)
    2131               tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j)
     2256!jyg<
     2257!!              tke(i,k,nsrf)    = ytke(j,k)
     2258!!              tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j)
     2259              tke_x(i,k,nsrf)    = ytke(j,k)
     2260              tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + ytke(j,k)*ypct(j)
     2261!>jyg
    21322262           END DO
    21332263        END DO
    … …  
    21382268            i = ni(j)
    21392269            wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)
    2140             tke(i,k,nsrf)   = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)
    2141             tke(i,k,is_ave) = tke(i,k,is_ave) + tke(i,k,nsrf)*ypct(j)
     2270!jyg<
     2271!!            tke(i,k,nsrf)   = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)
     2272!!            tke(i,k,is_ave) = tke(i,k,is_ave) + tke(i,k,nsrf)*ypct(j)
     2273            tke_x(i,k,nsrf)   = ytke_x(j,k)
     2274            tke_x(i,k,is_ave)   = tke_x(i,k,is_ave) + tke_x(i,k,nsrf)*ypct(j)
     2275            wake_dltke(i,k,is_ave)   = wake_dltke(i,k,is_ave) + wake_dltke(i,k,nsrf)*ypct(j)
     2276
     2277!>jyg
    21422278          ENDDO
    21432279        ENDDO
    … …  
    21632299       END IF
    21642300       
    2165        ftsoil(:,:,nsrf) = 0.
     2301!jyg<
     2302!!       ftsoil(:,:,nsrf) = 0.
     2303!>jyg
    21662304       DO k = 1, nsoilmx
    21672305          DO j = 1, knon
    … …  
    22212359!
    22222360!****************************************************************************************
    2223        t2m(:,nsrf)    = 0.
    2224        q2m(:,nsrf)    = 0.
    2225        ustar(:,nsrf)   = 0.
    2226        wstar(:,nsrf)   = 0.
    2227        u10m(:,nsrf)   = 0.
    2228        v10m(:,nsrf)   = 0.
    2229 
    2230        pblh(:,nsrf)   = 0.        ! Hauteur de couche limite
    2231        plcl(:,nsrf)   = 0.        ! Niveau de condensation de la CLA
    2232        capCL(:,nsrf)  = 0.        ! CAPE de couche limite
    2233        oliqCL(:,nsrf) = 0.        ! eau_liqu integree de couche limite
    2234        cteiCL(:,nsrf) = 0.        ! cloud top instab. crit. couche limite
    2235        pblt(:,nsrf)   = 0.        ! T a la Hauteur de couche limite
    2236        therm(:,nsrf)  = 0.
    2237        trmb1(:,nsrf)  = 0.        ! deep_cape
    2238        trmb2(:,nsrf)  = 0.        ! inhibition
    2239        trmb3(:,nsrf)  = 0.        ! Point Omega
    2240 !
    2241 !!! jyg le 07/02/2012
    2242        IF (iflag_split .eq.1) THEN
    2243        t2m_x(:,nsrf)    = 0.
    2244        q2m_x(:,nsrf)    = 0.
    2245        ustar_x(:,nsrf)   = 0.
    2246        wstar_x(:,nsrf)   = 0.
    2247        u10m_x(:,nsrf)   = 0.
    2248        v10m_x(:,nsrf)   = 0.
    2249                            
    2250        pblh_x(:,nsrf)   = 0.      ! Hauteur de couche limite
    2251        plcl_x(:,nsrf)   = 0.      ! Niveau de condensation de la CLA
    2252        capCL_x(:,nsrf)  = 0.      ! CAPE de couche limite
    2253        oliqCL_x(:,nsrf) = 0.      ! eau_liqu integree de couche limite
    2254        cteiCL_x(:,nsrf) = 0.      ! cloud top instab. crit. couche limite
    2255        pblt_x(:,nsrf)   = 0.      ! T a la Hauteur de couche limite
    2256        therm_x(:,nsrf)  = 0.     
    2257        trmb1_x(:,nsrf)  = 0.      ! deep_cape
    2258        trmb2_x(:,nsrf)  = 0.      ! inhibition
    2259        trmb3_x(:,nsrf)  = 0.      ! Point Omega
    2260 !
    2261        t2m_w(:,nsrf)    = 0.
    2262        q2m_w(:,nsrf)    = 0.
    2263        ustar_w(:,nsrf)   = 0.
    2264        wstar_w(:,nsrf)   = 0.
    2265        u10m_w(:,nsrf)   = 0.
    2266        v10m_w(:,nsrf)   = 0.
    2267                            
    2268        pblh_w(:,nsrf)   = 0.      ! Hauteur de couche limite
    2269        plcl_w(:,nsrf)   = 0.      ! Niveau de condensation de la CLA
    2270        capCL_w(:,nsrf)  = 0.      ! CAPE de couche limite
    2271        oliqCL_w(:,nsrf) = 0.      ! eau_liqu integree de couche limite
    2272        cteiCL_w(:,nsrf) = 0.      ! cloud top instab. crit. couche limite
    2273        pblt_w(:,nsrf)   = 0.      ! T a la Hauteur de couche limite
    2274        therm_w(:,nsrf)  = 0.     
    2275        trmb1_w(:,nsrf)  = 0.      ! deep_cape
    2276        trmb2_w(:,nsrf)  = 0.      ! inhibition
    2277        trmb3_w(:,nsrf)  = 0.      ! Point Omega
    2278 !!!     
    2279        ENDIF  ! (iflag_split .eq.1)
    22802361!!!
    22812362!
    … …  
    28762957    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: alb1, alb2
    28772958    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: ustar,u10m, v10m
    2878     REAL, DIMENSION(klon,klev+1,nbsrf), INTENT(INOUT) :: tke
     2959    REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(INOUT) :: tke
    28792960
    28802961! Local variables
    … …  
    29653046
    29663047END MODULE pbl_surface_mod
     3048

  • LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90

    r2159 r2181  
    278278      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zxfluxlat_x, zxfluxlat_w
    279279!$OMP THREADPRIVATE(zxfluxlat_x, zxfluxlat_w)
    280 ! Entrées supplémentaires couche-limite
    281       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: t_x, t_w
    282 !$OMP THREADPRIVATE(t_x, t_w)
    283       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: q_x, q_w
    284 !$OMP THREADPRIVATE(q_x, q_w)
    285 ! Sorties ferret
     280!jyg<
     281!!! Entrées supplémentaires couche-limite
     282!!      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: t_x, t_w
     283!!!$OMP THREADPRIVATE(t_x, t_w)
     284!!      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: q_x, q_w
     285!!!$OMP THREADPRIVATE(q_x, q_w)
     286!>jyg
     287!!! Sorties ferret
    286288      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: dtvdf_x, dtvdf_w
    287289!$OMP THREADPRIVATE(dtvdf_x, dtvdf_w)
    … …  
    516518      ALLOCATE(sens_x(klon), sens_w(klon))
    517519      ALLOCATE(zxfluxlat_x(klon), zxfluxlat_w(klon))
    518       ALLOCATE(t_x(klon,klev), t_w(klon,klev))
    519       ALLOCATE(q_x(klon,klev), q_w(klon,klev))
     520!jyg<
     521!!      ALLOCATE(t_x(klon,klev), t_w(klon,klev))
     522!!      ALLOCATE(q_x(klon,klev), q_w(klon,klev))
     523!>jyg
    520524      ALLOCATE(dtvdf_x(klon,klev), dtvdf_w(klon,klev))
    521525      ALLOCATE(dqvdf_x(klon,klev), dqvdf_w(klon,klev))
    … …  
    710714      DEALLOCATE(sens_x, sens_w)
    711715      DEALLOCATE(zxfluxlat_x, zxfluxlat_w)
    712       DEALLOCATE(t_x, t_w)
    713       DEALLOCATE(q_x, q_w)
     716!jyg<
     717!!      DEALLOCATE(t_x, t_w)
     718!!      DEALLOCATE(q_x, q_w)
     719!>jyg
    714720      DEALLOCATE(dtvdf_x, dtvdf_w)
    715721      DEALLOCATE(dqvdf_x, dqvdf_w)

  • LMDZ5/trunk/libf/phylmd/phys_state_var_mod.F90

    r2159 r2181  
    491491      ALLOCATE(dt_wake(klon,klev), dq_wake(klon,klev))
    492492!jyg<
    493       ALLOCATE(wake_delta_pbl_TKE(klon,klev+1,nbsrf))
     493      ALLOCATE(wake_delta_pbl_TKE(klon,klev+1,nbsrf+1))
    494494!>jyg
    495495      ALLOCATE(pfrac_impa(klon,klev), pfrac_nucl(klon,klev))

  • LMDZ5/trunk/libf/phylmd/physiq.F90

    r2180 r2181  
    18561856!
    18571857!  Add turbulent diffusion tendency to the wake difference variables
    1858     wake_deltat(:,:) = wake_deltat(:,:) + (d_t_vdf_w(:,:)-d_t_vdf_x(:,:))
    1859     wake_deltaq(:,:) = wake_deltaq(:,:) + (d_q_vdf_w(:,:)-d_q_vdf_x(:,:))
     1858    IF (mod(iflag_pbl_split,2) .NE. 0) THEN
     1859     wake_deltat(:,:) = wake_deltat(:,:) + (d_t_vdf_w(:,:)-d_t_vdf_x(:,:))
     1860     wake_deltaq(:,:) = wake_deltaq(:,:) + (d_q_vdf_w(:,:)-d_q_vdf_x(:,:))
     1861    ENDIF
    18601862
    18611863
Note: See TracChangeset for help on using the changeset viewer.