Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/Makefile =================================================================== --- /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/Makefile (revision 1754) +++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/Makefile (revision 1755) @@ -9,4 +9,5 @@ module_physics_addtendc.o \ module_physics_init.o \ + variables_mod.o \ update_inputs_physiq_mod.o \ update_outputs_physiq_mod.o \ Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/dynphy_wrf_mars_lmd_new/callphysiq_mod.F =================================================================== --- /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/dynphy_wrf_mars_lmd_new/callphysiq_mod.F (revision 1754) +++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/dynphy_wrf_mars_lmd_new/callphysiq_mod.F (revision 1755) @@ -13,13 +13,7 @@ SUBROUTINE call_physiq(planet_type, klon,llm,nqtot, & - debut_split,lafin_split, & - jD_cur,jH_cur_split,zdt_split, & - zplev_omp,zplay_omp, & - zpk_omp,zphi_omp,zphis_omp, & - presnivs_omp, & - zufi_omp,zvfi_omp,zrfi_omp,ztfi_omp,zqfi_omp, & - flxwfi_omp, & - zdufi_omp,zdvfi_omp,zdtfi_omp,zdqfi_omp,zdpsrf_omp) + debut_split,lafin_split) + USE variables_mod USE physiq_mod, ONLY: physiq IMPLICIT NONE @@ -32,25 +26,4 @@ LOGICAL,INTENT(IN) :: debut_split ! .true. if very first call to physics LOGICAL,INTENT(IN) :: lafin_split ! .true. if last call to physics - REAL,INTENT(IN) :: JD_cur ! Julian day - REAL,INTENT(IN) :: JH_cur_split ! Julian hour (fraction of day) - REAL,INTENT(IN) :: zdt_split ! time step over which the physics are evaluated - REAL,INTENT(IN) :: zplev_omp(klon,llm+1) ! interlayer pressure (Pa) - REAL,INTENT(IN) :: zplay_omp(klon,llm) ! mid-layer pressure (Pa) - REAL,INTENT(INOUT) :: zpk_omp(klon,llm) - REAL,INTENT(IN) :: zphi_omp(klon,llm) ! geopotential at midlayer - REAL,INTENT(INOUT) :: zphis_omp(klon) ! surface geopotential - REAL,INTENT(INOUT) :: presnivs_omp(llm) ! approximate pressure of atm. layers - REAL,INTENT(IN) :: zufi_omp(klon,llm) ! zonal wind (m/s) - REAL,INTENT(IN) :: zvfi_omp(klon,llm) ! meridional wind (m/s) - REAL,INTENT(INOUT) :: zrfi_omp(klon,llm) ! relative wind vorticity, in s-1 - REAL,INTENT(IN) :: ztfi_omp(klon,llm) ! temperature (K) - REAL,INTENT(IN) :: zqfi_omp(klon,llm,nqtot) ! tracers (*/kg of air) - REAL,INTENT(IN) :: flxwfi_omp(klon,llm) ! Vertical mass flux on lower mesh interfaces (kg/s) - ! tendencies (in */s) from the physics: - REAL,INTENT(OUT) :: zdufi_omp(klon,llm) ! tendency on zonal winds - REAL,INTENT(OUT) :: zdvfi_omp(klon,llm) ! tendency on meridional winds - REAL,INTENT(OUT) :: zdtfi_omp(klon,llm) ! tendency on temperature - REAL,INTENT(OUT) :: zdqfi_omp(klon,llm,nqtot) ! tendency on tracers - REAL,INTENT(OUT) :: zdpsrf_omp(klon) ! tendency on surface pressure ! ! Local variables Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/dynphy_wrf_mars_lmd_new/iniphysiq_mod.F =================================================================== --- /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/dynphy_wrf_mars_lmd_new/iniphysiq_mod.F (revision 1754) +++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/dynphy_wrf_mars_lmd_new/iniphysiq_mod.F (revision 1755) @@ -3,6 +3,6 @@ CONTAINS -subroutine iniphysiq(ngrid,nlayer,nq,phour_ini,piphysiq,& - punjours, pdayref,ptimestep, & +subroutine iniphysiq(ngrid,nlayer,nq,piphysiq,& + punjours, pdayref, & prad,pg,pr,pcpp,iflag_phys) @@ -11,4 +11,6 @@ dtphys,daysec,day_ini,hour_ini use update_inputs_physiq_mod, only: traceurs +USE variables_mod, only: phour_ini,zdt_split !! zdt_split <> pttimestep + !real*8,intent(in) :: ptimestep !physics time step (s) [dtphys] implicit none @@ -20,5 +22,5 @@ real,intent(in) :: punjours ! length (in s) of a standard day [daysec] integer,intent(in) :: pdayref ! reference day of for the simulation [day_ini] -real,intent(in) :: ptimestep !physics time step (s) [dtphys] +!real*8,intent(in) :: ptimestep !physics time step (s) [dtphys] integer,intent(in) :: iflag_phys ! type of physics to be called @@ -26,5 +28,5 @@ integer,intent(in) :: nlayer ! number of atmospheric layers integer,intent(in) :: nq ! number of tracers -real,intent(in) :: phour_ini ! start time (fraction of day) of the run 0= pday ! Julian day + !! JH_cur_split <> ptime ! Julian hour (fraction of day) + + implicit none INTEGER, INTENT(IN) :: JULDAY, JULYR REAL, INTENT(IN) :: GMT,elaps,lon_input,ls_input,lct_input - REAL,INTENT(OUT) :: pday,ptime,MY + REAL,INTENT(OUT) :: MY IF (JULYR .ne. 9999) THEN @@ -80,9 +86,9 @@ ! specified ! - ptime = (GMT + elaps/3700.) !! universal time (0 ptimestep + !real*8,intent(in) :: ptimestep !physics time step (s) [dtphys] + implicit none @@ -27,5 +30,4 @@ REAL,intent(in) :: punjours !DOUBLE PRECISION,intent(in) :: ptimestep -REAL*8,intent(in) :: phour_ini !real,intent(in) :: prad ! radius of the planet (m) @@ -35,5 +37,5 @@ !real,intent(in) :: punjours ! length (in s) of a standard day [daysec] integer,intent(in) :: pdayref ! reference day of for the simulation [day_ini] -real*8,intent(in) :: ptimestep !physics time step (s) [dtphys] +!real*8,intent(in) :: ptimestep !physics time step (s) [dtphys] integer,intent(in) :: iflag_phys ! type of physics to be called Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/dynphy_wrf_venus_lmd_new/update_inputs_physiq_mod.F =================================================================== --- /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/dynphy_wrf_venus_lmd_new/update_inputs_physiq_mod.F (revision 1754) +++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/dynphy_wrf_venus_lmd_new/update_inputs_physiq_mod.F (revision 1755) @@ -32,5 +32,9 @@ elaps,& lct_input,lon_input,ls_input,& - ptime,pday,MY) + MY) + + USE variables_mod, only: JD_cur,JH_cur_split,phour_ini + !! JD_cur <> pday ! Julian day + !! JH_cur_split <> ptime ! Julian hour (fraction of day) implicit none @@ -39,5 +43,4 @@ REAL, INTENT(IN) :: GMT,elaps,lon_input,ls_input,lct_input REAL,INTENT(OUT) :: MY - REAL*8,INTENT(OUT) :: ptime,pday ! @@ -45,17 +48,17 @@ ! IF (JULYR .ne. 9999) THEN - ptime = (GMT + elaps/420000.) !! universal time (0 inputs: INTEGER :: ngrid,nlayer,nq,nsoil - REAL*8 :: pday,ptime REAL :: MY REAL :: phisfi_val LOGICAL :: firstcall,lastcall ! ---------- - REAL*8,DIMENSION(:,:),ALLOCATABLE :: pplev,pplay,pphi,pu,pv,pt,flxw - REAL*8,DIMENSION(:,:,:),ALLOCATABLE :: pq ! <------ outputs: ! physical tendencies - REAL*8,DIMENSION(:),ALLOCATABLE :: pdpsrf - REAL*8,DIMENSION(:,:),ALLOCATABLE :: pdu,pdv,pdt,pdtheta - REAL*8,DIMENSION(:,:,:),ALLOCATABLE :: pdq + REAL*8,DIMENSION(:,:),ALLOCATABLE :: pdtheta ! ... intermediate arrays REAL, DIMENSION(:), ALLOCATABLE :: & @@ -160,5 +156,4 @@ INTEGER :: sponge_top,relax,ips,ipe,jps,jpe,kps,kpe REAL :: elaps - REAL*8 :: ptimestep INTEGER :: test REAL :: wappel_phys @@ -192,9 +187,4 @@ !!!IDEALIZED IDEALIZED - !! arguments to physiq - REAL*8,ALLOCATABLE :: zpk_omp(:,:) - REAL*8,ALLOCATABLE :: zphis_omp(:) ! surface geopotential - REAL*8,ALLOCATABLE :: presnivs_omp(:) ! approximate pressure of atm. layers - REAL*8,ALLOCATABLE :: zrfi_omp(:,:) ! relative wind vorticity, in s-1 REAL :: tk1,tk2 !================================================================== @@ -264,5 +254,5 @@ !day_ini = JULDAY - 1 !! GCM convention !! pday at firstcall is day_ini wappel_phys = RADT -ptimestep = dt*wappel_phys ! physical timestep (s) +zdt_split = dt*wappel_phys ! physical timestep (s) ngrid=(ipe-ips+1)*(jpe-jps+1) ! size of the horizontal grid nlayer = kpe-kps+1 ! number of vertical layers: nlayermx @@ -375,18 +365,19 @@ ! ALLOCATE ! !----------! -!-------------------------------------------------------------------------------! -! outputs: ! -! pdu(ngrid,nlayermx) \ ! -! pdv(ngrid,nlayermx) \ Temporal derivative of the corresponding ! -! pdt(ngrid,nlayermx) / variables due to physical processes. ! -! pdq(ngrid,nlayermx) / ! -! pdpsrf(ngrid) / ! -!-------------------------------------------------------------------------------! -ALLOCATE(pdpsrf(ngrid)) -ALLOCATE(pdu(ngrid,nlayer)) -ALLOCATE(pdv(ngrid,nlayer)) -ALLOCATE(pdt(ngrid,nlayer)) -ALLOCATE(pdtheta(ngrid,nlayer)) -ALLOCATE(pdq(ngrid,nlayer,nq)) +!!-------------------------------------------------------------------------------! +!! outputs: ! +!! pdu(ngrid,nlayermx) \ ! +!! pdv(ngrid,nlayermx) \ Temporal derivative of the corresponding ! +!! pdt(ngrid,nlayermx) / variables due to physical processes. ! +!! pdq(ngrid,nlayermx) / ! +!! pdpsrf(ngrid) / ! +!!-------------------------------------------------------------------------------! +!ALLOCATE(pdpsrf(ngrid)) +!ALLOCATE(pdu(ngrid,nlayer)) +!ALLOCATE(pdv(ngrid,nlayer)) +!ALLOCATE(pdt(ngrid,nlayer)) +!ALLOCATE(pdtheta(ngrid,nlayer)) +!ALLOCATE(pdq(ngrid,nlayer,nq)) +CALL allocate_interface(ngrid,nlayer,nq) !!! !!! BIG LOOP : 1. no call for physics, used saved values @@ -395,18 +386,18 @@ print *,'** ',planet_type,'** NO CALL FOR PHYSICS, go to next step...',test #ifdef SPECIAL_NEST_SAVE -pdpsrf(:)=dp_save(:,id) -pdu(:,:)=du_save(:,:,id) -pdv(:,:)=dv_save(:,:,id) -pdt(:,:)=dt_save(:,:,id) +zdpsrf_omp(:)=dp_save(:,id) +zdufi_omp(:,:)=du_save(:,:,id) +zdvfi_omp(:,:)=dv_save(:,:,id) +zdtfi_omp(:,:)=dt_save(:,:,id) pdtheta(:,:)=dtheta_save(:,:,id) -pdq(:,:,:)=dq_save(:,:,:,id) +zdqfi_omp(:,:,:)=dq_save(:,:,:,id) #else print*,'else' -pdpsrf(:)=dp_save(:) -pdu(:,:)=du_save(:,:) -pdv(:,:)=dv_save(:,:) -pdt(:,:)=dt_save(:,:) +zdpsrf_omp(:)=dp_save(:) +zdufi_omp(:,:)=du_save(:,:) +zdvfi_omp(:,:)=dv_save(:,:) +zdtfi_omp(:,:)=dt_save(:,:) pdtheta(:,:)=dtheta_save(:,:) -pdq(:,:,:)=dq_save(:,:,:) +zdqfi_omp(:,:,:)=dq_save(:,:,:) #endif !!! @@ -417,17 +408,4 @@ ! ALLOCATE ! !----------! -! inputs ... -ALLOCATE(pplev(ngrid,nlayer+1)) !!!!! -ALLOCATE(pplay(ngrid,nlayer)) !!!!! -ALLOCATE(pphi(ngrid,nlayer)) !!!!! -ALLOCATE(pu(ngrid,nlayer)) !!!!! -ALLOCATE(pv(ngrid,nlayer)) !!!!! -ALLOCATE(pt(ngrid,nlayer)) !!!!! -ALLOCATE(flxw(ngrid,nlayer)) !!!!! -ALLOCATE(pq(ngrid,nlayer,nq)) !!!!! -ALLOCATE(zpk_omp(ngrid,nlayer)) -ALLOCATE(zphis_omp(ngrid)) -ALLOCATE(presnivs_omp(nlayer)) -ALLOCATE(zrfi_omp(ngrid,nlayer)) ! interm ALLOCATE(dz8w_prof(nlayer)) @@ -461,12 +439,14 @@ elaps,& lct_input,lon_input,ls_input,& - ptime,pday,MY) + MY) !! Fill planetary parameters in modules !! Values defined in the module_model_constants.F WRF routine CALL update_inputs_physiq_constants !! Initialize physics - CALL iniphysiq(ngrid,nlayer,nq,ptime,wappel_phys,& - wdaysec,floor(pday),ptimestep, & + CALL iniphysiq(ngrid,nlayer,nq,wappel_phys,& + wdaysec,floor(JD_cur), & 1./reradius,g,r_d,cp,1) + !! Set up initial time + phour_ini = JH_cur_split ENDIF allocation_firstcall @@ -492,11 +472,11 @@ !--------------------------------------! dz8w_prof(:) = dz8w(i,kps:kpe,j) ! dz between full levels (m) -p8w_prof(:) = p8w(i,kps:kpe,j) ! pressure full level (Pa) >> pplev -p_prof(:) = p(i,kps:kpe,j) ! pressure half level (Pa) >> pplay +p8w_prof(:) = p8w(i,kps:kpe,j) ! pressure full level (Pa) >> zplev_omp +p_prof(:) = p(i,kps:kpe,j) ! pressure half level (Pa) >> zplay_omp t_prof(:) = t(i,kps:kpe,j) ! temperature half level (K) >> pt t8w_prof(:) = t8w(i,kps:kpe,j) ! temperature full level (K) -u_prof(:) = u(i,kps:kpe,j) ! zonal wind (A-grid: unstaggered) half level (m/s) >> pu +u_prof(:) = u(i,kps:kpe,j) ! zonal wind (A-grid: unstaggered) half level (m/s) >> zufi_omp v_prof(:) = v(i,kps:kpe,j) ! meridional wind (A-grid: unstaggered) half level (m/s) >> pv -z_prof(:) = z(i,kps:kpe,j) ! geopotential height half level (m) >> pphi/g +z_prof(:) = z(i,kps:kpe,j) ! geopotential height half level (m) >> zphi_omp/g !--------------------------------! @@ -544,22 +524,22 @@ ! expressed with respect to the local surface ! !---------------------------------------------! -pphi(subs,:) = g*( z_prof(:)-(z_prof(1)-dz8w_prof(1)/2.) ) +zphi_omp(subs,:) = g*( z_prof(:)-(z_prof(1)-dz8w_prof(1)/2.) ) !--------------------------------! ! Dynamic fields for LMD physics ! !--------------------------------! -pplev(subs,1:nlayer) = p8w_prof(1:nlayer) !! NB: last level: no data -pplay(subs,:) = p_prof(:) -pt(subs,:) = t_prof(:) -pu(subs,:) = u_prof(:) -pv(subs,:) = v_prof(:) -flxw(subs,:) = 0 !! NB: not used in the physics, only diagnostic... +zplev_omp(subs,1:nlayer) = p8w_prof(1:nlayer) !! NB: last level: no data +zplay_omp(subs,:) = p_prof(:) +ztfi_omp(subs,:) = t_prof(:) +zufi_omp(subs,:) = u_prof(:) +zvfi_omp(subs,:) = v_prof(:) +flxwfi_omp(subs,:) = 0 !! NB: not used in the physics, only diagnostic... !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! for IDEALIZED CASES ONLY -IF (JULYR .eq. 9999) pplev(subs,nlayer+1)=0. !! pplev(subs,nlayer+1)=ptop >> NO ! +IF (JULYR .eq. 9999) zplev_omp(subs,nlayer+1)=0. !! zplev_omp(subs,nlayer+1)=ptop >> NO ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! NOTE: -! IF ( pplev(subs,nlayer+1) .le. 0 ) pplev(subs,nlayer+1)=ptop +! IF ( zplev_omp(subs,nlayer+1) .le. 0 ) zplev_omp(subs,nlayer+1)=ptop ! cree des diagnostics delirants et aleatoires dans le transfert radiatif @@ -567,5 +547,5 @@ ! Tracers ! !---------! -pq(subs,:,:) = q_prof(:,:) !! traceurs generiques, seuls noms sont specifiques +zqfi_omp(subs,:,:) = q_prof(:,:) !! traceurs generiques, seuls noms sont specifiques ENDDO @@ -659,10 +639,10 @@ call_physics : IF (wappel_phys .ne. 0.) THEN !!! initialize tendencies (security) -pdpsrf(:)=0. -pdu(:,:)=0. -pdv(:,:)=0. -pdt(:,:)=0. +zdpsrf_omp(:)=0. +zdufi_omp(:,:)=0. +zdvfi_omp(:,:)=0. +zdtfi_omp(:,:)=0. pdtheta(:,:)=0. -pdq(:,:,:)=0. +zdqfi_omp(:,:,:)=0. print *, '** ',planet_type,'** CALL TO LMD PHYSICS' !!! @@ -671,15 +651,8 @@ elaps,& lct_input,lon_input,ls_input,& - ptime,pday,MY) + MY) !!! CALL call_physiq(planet_type,ngrid,nlayer,nq, & - firstcall,lastcall, & - pday,ptime,ptimestep, & - pplev,pplay, & - zpk_omp,pphi,zphis_omp, & - presnivs_omp, & - pu,pv,zrfi_omp,pt,pq, & - flxw, & - pdu,pdv,pdt,pdq,pdpsrf) + firstcall,lastcall) !!! @@ -687,5 +660,5 @@ #ifdef DUSTSTORM IF (firstcall .EQV. .true.) THEN - pdq(:,:,:) = pdq(:,:,:) / dt + zdqfi_omp(:,:,:) = zdqfi_omp(:,:,:) / dt ENDIF #endif @@ -696,6 +669,6 @@ do k=kps,kpe subs=(j-jps)*(ipe-ips+1)+(i-ips+1) - tk1=(pt(subs,k)**nu + nu*TT00**nu*log((p1000mb/pplay(subs,k))**rcp))**(1/nu) - tk2=((pt(subs,k) + pdt(subs,k))**nu + nu*TT00**nu*log((p1000mb/pplay(subs,k))**rcp))**(1/nu) + tk1=(ztfi_omp(subs,k)**nu + nu*TT00**nu*log((p1000mb/zplay_omp(subs,k))**rcp))**(1/nu) + tk2=((ztfi_omp(subs,k) + zdtfi_omp(subs,k))**nu + nu*TT00**nu*log((p1000mb/zplay_omp(subs,k))**rcp))**(1/nu) pdtheta(subs,k)=tk2-tk1 enddo @@ -705,16 +678,4 @@ print *, '** ',planet_type,'** CALL TO LMD PHYSICS DONE' -DEALLOCATE(pplev) -DEALLOCATE(pplay) -DEALLOCATE(pphi) -DEALLOCATE(pu) -DEALLOCATE(pv) -DEALLOCATE(pt) -DEALLOCATE(flxw) -DEALLOCATE(pq) -DEALLOCATE(zpk_omp) -DEALLOCATE(zphis_omp) -DEALLOCATE(presnivs_omp) -DEALLOCATE(zrfi_omp) !---------------------------------------------------------------------------------! @@ -722,17 +683,17 @@ !---------------------------------------------------------------------------------! #ifdef SPECIAL_NEST_SAVE -dp_save(:,id)=pdpsrf(:) -du_save(:,:,id)=pdu(:,:) -dv_save(:,:,id)=pdv(:,:) -dt_save(:,:,id)=pdt(:,:) +dp_save(:,id)=zdpsrf_omp(:) +du_save(:,:,id)=zdufi_omp(:,:) +dv_save(:,:,id)=zdvfi_omp(:,:) +dt_save(:,:,id)=zdtfi_omp(:,:) dtheta_save(:,:,id)=pdtheta(:,:) -dq_save(:,:,:,id)=pdq(:,:,:) +dq_save(:,:,:,id)=zdqfi_omp(:,:,:) #else -dp_save(:)=pdpsrf(:) -du_save(:,:)=pdu(:,:) -dv_save(:,:)=pdv(:,:) -dt_save(:,:)=pdt(:,:) +dp_save(:)=zdpsrf_omp(:) +du_save(:,:)=zdufi_omp(:,:) +dv_save(:,:)=zdvfi_omp(:,:) +dt_save(:,:)=zdtfi_omp(:,:) dtheta_save(:,:)=pdtheta(:,:) -dq_save(:,:,:)=pdq(:,:,:) +dq_save(:,:,:)=zdqfi_omp(:,:,:) #endif @@ -806,7 +767,7 @@ ! zonal wind - RUBLTEN(i,kps:kpe,j) = pdu(subs,kps:kpe) + RUBLTEN(i,kps:kpe,j) = zdufi_omp(subs,kps:kpe) ! meridional wind - RVBLTEN(i,kps:kpe,j) = pdv(subs,kps:kpe) + RVBLTEN(i,kps:kpe,j) = zdvfi_omp(subs,kps:kpe) ! potential temperature ! (dT = dtheta * exner for isobaric coordinates or if pressure variations are negligible) @@ -814,9 +775,9 @@ RTHBLTEN(i,kps:kpe,j) = pdtheta(subs,kps:kpe) ELSE - RTHBLTEN(i,kps:kpe,j) = pdt(subs,kps:kpe) / exner(i,kps:kpe,j) + RTHBLTEN(i,kps:kpe,j) = zdtfi_omp(subs,kps:kpe) / exner(i,kps:kpe,j) ENDIF ! update surface pressure (cf CO2 cycle in physics) ! here dt is needed - PSFC(i,j)=PSFC(i,j)+pdpsrf(subs)*dt + PSFC(i,j)=PSFC(i,j)+zdpsrf_omp(subs)*dt ! tracers SCALAR(i,kps:kpe,j,1)=0. @@ -831,14 +792,10 @@ SCALAR(i,1,j,2) = SCALAR(i,1,j,2) + 1. !! this tracer is emitted in the surface layer CASE DEFAULT - SCALAR(i,kps:kpe,j,2:nq+1)=SCALAR(i,kps:kpe,j,2:nq+1)+pdq(subs,kps:kpe,1:nq)*dt !!! here dt is needed + SCALAR(i,kps:kpe,j,2:nq+1)=SCALAR(i,kps:kpe,j,2:nq+1)+zdqfi_omp(subs,kps:kpe,1:nq)*dt !!! here dt is needed END SELECT ENDDO ENDDO -DEALLOCATE(pdpsrf) -DEALLOCATE(pdu) -DEALLOCATE(pdv) -DEALLOCATE(pdt) -DEALLOCATE(pdq) +CALL deallocate_interface DEALLOCATE(pdtheta) !!*****!!