Changeset 234 for trunk/LMDZ.MARS/libf


Ignore:
Timestamp:
Jul 19, 2011, 5:25:58 PM (13 years ago)
Author:
aslmd
Message:

LMDZ.MARS
+ MESOSCALE

  • An important change to merge physiq.F and inifis.F for GCM and mesoscale
  • This is mostly transparent to GCM users and developers (use of MESOSCALE precompiler flags)
  • Makes it easy (and mostly automatic!) for changes in GCM physics to be impacted in mesoscale physics
  • A few minor changes have followed in the GCM (slope scheme, one-point diagnostic).
  • Compilation + run is OK on both sides (GCM and mesoscale).
  • On the mesoscale side, call_meso_physiq?.inc and call_meso_inifis?.inc have been changed accordingly.

Here is the excerpt from README file:

19/07/2011 == AS

  • Finished converging meso_physiq.F and meso_inifis.F towards physiq.F and inifis.F --> see previous point 15/07/2011 --> meso_ routines no longer exist (everything is in meso_inc and transparent to GCM users) --> GCM routines include mesoscale parts within MESOSCALE precompiler commands --> MESOSCALE parts are as hidden as possible not to mess up with GCM users/developers
  • Cleaned inelegant or useless #ifdef [or] #ifndef MESOSCALE in physiq and inifis so that a minimum amount of such precompiler commands is now reached [mainly related to I/O]
  • Added the SF08 slope insolation model in the general physics parameterizations. Added a callslope keyword in inifis.F and callkeys.h --> This keyword is False by default / True if you use -DMESOSCALE
  • Removed the obsolete call to Viking Lander 1 diagnostic Replaced it with a diagnostic for opacity at the domain center [valid for GCM and mesoscale]
Location:
trunk/LMDZ.MARS/libf/phymars
Files:
1 deleted
5 edited
2 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/callkeys.h

    r161 r234  
    1212     &   ,callg2d,linear,rayleigh,tracer,active,doubleq,submicron       &
    1313     &   ,lifting,callddevil,scavenging,sedimentation,activice,water    &
    14      &   ,caps,photochem,calltherm,outptherm
     14     &   ,caps,photochem,calltherm,outptherm,callslope
    1515     
    1616      COMMON/callkeys_i/iradia,iaervar,iddist,ilwd,ilwb,ilwn,ncouche    &
    … …  
    2424     &   ,callnirco2,callnlte,callthermos,callconduct,                  &
    2525     &    calleuv,callmolvis,callmoldiff,thermochem,thermoswater        &
    26      &   ,calltherm,outptherm
     26     &   ,calltherm,outptherm,callslope
    2727
    2828

  • trunk/LMDZ.MARS/libf/phymars/inifis.F

    r233 r234  
    1       SUBROUTINE meso_inifis(
     1      SUBROUTINE inifis(
    22     $           ngrid,nlayer
    3 #ifdef MESOSCALE
    4      $           ,nq,wdt,wday_ini,wdaysec,wappel_phys
    5 #else
    63     $           ,day_ini,pdaysec,ptimestep
    7 #endif
    84     $           ,plat,plon,parea
    95     $           ,prad,pg,pr,pcpp
    … …  
    2824!              stored in the q(:,:,:,:) array
    2925!             E.M. (june 2009) use getin routine to load parameters
    30 !             adapted to the WRF use - Aymeric Spiga - Jan 2009 - Jan 2007
     26!             adapted to the mesoscale use - Aymeric Spiga - 01/2007-07/2011
    3127!
    3228!
    … …  
    6460#include "yomaer.h"
    6561#include "datafile.h"
     62#include "slope.h"
    6663#ifdef MESOSCALE
    67 #include "slope.h"
    6864#include "comsoil.h"     !!MESOSCALE -- needed to fill volcapa
    6965#include "meso_inc/meso_inc_inifisvar.F"
    7066#endif
    7167      REAL prad,pg,pr,pcpp,pdaysec
    72 #ifndef MESOSCALE
     68
    7369      REAL ptimestep
    7470      INTEGER day_ini
    75 #endif
     71
    7672      INTEGER ngrid,nlayer
    7773      REAL plat(ngrid),plon(ngrid),parea(ngridmx)
    … …  
    9591      r=pr
    9692      rcp=r/cpp
    97 #ifndef MESOSCALE
    9893      daysec=pdaysec
    9994      dtphys=ptimestep
    100 #else
     95#ifdef MESOSCALE
    10196#include "meso_inc/meso_inc_inifisini.F"
    10297#endif
    … …  
    210205         call getin("callrad",callrad)
    211206         write(*,*) " callrad = ",callrad
     207
     208         write(*,*) "call slope insolation scheme ?",
     209     &              "(matters only if callrad=T)"
     210#ifdef MESOSCALE
     211         callslope=.true. ! default value
     212#else
     213         callslope=.false. ! default value (not supported yet)
     214#endif
     215         call getin("callslope",callslope)
     216         write(*,*) " callslope = ",callslope
    212217
    213218         write(*,*) "call NLTE radiative schemes ?",

  • trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_inifisini.F

    r227 r234  
    88c     
    99      ! in 'comcstfi.h'
    10       daysec=wdaysec 
    1110      omeg=womeg               
    1211      mugaz=wmugaz 
    … …  
    8988c It must be set now, because it is used afterwards
    9089c*****************************************************
    91         dtphys=wdt*float(wappel_phys)
     90        dtphys=wdt*float(ptimestep)
    9291        print*,'Physical timestep (s) ',dtphys

  • trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_inifisinvar.F

    r226 r234  
     1     $           ,nq,wdt
    12     $           ,womeg,wmugaz
    23     $           ,wyear_day,wperiheli,waphelie,wperi_day,wobliquit

  • trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_inifisvar.F

    r226 r234  
    1       INTEGER nq,wday_ini
     1      INTEGER nq
     2      REAL wdt
    23
    3       REAL womeg,wmugaz,wdaysec
     4      REAL womeg,wmugaz
    45      REAL wyear_day,wperiheli,waphelie,wperi_day,wobliquit
    56      REAL wz0,wemin_turb,wlmixmin
    … …  
    1213      REAL wtheta(ngrid),wpsi(ngrid)
    1314      REAL wvolcapa
    14       REAL wdt
    15       INTEGER wappel_phys

  • trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_var.F

    r226 r234  
    1010      REAL output_tab2d(ngridmx,n2d)
    1111      REAL output_tab3d(ngridmx,nlayer,n3d)
    12       REAL sl_ls, sl_lct, sl_lat, sl_tau, sl_alb, sl_the, sl_psi
    13       REAL sl_fl0, sl_flu
    14       REAL sl_ra, sl_di0
    15       REAL sky
    1612      REAL hfx(ngridmx)    !! pour LES avec isfflx!=0
    1713      REAL ust(ngridmx)    !! pour LES avec isfflx!=0

  • trunk/LMDZ.MARS/libf/phymars/physiq.F

    r233 r234  
    1       SUBROUTINE meso_physiq(
     1      SUBROUTINE physiq(
    22     $            ngrid,nlayer,nq
    33     $            ,firstcall,lastcall
    … …  
    6464c           Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
    6565c             Nb: See callradite.F for more information.
    66 c           Mesoscale version: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags
     66c           Mesoscale lines: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags
    6767c           
    6868c   arguments:
    … …  
    135135#include "netcdf.inc"
    136136
     137#include "slope.h"
     138
    137139#ifdef MESOSCALE
    138 #include "slope.h"
    139140#include "wrf_output_2d.h"
    140141#include "wrf_output_3d.h"
    … …  
    191192      REAL qsurf(ngridmx,nqmx)       ! tracer on surface (e.g. kg.m-2)
    192193      REAL q2(ngridmx,nlayermx+1)    ! Turbulent Kinetic Energy
    193       INTEGER ig_vl1                 ! Grid Point near VL1   (for diagnostic)
    194194
    195195c     Variables used by the water ice microphysical scheme:
    … …  
    201201      REAL, PARAMETER :: alb_surfice = 0.45 !!TESTS_JB
    202202
     203c     Variables used by the slope model
     204      REAL sl_ls, sl_lct, sl_lat
     205      REAL sl_tau, sl_alb, sl_the, sl_psi
     206      REAL sl_fl0, sl_flu
     207      REAL sl_ra, sl_di0
     208      REAL sky
     209
    203210      SAVE day_ini, icount
    204211      SAVE aerosol, tsurf,tsoil
    205212      SAVE co2ice,albedo,emis, q2
    206213      SAVE capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
    207       SAVE ig_vl1
    208214
    209215      REAL stephan   
    … …  
    316322      REAL lmax_th_out(ngridmx),zmax_th(ngridmx)
    317323      REAL wmax_th(ngridmx)
    318       REAL ,SAVE :: hfmax_th(ngridmx)
     324      REAL, SAVE :: hfmax_th(ngridmx)
    319325      REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx)
    320326      REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nqmx)
    … …  
    390396#ifdef MESOSCALE
    391397#include "meso_inc/meso_inc_caps.F"
    392 #endif
    393 
    394 #ifndef MESOSCALE
    395 c        Determining gridpoint near Viking Lander 1 (used for diagnostic only)
    396 c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    397 
    398          if(ngrid.ne.1) then
    399            latvl1= 22.27
    400            lonvl1= -47.94
    401            ig_vl1= 1+ int( (1.5-(latvl1-90.)*jjm/180.)  -2 )*iim +
    402      &              int(1.5+(lonvl1+180)*iim/360.)
    403            write(*,*) 'Viking Lander 1 GCM point: lat,lon',
    404      &              lati(ig_vl1)*180/pi, long(ig_vl1)*180/pi
    405          end if
    406398#endif
    407399
    … …  
    515507         IF( MOD(icount-1,iradia).EQ.0) THEN
    516508
    517 #ifdef MESOSCALE
    518            write (*,*) 'call radiative transfer'
    519 #endif
    520509c          Local Solar zenith angle
    521510c          ~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    558547     &     tauref,tau,aerosol,ccn,rdust,rice,nuice)
    559548
    560 #ifdef MESOSCALE
    561 #include "meso_inc/meso_inc_slope.F"
    562 #endif
     549c          Outputs for basic check (middle of domain)
     550c          ------------------------------------------
     551           print*, 'check lat lon', lati(igout)*180/pi,
     552     .                              long(igout)*180/pi
     553           print*, 'Ls =',zls*180./pi
     554           print*, 'tauref(700 Pa) =',tauref(igout)
     555           print*, 'tau(700 Pa) =',tau(igout,1)*700./pplev(igout,1)
     556
     557c          ---------------------------------------------------------
     558c          Call slope parameterization for direct and scattered flux
     559c          ---------------------------------------------------------
     560           IF(callslope) THEN
     561            print *, 'Slope scheme is on and computing...'
     562            DO ig=1,ngrid 
     563              sl_the = theta_sl(ig)
     564              IF (sl_the .ne. 0.) THEN
     565                ztim1=fluxsurf_sw(ig,1)+fluxsurf_sw(ig,2)
     566                DO l=1,2
     567                 sl_lct = ptime*24. + 180.*long(ig)/pi/15.
     568                 sl_ra  = pi*(1.0-sl_lct/12.)
     569                 sl_lat = 180.*lati(ig)/pi
     570                 sl_tau = tau(ig,1)
     571                 sl_alb = albedo(ig,l)
     572                 sl_psi = psi_sl(ig)
     573                 sl_fl0 = fluxsurf_sw(ig,l)
     574                 sl_di0 = 0.
     575                 if (mu0(ig) .gt. 0.) then
     576                  sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
     577                  sl_di0 = sl_di0*1370./dist_sol/dist_sol
     578                  sl_di0 = sl_di0/ztim1
     579                  sl_di0 = fluxsurf_sw(ig,l)*sl_di0
     580                 endif
     581                 ! you never know (roundup concern...)
     582                 if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
     583                 !!!!!!!!!!!!!!!!!!!!!!!!!!
     584                 CALL param_slope( mu0(ig), declin, sl_ra, sl_lat,
     585     &                             sl_tau, sl_alb, sl_the, sl_psi,
     586     &                             sl_di0, sl_fl0, sl_flu )
     587                 !!!!!!!!!!!!!!!!!!!!!!!!!!
     588                 fluxsurf_sw(ig,l) = sl_flu
     589                ENDDO
     590              !!! compute correction on IR flux as well
     591              sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
     592              fluxsurf_lw(ig)= fluxsurf_lw(ig)*sky
     593              ENDIF
     594            ENDDO
     595           ENDIF
    563596
    564597c          CO2 near infrared absorption
    … …  
    606639     $         stephan*zplanck(ig)*zplanck(ig)
    607640               fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
    608 #ifdef MESOSCALE
    609                !!!! param slope
    610                sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
    611                fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig)
    612 #endif
     641               IF(callslope) THEN
     642                 sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
     643                 fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig)
     644               ENDIF
    613645           ENDDO
    614646
    … …  
    621653      ENDIF ! of IF (callrad)
    622654
    623 #ifndef MESOSCALE
    624655c-----------------------------------------------------------------------
    625656c    3. Gravity wave and subgrid scale topography drag :
    … …  
    640671        ENDDO
    641672      ENDIF
    642 #endif
     673
    643674c-----------------------------------------------------------------------
    644675c    4. Vertical diffusion (turbulent mixing):
    … …  
    11451176        ENDDO
    11461177      ENDDO
    1147       if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then       
     1178      if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then
    11481179        write(*,*) 'PHYSIQ: stability WARNING :'
    11491180        write(*,*) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
    … …  
    11721203
    11731204      IF (ngrid.NE.1) THEN
    1174          print*,'Ls =',zls*180./pi
    1175      &   ,' tauref(700 Pa,lat=0) =',tauref(ngrid/2)
    1176 #ifndef MESOSCALE
    1177      &   ,' tau(Viking1) =',tau(ig_vl1,1)
    1178 #endif
    11791205
    11801206#ifndef MESOSCALE
    … …  
    12411267         endif ! of if (tracer)
    12421268
    1243 #ifndef MESOSCALE
    12441269c        -----------------------------------------------------------------
    12451270c        WSTATS: Saving statistics
    … …  
    13541379            CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
    13551380         ENDIF
    1356 #endif
     1381
    13571382
    13581383#ifdef MESOSCALE
     1384      !!!
     1385      !!! OUTPUT FIELDS
     1386      !!!
    13591387      wtsurf(1:ngrid) = tsurf(1:ngrid)    !! surface temperature
    13601388      wco2ice(1:ngrid) = co2ice(1:ngrid)  !! co2 ice
    … …  
    15791607      ELSE     ! if(ngrid.eq.1)
    15801608
    1581 #ifndef MESOSCALE
    15821609         print*,'Ls =',zls*180./pi,
    15831610     &  '  tauref(700 Pa) =',tauref
    … …  
    16691696     &                   log(pplay(1,nlayer)/pplay(1,nlayer-1))*
    16701697     &                   rnew(1,nlayer)*tmean/g
    1671 #endif
    16721698
    16731699      END IF       ! if(ngrid.ne.1)
    16741700
    16751701      icount=icount+1
    1676 #ifdef MESOSCALE
    1677       write(*,*) 'now, back to the dynamical core...'
    1678 #endif
    16791702      RETURN
    16801703      END
Note: See TracChangeset for help on using the changeset viewer.