Changeset 899 for trunk/LMDZ.MARS/libf


Ignore:
Timestamp:
Mar 12, 2013, 1:15:20 PM (12 years ago)
Author:
emillour
Message:

Mars GCM:

  • added possibility to set slope parameters in testphys1d

EM

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r720 r899  
    2929#include "comgeomfi.h"
    3030#include "surfdat.h"
     31#include "slope.h"
    3132#include "comsoil.h"
    3233#include "comdiurn.h"
    … …  
    100101
    101102c ------------------------------------------------------
    102 c  Constantes prescrites ICI
     103c  Prescribed constants to be set here
    103104c ------------------------------------------------------
    104105
    105106      pi=2.E+0*asin(1.E+0)
    106107
    107 c     Constante de la Planete Mars
     108c     Mars planetary constants
    108109c     ----------------------------
    109       rad=3397200.               ! rayon de mars (m)  ~3397200 m
    110       daysec=88775.              ! duree du sol (s)  ~88775 s
    111       omeg=4.*asin(1.)/(daysec)  ! vitesse de rotation (rad.s-1)
    112       g=3.72                     ! gravite (m.s-2) ~3.72 
    113       mugaz=43.49                ! Masse molaire de l'atm (g.mol-1) ~43.49
     110      rad=3397200.               ! mars radius (m)  ~3397200 m
     111      daysec=88775.              ! length of a sol (s)  ~88775 s
     112      omeg=4.*asin(1.)/(daysec)  ! rotation rate (rad.s-1)
     113      g=3.72                     ! gravity (m.s-2) ~3.72 
     114      mugaz=43.49                ! atmosphere mola mass (g.mol-1) ~43.49
    114115      rcp=.256793                ! = r/cp  ~0.256793
    115116      r= 8.314511E+0 *1000.E+0/mugaz
    116117      cpp= r/rcp
    117       year_day = 669           ! duree de l'annee (sols) ~668.6
    118       periheli = 206.66          ! dist.min. soleil-mars (Mkm) ~206.66
    119       aphelie = 249.22           ! dist.max. soleil-mars (Mkm) ~249.22
    120       peri_day =  485.           ! date du perihelie (sols depuis printemps)
    121       obliquit = 25.2           ! Obliquite de la planete (deg) ~25.2         
    122  
    123 c     Parametres Couche limite et Turbulence
    124 c     --------------------------------------
    125       z0_default =  1.e-2                ! surface roughness (m) ~0.01
    126       emin_turb = 1.e-6          ! energie minimale ~1.e-8
    127       lmixmin = 30               ! longueur de melange ~100
    128  
    129 c     propriete optiques des calottes et emissivite du sol
     118      year_day = 669             ! lenght of year (sols) ~668.6
     119      periheli = 206.66          ! minimum sun-mars distance (Mkm) ~206.66
     120      aphelie = 249.22           ! maximum sun-mars distance (Mkm) ~249.22
     121      peri_day =  485.           ! perihelion date (sols since N. Spring)
     122      obliquit = 25.2            ! Obliquity (deg) ~25.2         
     123 
     124c     Planetary Boundary Layer and Turbulence parameters
     125c     --------------------------------------------------
     126      z0_default =  1.e-2        ! surface roughness (m) ~0.01
     127      emin_turb = 1.e-6          ! minimal turbulent energy ~1.e-8
     128      lmixmin = 30               ! mixing length ~100
     129 
     130c     cap properties and surface emissivities
    130131c     ----------------------------------------------------
    131       emissiv= 0.95              ! Emissivite du sol martien ~.95
    132       emisice(1)=0.95            ! Emissivite calotte nord
    133       emisice(2)=0.95            ! Emissivite calotte sud 
    134       albedice(1)=0.5           ! Albedo calotte nord
    135       albedice(2)=0.5            ! Albedo calotte sud
     132      emissiv= 0.95              ! Bare ground emissivity ~.95
     133      emisice(1)=0.95            ! Northern cap emissivity
     134      emisice(2)=0.95            ! Southern cap emisssivity
     135      albedice(1)=0.5            ! Northern cap albedo
     136      albedice(2)=0.5            ! Southern cap albedo
    136137      iceradius(1) = 100.e-6     ! mean scat radius of CO2 snow (north)
    137138      iceradius(2) = 100.e-6     ! mean scat radius of CO2 snow (south)
    138139      dtemisice(1) = 2.          ! time scale for snow metamorphism (north)
    139140      dtemisice(2) = 2.          ! time scale for snow metamorphism (south
    140       hybrid=.false.
     141
    141142
    142143c ------------------------------------------------------
    143 c  Lecture des parametres dans "run.def"
     144c  Loading run parameters from "run.def" file
    144145c ------------------------------------------------------
    145146
    … …  
    354355          tnom(iq)=str7
    355356        enddo
     357      ! and just to be clean, also initialize tracers to zero for physdem1
     358        q(:,:)=0
     359        qsurf(:)=0     
    356360      endif ! of if (tracer)
    357361     
    … …  
    359363      !write(*,*) "testphys1d qsurf", qsurf
    360364
    361 c  Date et heure locale du debut du run
    362 c  ------------------------------------
    363 c    Date (en sols depuis le solstice de printemps) du debut du run
     365c  Date and local time at beginning of run
     366c  ---------------------------------------
     367c    Date (in sols since spring solstice) at beginning of run
    364368      day0 = 0 ! default value for day0
    365369      write(*,*) 'Initial date (in martian sols ; =0 at Ls=0)?'
    … …  
    367371      day=float(day0)
    368372      write(*,*) " day0 = ",day0
    369 c  Heure de demarrage
     373c  Local time at beginning of run
    370374      time=0 ! default value for time
    371375      write(*,*)'Initial local time (in hours, between 0 and 24)?'
    … …  
    374378      time=time/24.E+0 ! convert time (hours) to fraction of sol
    375379
    376 c  Discretisation (Definition de la grille et des pas de temps)
     380c  Discretization (Definition of grid and time steps)
    377381c  --------------
    378382c
    … …  
    393397      ndt=ndt*day_step     
    394398      dtphys=daysec/day_step 
    395 c Pression de surface sur la planete
     399
     400c Imposed surface pressure
    396401c ------------------------------------
    397402c
    … …  
    400405      call getin("psurf",psurf)
    401406      write(*,*) " psurf = ",psurf
    402 c Pression de reference
    403       pa=20.
    404       preff=610.     
    405  
    406 c Proprietes des poussiere aerosol
     407c Reference pressures
     408      pa=20.   ! transition pressure (for hybrid coord.)
     409      preff=610.      ! reference surface pressure
     410 
     411c Aerosol properties
    407412c --------------------------------
    408       tauvis=0.2 ! default value for tauvis
     413      tauvis=0.2 ! default value for tauvis (dust opacity)
    409414      write(*,'("Reference dust opacity at ",f4.0," Pa ?")')odpref
    410415      call getin("tauvis",tauvis)
    … …  
    439444      long(1)=long(1)*pi/180.E+0
    440445
    441 c  Initialisation albedo / inertie du sol
    442 c  --------------------------------------
     446c  Initialize albedo / soil thermal inertia
     447c  ----------------------------------------
    443448c
    444449      albedodat(1)=0.2 ! default value for albedodat
    … …  
    456461      call getin("z0",z0(1))
    457462      write(*,*) " z0 = ",z0(1)
    458 c
    459 c  pour le schema d'ondes de gravite
     463
     464! Initialize local slope parameters (only matters if "callslope"
     465! is .true. in callphys.def)
     466      ! slope inclination angle (deg) 0: horizontal, 90: vertical
     467      theta_sl(1)=0.0 ! default: no inclination
     468      call getin("slope_inclination",theta_sl(1))
     469      ! slope orientation (deg)
     470      ! 0 == Northward, 90 == Eastward, 180 == Southward, 270 == Westward
     471      psi_sl(1)=0.0 ! default value
     472      call getin("slope_orientation",psi_sl(1))
     473     
     474c
     475c  for the gravity wave scheme
    460476c  ---------------------------------
    461477c
    … …  
    467483
    468484
    469 
    470 c   Initialisation speciales "physiq"
    471 c   ---------------------------------
    472 c   la surface de chaque maille est inutile en 1D --->
     485c   Specific initializations for "physiq"
     486c   -------------------------------------
     487c   mesh surface (not a very usefull quantity in 1D)
    473488      area(1)=1.E+0
    474489
    475 c   le geopotentiel au sol est inutile en 1D car tout est controle
    476 c   par la pression de surface --->
     490c   surface geopotential is not used (or useful) since in 1D
     491c   everything is controled by surface pressure
    477492      phisfi(1)=0.E+0
    478493
    479 c  "inifis" reproduit un certain nombre d'initialisations deja faites
    480 c  + lecture des clefs de callphys.def
    481 c  + calcul de la frequence d'appel au rayonnement
    482 c  + calcul des sinus et cosinus des longitude latitude
     494c  "inifis" does some initializations (some of which have already been
     495c  done above!) and loads parameters set in callphys.def
    483496
    484497!Mars possible matter with dtphys in input and include!!!
    485498      CALL inifis(1,llm,day0,daysec,dtphys,
    486499     .            lati,long,area,rad,g,r,cpp)
    487 c   Initialisation pour prendre en compte les vents en 1-D
     500
     501c   Initialization to take into account prescribed winds
    488502c   ------------------------------------------------------
    489503      ptif=2.E+0*omeg*sinlat(1)
    490504 
    491 c    vent geostrophique
     505c    geostrophic wind
    492506      gru=10. ! default value for gru
    493507      PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
    … …  
    500514      write(*,*) " v = ",grv
    501515
    502 c     Initialisation des vents  au premier pas de temps
     516c     Initialize winds  for first time step
    503517      DO ilayer=1,nlayer
    504518         u(ilayer)=gru
    … …  
    506520      ENDDO
    507521
    508 c     energie cinetique turbulente
     522c     Initialize turbulente kinetic energy
    509523      DO ilevel=1,nlevel
    510524         q2(ilevel)=0.E+0
    511525      ENDDO
    512526
    513 c  glace de CO2 au sol
     527c  CO2 ice on the surface
    514528c  -------------------
    515529      co2ice=0.E+0 ! default value for co2ice
    … …  
    519533
    520534c
    521 c  emissivite
     535c  emissivity
    522536c  ----------
    523537      emis=emissiv
    … …  
    529543 
    530544
    531 c  calcul des pressions et altitudes en utilisant les niveaux sigma
     545c  Compute pressures and altitudes of atmospheric levels
    532546c  ----------------------------------------------------------------
    533547
    … …  
    555569
    556570
    557 c  profil de temperature au premier appel
     571c  Initialize temperature profile
    558572c  --------------------------------------
    559573      pks=psurf**rcp
    560574
    561 c altitude en km dans profile: on divise zlay par 1000
     575c altitude in km in profile: divide zlay by 1000
    562576      tmp1(0)=0.E+0
    563577      DO ilayer=1,nlayer
    … …  
    591605      enddo
    592606
    593 c    Initialisation des traceurs
     607c    Initialize traceurs
    594608c    ---------------------------
    595609
    596610      if (photochem.or.callthermos) then
    597          write(*,*) 'Especes chimiques initialisees'
    598          ! thermo=0: initialisation pour toutes les couches
     611         write(*,*) 'Initializing chemical species'
     612         ! thermo=0: initialize over all atmospheric layers
    599613         thermo=0
    600614         call inichim_newstart(q,psurf,sig,nqmx,lati,long,area,
    … …  
    602616      endif
    603617
    604 c Regarder si le sol est un reservoir de glace d'eau
     618c Check if the surface is a water ice reservoir
    605619c --------------------------------------------------
    606       watercaptag(ngridmx)=.false. ! Par defaut il n'y pas de glace au sol
     620      watercaptag(ngridmx)=.false. ! Default: no water ice reservoir
    607621      print *,'Water ice cap on ground ?'
    608622      call getin("watercaptag",watercaptag)
    … …  
    610624     
    611625
    612 c    Initialisation pour les sorties GRADS dans "g1d.dat" et "g1d.ctl"
     626c    Initialization for GRADS outputs in "g1d.dat" and "g1d.ctl"
    613627c    ----------------------------------------------------------------
    614 c    (effectuee avec "writeg1d" a partir de "physiq.F"
    615 c    ou tout autre subroutine physique, y compris celle ci).
     628c    (output done in "writeg1d", typically called by "physiq.F")
    616629
    617630        g1d_nlayer=nlayer
    … …  
    623636        g2d_premier=.true.
    624637
    625 c  Ecriture de "startfi"
     638c  Write a "startfi" file
    626639c  --------------------
    627 c  (Ce fichier sera aussitot relu au premier
    628 c   appel de "physiq", mais il est necessaire pour passer
    629 c   les variables purement physiques a "physiq"...
     640c  This file will be read during the first call to "physiq".
     641c  It is needed to transfert physics variables to "physiq"...
    630642
    631643      call physdem1("startfi.nc",long,lati,nsoilmx,nqmx,
    … …  
    633645     .              tsoil,co2ice,emis,q2,qsurf,area,albedodat,
    634646     .              inertiedat,zmea,zstd,zsig,zgam,zthe)
     647
    635648c=======================================================================
    636 c  BOUCLE TEMPORELLE DU MODELE 1D
     649c  1D MODEL TIME STEPPING LOOP
    637650c=======================================================================
    638651c
    … …  
    656669        ENDIF
    657670
    658 c    calcul du geopotentiel
     671c     compute geopotential
    659672c     ~~~~~~~~~~~~~~~~~~~~~
    660673      DO ilayer=1,nlayer
    … …  
    670683      ENDDO
    671684
    672 c       appel de la physique
     685c       call physics
    673686c       --------------------
    674687!      write(*,*) "testphys1d avant q", q(1,:)
    … …  
    679692     ,     u, v,temp, q, 
    680693     ,     w,
    681 C - sorties
     694C - outputs
    682695     s     du, dv, dtemp, dq,dpsurf,tracerdyn)
    683696!      write(*,*) "testphys1d apres q", q(1,:)
    684697
    685698
    686 c       evolution du vent : modele 1D
    687 c       -----------------------------
    688  
    689 c       la physique calcule les derivees temporelles de u et v.
    690 c       on y rajoute betement un effet Coriolis.
     699c       wind increment : specific for 1D
     700c       --------------------------------
     701 
     702c       The physics compute the tendencies on u and v,
     703c       here we just add Coriolos effect
    691704c
    692705c       DO ilayer=1,nlayer
    … …  
    695708c       ENDDO
    696709
    697 c       Pour certain test : pas de coriolis a l'equateur
     710c       For some tests : No coriolis force at equator
    698711c       if(lati(1).eq.0.) then
    699712          DO ilayer=1,nlayer
    … …  
    704717c     
    705718c
    706 c       Calcul du temps au pas de temps suivant
     719c       Compute time for next time step
    707720c       ---------------------------------------
    708721        firstcall=.false.
    … …  
    713726        ENDIF
    714727
    715 c       calcul des vitesses et temperature au pas de temps suivant
     728c       compute winds and temperature for next time step
    716729c       ----------------------------------------------------------
    717730
    … …  
    722735        ENDDO
    723736
    724 c       calcul des pressions au pas de temps suivant
     737c       compute pressure for next time step
    725738c       ----------------------------------------------------------
    726739
    727            psurf=psurf+dtphys*dpsurf   ! evolution de la pression de surface
     740           psurf=psurf+dtphys*dpsurf   ! surface pressure change
    728741           DO ilevel=1,nlevel
    729742             plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
    … …  
    733746           ENDDO
    734747
    735 !       increment tracer values
     748!       increment tracers
    736749        DO iq = 1, nqmx
    737750          DO ilayer=1,nlayer
    … …  
    740753        ENDDO
    741754
    742       ENDDO   ! fin de la boucle temporelle
     755      ENDDO   ! of idt=1,ndt ! end of time stepping loop
    743756
    744757c    ========================================================
    745 c    GESTION DES SORTIE
     758c    OUTPUTS
    746759c    ========================================================
    747760
    748 c    fermeture pour conclure les sorties format grads dans "g1d.dat"
    749 c    et "g1d.ctl" (effectuee avec "writeg1d" a partir de "physiq.F"
    750 c    ou tout autre subroutine physique, y compris celle ci).
     761c    finalize and close grads files "g1d.dat" and "g1d.ctl"
    751762
    752763c        CALL endg1d(1,nlayer,zphi/(g*1000.),ndt)
    … …  
    758769c***********************************************************************
    759770c***********************************************************************
    760 c     Subroutines Bidons utilise seulement en 3D, mais
    761 c     necessaire a la compilation de testphys1d en 1-D
    762 
    763       subroutine gr_fi_dyn
    764       RETURN
    765       END
     771c     Dummy subroutines used only in 3D, but required to
     772c     compile testphys1d (to cleanly use writediagfi)
     773
     774      subroutine gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
     775
     776      IMPLICIT NONE
     777
     778      INTEGER im,jm,ngrid,nfield
     779      REAL pdyn(im,jm,nfield)
     780      REAL pfi(ngrid,nfield)
     781     
     782      if (ngrid.ne.1) then
     783        write(*,*) "gr_fi_dyn error: in 1D ngrid should be 1!!!"
     784        stop
     785      endif
     786     
     787      pdyn(1,1,1:nfield)=pfi(1,1:nfield)
     788     
     789      end
    766790 
    767791c***********************************************************************
Note: See TracChangeset for help on using the changeset viewer.