Changeset 1671 for LMDZ5


Ignore:
Timestamp:
Oct 24, 2012, 9:10:10 AM (12 years ago)
Author:
Ehouarn Millour
Message:
  • fixed "aquaplanet case" so that initializations (creation of files startphy.nc and limit.nc) now also works in parallel (mpi,omp,mixed).
  • call to "iniaqua" is now done from within "iniphysiq" ; also added some tests to check consistency between essential variables shared by dynamics and physics (planetary radius, gravity, Cp, ...)
  • simillarily adapted "phydev" routines, and added necessary routines to also be able to generate startphy.nc/restartphy.nc files there. Also removed common file "comcstphy.h" and replaced it with a module "comcstphy.F90"

EM

Location:
LMDZ5/trunk/libf
Files:
4 added
1 deleted
10 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/dyn3d/comconst.h

    r1505 r1671  
    2121      REAL dtdiss ! (s) time step for the dissipation
    2222      REAL rad ! (m) radius of the planet
    23       REAL r ! Gas constant R=8.31 J.K-1.mol-1
    24       REAL cpp   ! Cp
     23      REAL r ! Reduced Gas constant r=R/mu
     24             ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
     25      REAL cpp   ! Specific heat Cp (J.kg-1.K-1)
    2526      REAL kappa ! kappa=R/Cp
    2627      REAL cotot

  • LMDZ5/trunk/libf/dyn3d/gcm.F

    r1615 r1671  
    433433! Physics:
    434434#ifdef CPP_PHYS
    435          CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys ,
    436      ,                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp     )
     435         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys,
     436     &                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp,
     437     &                iflag_phys)
    437438#endif
    438439         call_iniphys=.false.
    … …  
    457458 301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
    458459 302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
    459 #endif
    460 
    461 #ifdef CPP_PHYS
    462 ! Create start file (startphy.nc) and boundary conditions (limit.nc)
    463 ! for the Earth verstion
    464        if (iflag_phys>=100) then
    465           call iniaqua(ngridmx,latfi,lonfi,iflag_phys)
    466        endif
    467460#endif
    468461

  • LMDZ5/trunk/libf/dyn3dpar/comconst.h

    r1505 r1671  
    2121      REAL dtdiss ! (s) time step for the dissipation
    2222      REAL rad ! (m) radius of the planet
    23       REAL r ! Gas constant R=8.31 J.K-1.mol-1
    24       REAL cpp   ! Cp
     23      REAL r ! Reduced Gas constant r=R/mu
     24             ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
     25      REAL cpp   ! Specific heat Cp (J.kg-1.K-1)
    2526      REAL kappa ! kappa=R/Cp
    2627      REAL cotot

  • LMDZ5/trunk/libf/dyn3dpar/gcm.F

    r1615 r1671  
    423423c   Initialisation de la physique :
    424424c   -------------------------------
    425       IF (call_iniphys.and.iflag_phys.eq.1) THEN
     425      IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN
    426426         latfi(1)=rlatu(1)
    427427         lonfi(1)=0.
    … …  
    446446! Physics:
    447447#ifdef CPP_PHYS
    448          CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys ,
    449      ,                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp     )
     448         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys,
     449     &                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp,
     450     &                iflag_phys)
    450451#endif
    451452         call_iniphys=.false.
    452       ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
     453      ENDIF ! of IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100))
    453454
    454455
    … …  
    481482 301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
    482483 302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
    483 #endif
    484 
    485 #ifdef CPP_PHYS
    486 ! Create start file (startphy.nc) and boundary conditions (limit.nc)
    487 ! for the Earth verstion
    488        if (iflag_phys>=100) then
    489           call iniaqua(ngridmx,latfi,lonfi,iflag_phys)
    490        endif
    491484#endif
    492485

  • LMDZ5/trunk/libf/dyn3dpar/leapfrog_p.F

    r1626 r1671  
    139139      REAL :: secondes
    140140
     141      logical :: physic
    141142      LOGICAL first,callinigrads
    142143
    … …  
    208209
    209210      itau = 0
     211      physic=.true.
     212      if (iflag_phys==0.or.iflag_phys==2) physic=.false.
    210213!      iday = day_ini+itau/day_step
    211214!      time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
    … …  
    364367     s        apdiss = .TRUE.
    365368         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
    366      s          .and. iflag_phys.EQ.1                 ) apphys = .TRUE.
     369     s          .and. physic                        ) apphys = .TRUE.
    367370      ELSE
    368371      ! Leapfrog/Matsuno time stepping
    … …  
    370373         IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
    371374     s        apdiss = .TRUE.
    372          IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
     375         IF( MOD(itau+1,iphysiq).EQ.0.AND.physic) apphys=.TRUE.
    373376      END IF
    374377

  • LMDZ5/trunk/libf/phydev/iniphysiq.F

    r1615 r1671  
    22! $Id: iniphysiq.F 1403 2010-07-01 09:02:53Z fairhead $
    33!
    4 c
    5 c
    64      SUBROUTINE iniphysiq(ngrid,nlayer,
    75     $           punjours,
    86     $           pdayref,ptimestep,
    97     $           plat,plon,parea,pcu,pcv,
    10      $           prad,pg,pr,pcpp)
    11       USE dimphy
    12       USE mod_grid_phy_lmdz
    13       USE mod_phys_lmdz_para
    14       USE comgeomphy
     8     $           prad,pg,pr,pcpp,iflag_phys)
     9      USE dimphy, only : klev
     10      USE mod_grid_phy_lmdz, only : klon_glo
     11      USE mod_phys_lmdz_para, only : klon_omp,klon_omp_begin,
     12     &                               klon_omp_end,klon_mpi_begin
     13      USE comgeomphy, only : airephy,cuphy,cvphy,rlond,rlatd
     14      USE comcstphy, only : rradius,rg,rr,rcpp
    1515
    1616      IMPLICIT NONE
    … …  
    1818c=======================================================================
    1919c
    20 c   subject:
    21 c   --------
     20c   Initialisation of the physical constants and some positional and
     21c   geometrical arrays for the physics
    2222c
    23 c   Initialisation for the physical parametrisations of the LMD
    24 c   martian atmospheric general circulation modele.
    25 c
    26 c   author: Frederic Hourdin 15 / 10 /93
    27 c   -------
    28 c
    29 c   arguments:
    30 c   ----------
    31 c
    32 c   input:
    33 c   ------
    3423c
    3524c    ngrid                 Size of the horizontal grid.
    … …  
    3726c    nlayer                Number of vertical layers.
    3827c    pdayref               Day of reference for the simulation
    39 c    firstcall             True at the first call
    40 c    lastcall              True at the last call
    41 c    pday                  Number of days counted from the North. Spring
    42 c                          equinoxe.
    4328c
    4429c=======================================================================
    45 c
    46 c-----------------------------------------------------------------------
    47 c   declarations:
    48 c   -------------
     30 
    4931 
    5032cym#include "dimensions.h"
    5133cym#include "dimphy.h"
    5234cym#include "comgeomphy.h"
    53 #include "comcstphy.h"
    54       REAL prad,pg,pr,pcpp,punjours
    55  
    56       INTEGER ngrid,nlayer
    57       REAL plat(ngrid),plon(ngrid),parea(klon_glo)
    58       REAL pcu(klon_glo),pcv(klon_glo)
    59       INTEGER pdayref
     35#include "iniprint.h"
     36
     37      REAL,INTENT(IN) :: prad ! radius of the planet (m)
     38      REAL,INTENT(IN) :: pg ! gravitational acceleration (m/s2)
     39      REAL,INTENT(IN) :: pr ! ! reduced gas constant R/mu
     40      REAL,INTENT(IN) :: pcpp ! specific heat Cp
     41      REAL,INTENT(IN) :: punjours ! length (in s) of a standard day
     42      INTEGER,INTENT(IN) :: ngrid ! number of horizontal grid points in the physics
     43      INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
     44      REAL,INTENT(IN) :: plat(ngrid) ! latitudes of the physics grid
     45      REAL,INTENT(IN) :: plon(ngrid) ! longitudes of the physics grid
     46      REAL,INTENT(IN) :: parea(klon_glo) ! area (m2)
     47      REAL,INTENT(IN) :: pcu(klon_glo) ! cu coeff. (u_covariant = cu * u)
     48      REAL,INTENT(IN) :: pcv(klon_glo) ! cv coeff. (v_covariant = cv * v)
     49      INTEGER,INTENT(IN) :: pdayref ! reference day of for the simulation
     50      REAL,INTENT(IN) :: ptimestep !physics time step (s)
     51      INTEGER,INTENT(IN) :: iflag_phys ! type of physics to be called
     52
    6053      INTEGER :: ibegin,iend,offset
    61  
    62       REAL ptimestep
    6354      CHARACTER (LEN=20) :: modname='iniphysiq'
    6455      CHARACTER (LEN=80) :: abort_message
    6556 
    6657      IF (nlayer.NE.klev) THEN
    67          PRINT*,'STOP in inifis'
    68          PRINT*,'Probleme de dimensions :'
    69          PRINT*,'nlayer     = ',nlayer
    70          PRINT*,'klev   = ',klev
     58         write(lunout,*) 'STOP in ',trim(modname)
     59         write(lunout,*) 'Problem with dimensions :'
     60         write(lunout,*) 'nlayer     = ',nlayer
     61         write(lunout,*) 'klev   = ',klev
    7162         abort_message = ''
    7263         CALL abort_gcm (modname,abort_message,1)
    … …  
    7465
    7566      IF (ngrid.NE.klon_glo) THEN
    76          PRINT*,'STOP in inifis'
    77          PRINT*,'Probleme de dimensions :'
    78          PRINT*,'ngrid     = ',ngrid
    79          PRINT*,'klon   = ',klon_glo
     67         write(lunout,*) 'STOP in ',trim(modname)
     68         write(lunout,*) 'Problem with dimensions :'
     69         write(lunout,*) 'ngrid     = ',ngrid
     70         write(lunout,*) 'klon   = ',klon_glo
    8071         abort_message = ''
    8172         CALL abort_gcm (modname,abort_message,1)
    8273      ENDIF
    83 c$OMP PARALLEL PRIVATE(ibegin,iend)
    84 c$OMP+         SHARED(parea,pcu,pcv,plon,plat)
     74
     75!$OMP PARALLEL PRIVATE(ibegin,iend)
     76!$OMP+         SHARED(parea,pcu,pcv,plon,plat)
    8577     
    8678      offset=klon_mpi_begin-1
    … …  
    9284      rlatd(1:klon_omp)=plat(offset+klon_omp_begin:offset+klon_omp_end)
    9385
    94 !     call suphel
    95 !     prad,pg,pr,pcpp
     86! copy some fundamental parameters to physics
    9687      rradius=prad
    9788      rg=pg
    … …  
    9990      rcpp=pcpp
    10091
    101      
     92!$OMP END PARALLEL
    10293
    103 c$OMP END PARALLEL
     94!      print*,'ATTENTION !!! TRAVAILLER SUR INIPHYSIQ'
     95!      print*,'CONTROLE DES LATITUDES, LONGITUDES, PARAMETRES ...'
    10496
    105       print*,'ATTENTION !!! TRAVAILLER SUR INIPHYSIQ'
    106       print*,'CONTROLE DES LATITUDES, LONGITUDES, PARAMETRES ...'
     97! Additional initializations for aquaplanets
     98!$OMP PARALLEL
     99      if (iflag_phys>=100) then
     100        call iniaqua(klon_omp,rlatd,rlond,iflag_phys)
     101      endif
     102!$OMP END PARALLEL
    107103
    108       RETURN
    109 9999  CONTINUE
    110       abort_message ='Cette version demande les fichier rnatur.dat
    111      & et surf.def'
    112       CALL abort_gcm (modname,abort_message,1)
     104!      RETURN
     105!9999  CONTINUE
     106!      abort_message ='Cette version demande les fichier rnatur.dat
     107!     & et surf.def'
     108!      CALL abort_gcm (modname,abort_message,1)
    113109
    114110      END

  • LMDZ5/trunk/libf/phydev/phyaqua.F

    r1615 r1671  
    1 ! Routines complementaires pour la physique planetaire.
    2 
     1!
     2! $Id: $
     3!
    34
    45      subroutine iniaqua(nlon,latfi,lonfi,iflag_phys)
    56
    67!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    7 !  Creation d'un etat initial et de conditions aux limites
    8 !  (resp startphy.nc et limit.nc) pour des configurations idealisees
    9 ! du modele LMDZ dans sa version terrestre.
    10 !  iflag_phys est un parametre qui controle
    11 !  iflag_phys = N 
    12 !    de 100 a 199 : aqua planetes avec SST forcees
    13 !                 N-100 determine le type de SSTs
    14 !    de 200 a 299 : terra planetes avec Ts calcule
    15 !       
     8!  Create an initial state (startphy.nc) for the physics
     9!  Usefull for idealised cases (e.g. aquaplanets or testcases)
    1610!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1711
     12      use phys_state_var_mod, only : rlat,rlon,
     13     &                               phys_state_var_init
     14      use mod_phys_lmdz_para, only : klon_omp
     15      use comgeomphy, only : rlond,rlatd
     16      implicit none
     17     
     18      integer,intent(in) :: nlon,iflag_phys
     19      real,intent(in) :: lonfi(nlon),latfi(nlon)
    1820
    19       integer nlon,iflag_phys
    20 cIM ajout latfi, lonfi
    21       REAL, DIMENSION (nlon) :: lonfi, latfi
     21! local variables
     22      real :: pi
     23
     24! initializations:
     25      pi=2.*asin(1.)
     26
     27      call phys_state_var_init()
     28
     29      rlat(1:klon_omp)=rlatd(1:klon_omp)*180./pi
     30      rlon(1:klon_omp)=rlond(1:klon_omp)*180./pi
    2231
    2332
    24       return
     33! Here you could create an initial condition for the physics
     34! ...
     35! ... fill in the fields...
     36! ...
     37! ... and create a "startphy.nc" file
     38!      CALL phyredem ("startphy.nc")
     39
    2540      end
    2641

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

    r1615 r1671  
    1111     &            , PVteta)
    1212
    13       USE dimphy
    14       USE infotrac
    15       USE comgeomphy
     13      USE dimphy, only : klon,klev
     14      USE infotrac, only : nqtot
     15      USE comgeomphy, only : rlatd
     16      USE comcstphy, only : rg
    1617
    1718      IMPLICIT none
    … …  
    5051!======================================================================
    5152#include "dimensions.h"
    52 #include "comcstphy.h"
     53!#include "comcstphy.h"
    5354
    5455      integer jjmp1
    … …  
    101102    PARAMETER    ( longcles = 20 )
    102103
    103 real temp_newton(klon,klev)
     104real :: temp_newton(klon,klev)
    104105integer k
    105106logical, save :: first=.true.
     107!$OMP THREADPRIVATE(first)
    106108
    107109      REAL clesphy0( longcles      )
    108110
    109 d_u=0.
    110 d_v=0.
    111 d_t=0.
    112 d_qx=0.
    113 d_ps=0.
     111! initializations
     112if (first) then
     113! ...
    114114
    115      d_u(:,1)=-u(:,1)/86400.
    116      do k=1,klev
    117         temp_newton(:,k)=280.+cos(rlatd(:))*40.-pphi(:,k)/rg*6.e-3
    118         d_t(:,k)=(temp_newton(:,k)-t(:,k))/1.e5
    119      enddo
     115  first=.false.
     116endif
     117
     118! set all tendencies to zero
     119d_u(:,:)=0.
     120d_v(:,:)=0.
     121d_t(:,:)=0.
     122d_qx(:,:,:)=0.
     123d_ps(:)=0.
     124
     125! compute tendencies to return to the dynamics:
     126! "friction" on the first layer
     127d_u(:,1)=-u(:,1)/86400.
     128! newtonian rlaxation towards temp_newton()
     129do k=1,klev
     130  temp_newton(:,k)=280.+cos(rlatd(:))*40.-pphi(:,k)/rg*6.e-3
     131  d_t(:,k)=(temp_newton(:,k)-t(:,k))/1.e5
     132enddo
    120133
    121134
    122       print*,'COUCOU PHYDEV'
    123       return
    124       end
     135print*,'COUCOU PHYDEV'
     136
     137! if lastcall, then it is time to write "restartphy.nc" file
     138if (lafin) then
     139  call phyredem("restartphy.nc")
     140endif
     141
     142end

  • LMDZ5/trunk/libf/phylmd/iniphysiq.F

    r1403 r1671  
    88     $           pdayref,ptimestep,
    99     $           plat,plon,parea,pcu,pcv,
    10      $           prad,pg,pr,pcpp)
    11       USE dimphy
    12       USE mod_grid_phy_lmdz
    13       USE mod_phys_lmdz_para
    14       USE comgeomphy
     10     $           prad,pg,pr,pcpp,iflag_phys)
     11      USE dimphy, only : klev
     12      USE mod_grid_phy_lmdz, only : klon_glo
     13      USE mod_phys_lmdz_para, only : klon_omp,klon_omp_begin,
     14     &                               klon_omp_end,klon_mpi_begin
     15      USE comgeomphy, only : airephy,cuphy,cvphy,rlond,rlatd
    1516
    1617      IMPLICIT NONE
    … …  
    1819c=======================================================================
    1920c
    20 c   subject:
    21 c   --------
     21c   Initialisation of the physical constants and some positional and
     22c   geometrical arrays for the physics
    2223c
    23 c   Initialisation for the physical parametrisations of the LMD
    24 c   martian atmospheric general circulation modele.
    25 c
    26 c   author: Frederic Hourdin 15 / 10 /93
    27 c   -------
    28 c
    29 c   arguments:
    30 c   ----------
    31 c
    32 c   input:
    33 c   ------
    3424c
    3525c    ngrid                 Size of the horizontal grid.
    … …  
    3727c    nlayer                Number of vertical layers.
    3828c    pdayref               Day of reference for the simulation
    39 c    firstcall             True at the first call
    40 c    lastcall              True at the last call
    41 c    pday                  Number of days counted from the North. Spring
    42 c                          equinoxe.
    4329c
    4430c=======================================================================
    45 c
    46 c-----------------------------------------------------------------------
    47 c   declarations:
    48 c   -------------
    4931 
    5032cym#include "dimensions.h"
    … …  
    5234cym#include "comgeomphy.h"
    5335#include "YOMCST.h"
    54       REAL prad,pg,pr,pcpp,punjours
    55  
    56       INTEGER ngrid,nlayer
    57       REAL plat(ngrid),plon(ngrid),parea(klon_glo)
    58       REAL pcu(klon_glo),pcv(klon_glo)
    59       INTEGER pdayref
    60       INTEGER :: ibegin,iend,offset
    61  
    62       REAL ptimestep
     36#include "iniprint.h"
     37
     38      REAL,INTENT(IN) :: prad ! radius of the planet (m)
     39      REAL,INTENT(IN) :: pg ! gravitational acceleration (m/s2)
     40      REAL,INTENT(IN) :: pr ! ! reduced gas constant R/mu
     41      REAL,INTENT(IN) :: pcpp ! specific heat Cp
     42      REAL,INTENT(IN) :: punjours ! length (in s) of a standard day
     43      INTEGER,INTENT(IN) :: ngrid ! number of horizontal grid points in the physics
     44      INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
     45      REAL,INTENT(IN) :: plat(ngrid) ! latitudes of the physics grid
     46      REAL,INTENT(IN) :: plon(ngrid) ! longitudes of the physics grid
     47      REAL,INTENT(IN) :: parea(klon_glo) ! area (m2)
     48      REAL,INTENT(IN) :: pcu(klon_glo) ! cu coeff. (u_covariant = cu * u)
     49      REAL,INTENT(IN) :: pcv(klon_glo) ! cv coeff. (v_covariant = cv * v)
     50      INTEGER,INTENT(IN) :: pdayref ! reference day of for the simulation
     51      REAL,INTENT(IN) :: ptimestep !physics time step (s)
     52      INTEGER,INTENT(IN) :: iflag_phys ! type of physics to be called
     53
     54      INTEGER :: ibegin,iend,offset
    6355      CHARACTER (LEN=20) :: modname='iniphysiq'
    6456      CHARACTER (LEN=80) :: abort_message
    6557 
    6658      IF (nlayer.NE.klev) THEN
    67          PRINT*,'STOP in inifis'
    68          PRINT*,'Probleme de dimensions :'
    69          PRINT*,'nlayer     = ',nlayer
    70          PRINT*,'klev   = ',klev
     59         write(lunout,*) 'STOP in ',trim(modname)
     60         write(lunout,*) 'Problem with dimensions :'
     61         write(lunout,*) 'nlayer     = ',nlayer
     62         write(lunout,*) 'klev   = ',klev
    7163         abort_message = ''
    7264         CALL abort_gcm (modname,abort_message,1)
    … …  
    7466
    7567      IF (ngrid.NE.klon_glo) THEN
    76          PRINT*,'STOP in inifis'
    77          PRINT*,'Probleme de dimensions :'
    78          PRINT*,'ngrid     = ',ngrid
    79          PRINT*,'klon   = ',klon_glo
     68         write(lunout,*) 'STOP in ',trim(modname)
     69         write(lunout,*) 'Problem with dimensions :'
     70         write(lunout,*) 'ngrid     = ',ngrid
     71         write(lunout,*) 'klon   = ',klon_glo
    8072         abort_message = ''
    8173         CALL abort_gcm (modname,abort_message,1)
    8274      ENDIF
    83 c$OMP PARALLEL PRIVATE(ibegin,iend)
    84 c$OMP+         SHARED(parea,pcu,pcv,plon,plat)
     75
     76!$OMP PARALLEL PRIVATE(ibegin,iend)
     77!$OMP+         SHARED(parea,pcu,pcv,plon,plat)
    8578     
    8679      offset=klon_mpi_begin-1
    … …  
    9285      rlatd(1:klon_omp)=plat(offset+klon_omp_begin:offset+klon_omp_end)
    9386
     87      ! suphel => initialize some physical constants (orbital parameters,
     88      !           geoid, gravity, thermodynamical constants, etc.) in the
     89      !           physics
    9490      call suphel
     91     
     92!$OMP END PARALLEL
    9593
    96 c$OMP END PARALLEL
     94      ! check that physical constants set in 'suphel' are coherent
     95      ! with values set in the dynamics:
     96      if (RDAY.ne.punjours) then
     97        write(lunout,*) "iniphysiq: length of day discrepancy!!!"
     98        write(lunout,*) "  in the dynamics punjours=",punjours
     99        write(lunout,*) "   but in the physics RDAY=",RDAY
     100        if (abs(RDAY-punjours).gt.0.01) then
     101          ! stop here if the relative difference is more than 1%
     102          abort_message = 'length of day discrepancy'
     103          CALL abort_gcm (modname,abort_message,1)
     104        endif
     105      endif
     106      if (RG.ne.pg) then
     107        write(lunout,*) "iniphysiq: gravity discrepancy !!!"
     108        write(lunout,*) "     in the dynamics pg=",pg
     109        write(lunout,*) "  but in the physics RG=",RG
     110        if (abs(RG-pg).gt.0.01) then
     111          ! stop here if the relative difference is more than 1%
     112          abort_message = 'gravity discrepancy'
     113          CALL abort_gcm (modname,abort_message,1)
     114        endif
     115      endif
     116      if (RA.ne.prad) then
     117        write(lunout,*) "iniphysiq: planet radius discrepancy !!!"
     118        write(lunout,*) "   in the dynamics prad=",prad
     119        write(lunout,*) "  but in the physics RA=",RA
     120        if (abs(RA-prad).gt.0.01) then
     121          ! stop here if the relative difference is more than 1%
     122          abort_message = 'planet radius discrepancy'
     123          CALL abort_gcm (modname,abort_message,1)
     124        endif
     125      endif
     126      if (RD.ne.pr) then
     127        write(lunout,*)"iniphysiq: reduced gas constant discrepancy !!!"
     128        write(lunout,*)"     in the dynamics pr=",pr
     129        write(lunout,*)"  but in the physics RD=",RD
     130        if (abs(RD-pr).gt.0.01) then
     131          ! stop here if the relative difference is more than 1%
     132          abort_message = 'reduced gas constant discrepancy'
     133          CALL abort_gcm (modname,abort_message,1)
     134        endif
     135      endif
     136      if (RCPD.ne.pcpp) then
     137        write(lunout,*)"iniphysiq: specific heat discrepancy !!!"
     138        write(lunout,*)"     in the dynamics pcpp=",pcpp
     139        write(lunout,*)"  but in the physics RCPD=",RCPD
     140        if (abs(RCPD-pcpp).gt.0.01) then
     141          ! stop here if the relative difference is more than 1%
     142          abort_message = 'specific heat discrepancy'
     143          CALL abort_gcm (modname,abort_message,1)
     144        endif
     145      endif
    97146
    98       print*,'ATTENTION !!! TRAVAILLER SUR INIPHYSIQ'
    99       print*,'CONTROLE DES LATITUDES, LONGITUDES, PARAMETRES ...'
     147! Additional initializations for aquaplanets
     148!$OMP PARALLEL
     149      if (iflag_phys>=100) then
     150        call iniaqua(klon_omp,rlatd,rlond,iflag_phys)
     151      endif
     152!$OMP END PARALLEL
    100153
    101       RETURN
    102 9999  CONTINUE
    103       abort_message ='Cette version demande les fichier rnatur.dat
    104      & et surf.def'
    105       CALL abort_gcm (modname,abort_message,1)
     154!      RETURN
     155!9999  CONTINUE
     156!      abort_message ='Cette version demande les fichier rnatur.dat
     157!     & et surf.def'
     158!      CALL abort_gcm (modname,abort_message,1)
    106159
    107160      END

  • LMDZ5/trunk/libf/phylmd/phyaqua.F

    r1530 r1671  
    1616!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1717
    18       use comgeomphy
    19       use dimphy
     18      use comgeomphy, only : rlatd,rlond
     19      use dimphy, only : klon
    2020      use surface_data, only : type_ocean,ok_veget
    2121      use pbl_surface_mod, only : pbl_surface_init
    2222      USE fonte_neige_mod, only : fonte_neige_init
    2323      use phys_state_var_mod
    24       use control_mod
    25 
     24      use control_mod, only : dayref,nday,iphysiq
    2625
    2726      USE IOIPSL
    … …  
    3534#include "dimsoil.h"
    3635#include "indicesol.h"
    37 
    38       integer nlon,iflag_phys
     36#include "temps.h"
     37
     38      integer,intent(in) :: nlon,iflag_phys
    3939cIM ajout latfi, lonfi
    40       REAL, DIMENSION (nlon) :: lonfi, latfi
     40      real,intent(in) :: lonfi(nlon),latfi(nlon)
     41
    4142      INTEGER type_profil,type_aqua
    4243
    … …  
    7172!      integer demih_pas
    7273
    73       integer day_ini
    74 
    7574      CHARACTER*80 ans,file_forctl, file_fordat, file_start
    7675      character*100 file,var
    … …  
    8887      REAL phy_flic(nlon,360)
    8988
    90       integer, save::  read_climoz ! read ozone climatology
     89      integer, save::  read_climoz=0 ! read ozone climatology
    9190
    9291
    … …  
    131130      type_aqua=iflag_phys/100
    132131      type_profil=iflag_phys-type_aqua*100
    133       print*,'type_aqua, type_profil',type_aqua, type_profil
    134 
    135       if (klon.ne.nlon) stop'probleme de dimensions dans iniaqua'
     132      print*,'iniaqua:type_aqua, type_profil',type_aqua, type_profil
     133
     134      if (klon.ne.nlon) then
     135        write(*,*)"iniaqua: klon=",klon," nlon=",nlon
     136        stop'probleme de dimensions dans iniaqua'
     137      endif
    136138      call phys_state_var_init(read_climoz)
    137139
    … …  
    154156
    155157         day_ini=dayref
     158         day_end=day_ini+nday
    156159         airefi=1.
    157160         zcufi=1.
    … …  
    171174      radsol=0.
    172175      qsol_f=10.
    173       CALL getin('albedo',albedo)
     176!      CALL getin('albedo',albedo) ! albedo is set below, depending on type_aqua
    174177      alb_ocean=.true.
    175178      CALL getin('alb_ocean',alb_ocean)
    … …  
    180183      qsol(:)    = qsol_f
    181184      rugsrel = 0.0    ! (rugsrel = rugoro)
     185      rugoro = 0.0
     186      u_ancien = 0.0
     187      v_ancien = 0.0
    182188      agesno  = 50.0
    183189! Relief plat
    … …  
    308314     .     evap, frugs, agesno, tsoil)
    309315
    310         print*,'avant phyredem dans iniaqua'
     316        print*,'iniaqua: before phyredem'
    311317
    312318      falb1=albedo
    … …  
    329335      CALL phyredem ("startphy.nc")
    330336
    331         print*,'apres phyredem'
     337        print*,'iniaqua: after phyredem'
    332338      call phys_state_var_end
    333339
    … …  
    450456      RETURN
    451457      END
     458
     459!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     460
    452461      subroutine writelim
    453462     s   (klon,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
    454463     s    phy_fter,phy_foce,phy_flic,phy_fsic)
    455464c
     465      use mod_phys_lmdz_para, only: is_mpi_root,is_omp_root
     466      use mod_grid_phy_lmdz, only : klon_glo
     467      use mod_phys_lmdz_transfert_para, only : gather
    456468!#include "dimensions.h"
    457469!#include "dimphy.h"
    458470#include "netcdf.inc"
    459471 
    460       integer klon
    461       REAL phy_nat(klon,360)
    462       REAL phy_alb(klon,360)
    463       REAL phy_sst(klon,360)
    464       REAL phy_bil(klon,360)
    465       REAL phy_rug(klon,360)
    466       REAL phy_ice(klon,360)
    467       REAL phy_fter(klon,360)
    468       REAL phy_foce(klon,360)
    469       REAL phy_flic(klon,360)
    470       REAL phy_fsic(klon,360)
    471  
     472      integer,intent(in) :: klon
     473      real,intent(in) :: phy_nat(klon,360)
     474      real,intent(in) :: phy_alb(klon,360)
     475      real,intent(in) :: phy_sst(klon,360)
     476      real,intent(in) :: phy_bil(klon,360)
     477      real,intent(in) :: phy_rug(klon,360)
     478      real,intent(in) :: phy_ice(klon,360)
     479      real,intent(in) :: phy_fter(klon,360)
     480      real,intent(in) :: phy_foce(klon,360)
     481      real,intent(in) :: phy_flic(klon,360)
     482      real,intent(in) :: phy_fsic(klon,360)
     483
     484      real :: phy_glo(klon_glo,360) ! temporary variable, to store phy_***(:)
     485                                    ! on the whole physics grid
    472486      INTEGER ierr
    473487      INTEGER dimfirst(3)
    … …  
    480494      INTEGER id_FTER,id_FOCE,id_FSIC,id_FLIC
    481495 
    482       PRINT*, 'Ecriture du fichier limit'
    483 c
    484       ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid)
    485 c
    486       ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30,
     496      if (is_mpi_root.and.is_omp_root) then
     497     
     498        PRINT*, 'writelim: Ecriture du fichier limit'
     499c
     500        ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid)
     501c
     502        ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30,
    487503     .                       "Fichier conditions aux limites")
    488       ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
    489       ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim)
    490 c
    491       dims(1) = ndim
    492       dims(2) = ntim
     504!!        ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
     505        ierr = NF_DEF_DIM (nid, "points_physiques", klon_glo, ndim)
     506        ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim)
     507c
     508        dims(1) = ndim
     509        dims(2) = ntim
    493510c
    494511ccc      ierr = NF_DEF_VAR (nid, "TEMPS", NF_DOUBLE, 1,ntim, id_tim)
    495       ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim)
    496       ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17,
     512        ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim)
     513        ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17,
    497514     .                        "Jour dans l annee")
    498515ccc      ierr = NF_DEF_VAR (nid, "NAT", NF_DOUBLE, 2,dims, id_NAT)
    499       ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT)
    500       ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23,
     516        ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT)
     517        ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23,
    501518     .                        "Nature du sol (0,1,2,3)")
    502519ccc      ierr = NF_DEF_VAR (nid, "SST", NF_DOUBLE, 2,dims, id_SST)
    503       ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST)
    504       ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35,
     520        ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST)
     521        ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35,
    505522     .                        "Temperature superficielle de la mer")
    506523ccc      ierr = NF_DEF_VAR (nid, "BILS", NF_DOUBLE, 2,dims, id_BILS)
    507       ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS)
    508       ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32,
     524        ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS)
     525        ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32,
    509526     .                        "Reference flux de chaleur au sol")
    510527ccc      ierr = NF_DEF_VAR (nid, "ALB", NF_DOUBLE, 2,dims, id_ALB)
    511       ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB)
    512       ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19,
     528        ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB)
     529        ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19,
    513530     .                        "Albedo a la surface")
    514531ccc      ierr = NF_DEF_VAR (nid, "RUG", NF_DOUBLE, 2,dims, id_RUG)
    515       ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG)
    516       ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8,
     532        ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG)
     533        ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8,
    517534     .                        "Rugosite")
    518535
    519       ierr = NF_DEF_VAR (nid, "FTER", NF_FLOAT, 2,dims, id_FTER)
    520       ierr = NF_PUT_ATT_TEXT (nid, id_FTER, "title", 8,"Frac. Terre")
    521       ierr = NF_DEF_VAR (nid, "FOCE", NF_FLOAT, 2,dims, id_FOCE)
    522       ierr = NF_PUT_ATT_TEXT (nid, id_FOCE, "title", 8,"Frac. Terre")
    523       ierr = NF_DEF_VAR (nid, "FSIC", NF_FLOAT, 2,dims, id_FSIC)
    524       ierr = NF_PUT_ATT_TEXT (nid, id_FSIC, "title", 8,"Frac. Terre")
    525       ierr = NF_DEF_VAR (nid, "FLIC", NF_FLOAT, 2,dims, id_FLIC)
    526       ierr = NF_PUT_ATT_TEXT (nid, id_FLIC, "title", 8,"Frac. Terre")
    527 c
    528       ierr = NF_ENDDEF(nid)
    529 c
    530       DO k = 1, 360
    531 c
    532       debut(1) = 1
    533       debut(2) = k
    534       epais(1) = klon
    535       epais(2) = 1
    536 c
    537       print*,'Instant ',k
    538 #ifdef NC_DOUBLE
    539       print*,'NC DOUBLE'
    540       ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k))
    541       ierr = NF_PUT_VARA_DOUBLE (nid,id_NAT,debut,epais,phy_nat(1,k))
    542       ierr = NF_PUT_VARA_DOUBLE (nid,id_SST,debut,epais,phy_sst(1,k))
    543       ierr = NF_PUT_VARA_DOUBLE (nid,id_BILS,debut,epais,phy_bil(1,k))
    544       ierr = NF_PUT_VARA_DOUBLE (nid,id_ALB,debut,epais,phy_alb(1,k))
    545       ierr = NF_PUT_VARA_DOUBLE (nid,id_RUG,debut,epais,phy_rug(1,k))
    546       ierr = NF_PUT_VARA_DOUBLE (nid,id_FTER,debut,epais,phy_fter(1,k))
    547       ierr = NF_PUT_VARA_DOUBLE (nid,id_FOCE,debut,epais,phy_foce(1,k))
    548       ierr = NF_PUT_VARA_DOUBLE (nid,id_FSIC,debut,epais,phy_fsic(1,k))
    549       ierr = NF_PUT_VARA_DOUBLE (nid,id_FLIC,debut,epais,phy_flic(1,k))
    550 #else
    551       print*,'NC PAS DOUBLE'
    552       ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k))
    553       ierr = NF_PUT_VARA_REAL (nid,id_NAT,debut,epais,phy_nat(1,k))
    554       ierr = NF_PUT_VARA_REAL (nid,id_SST,debut,epais,phy_sst(1,k))
    555       ierr = NF_PUT_VARA_REAL (nid,id_BILS,debut,epais,phy_bil(1,k))
    556       ierr = NF_PUT_VARA_REAL (nid,id_ALB,debut,epais,phy_alb(1,k))
    557       ierr = NF_PUT_VARA_REAL (nid,id_RUG,debut,epais,phy_rug(1,k))
    558       ierr = NF_PUT_VARA_REAL (nid,id_FTER,debut,epais,phy_fter(1,k))
    559       ierr = NF_PUT_VARA_REAL (nid,id_FOCE,debut,epais,phy_foce(1,k))
    560       ierr = NF_PUT_VARA_REAL (nid,id_FSIC,debut,epais,phy_fsic(1,k))
    561       ierr = NF_PUT_VARA_REAL (nid,id_FLIC,debut,epais,phy_flic(1,k))
    562 
    563 #endif
    564 c
    565       ENDDO
    566 c
    567       ierr = NF_CLOSE(nid)
    568 c
    569       return
     536        ierr = NF_DEF_VAR (nid, "FTER", NF_FLOAT, 2,dims, id_FTER)
     537        ierr = NF_PUT_ATT_TEXT (nid, id_FTER, "title", 8,"Frac. Terre")
     538        ierr = NF_DEF_VAR (nid, "FOCE", NF_FLOAT, 2,dims, id_FOCE)
     539        ierr = NF_PUT_ATT_TEXT (nid, id_FOCE, "title", 8,"Frac. Terre")
     540        ierr = NF_DEF_VAR (nid, "FSIC", NF_FLOAT, 2,dims, id_FSIC)
     541        ierr = NF_PUT_ATT_TEXT (nid, id_FSIC, "title", 8,"Frac. Terre")
     542        ierr = NF_DEF_VAR (nid, "FLIC", NF_FLOAT, 2,dims, id_FLIC)
     543        ierr = NF_PUT_ATT_TEXT (nid, id_FLIC, "title", 8,"Frac. Terre")
     544c
     545        ierr = NF_ENDDEF(nid)
     546c
     547
     548! write the 'times'
     549        do k=1,360
     550#ifdef NC_DOUBLE
     551          ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k))
     552#else
     553          ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k))
     554#endif
     555        enddo
     556
     557      endif ! of if (is_mpi_root.and.is_omp_root)
     558
     559! write the fields, after having collected them on master
     560
     561      call gather(phy_nat,phy_glo)
     562      if (is_mpi_root.and.is_omp_root) then
     563#ifdef NC_DOUBLE
     564        ierr=NF_PUT_VAR_DOUBLE(nid,id_NAT,phy_glo)
     565#else
     566        ierr=NF_PUT_VAR_REAL(nid,id_NAT,phy_glo)
     567#endif
     568        if(ierr.ne.NF_NOERR) then
     569          write(*,*) "writelim error with phy_nat"
     570          write(*,*) NF_STRERROR(ierr)
     571        endif
     572      endif
     573
     574      call gather(phy_sst,phy_glo)
     575      if (is_mpi_root.and.is_omp_root) then
     576#ifdef NC_DOUBLE
     577        ierr=NF_PUT_VAR_DOUBLE(nid,id_SST,phy_glo)
     578#else
     579        ierr=NF_PUT_VAR_REAL(nid,id_SST,phy_glo)
     580#endif
     581        if(ierr.ne.NF_NOERR) then
     582          write(*,*) "writelim error with phy_sst"
     583          write(*,*) NF_STRERROR(ierr)
     584        endif
     585      endif
     586
     587      call gather(phy_bil,phy_glo)
     588      if (is_mpi_root.and.is_omp_root) then
     589#ifdef NC_DOUBLE
     590        ierr=NF_PUT_VAR_DOUBLE(nid,id_BILS,phy_glo)
     591#else
     592        ierr=NF_PUT_VAR_REAL(nid,id_BILS,phy_glo)
     593#endif
     594        if(ierr.ne.NF_NOERR) then
     595          write(*,*) "writelim error with phy_bil"
     596          write(*,*) NF_STRERROR(ierr)
     597        endif
     598      endif
     599
     600      call gather(phy_alb,phy_glo)
     601      if (is_mpi_root.and.is_omp_root) then
     602#ifdef NC_DOUBLE
     603        ierr=NF_PUT_VAR_DOUBLE(nid,id_ALB,phy_glo)
     604#else
     605        ierr=NF_PUT_VAR_REAL(nid,id_ALB,phy_glo)
     606#endif
     607        if(ierr.ne.NF_NOERR) then
     608          write(*,*) "writelim error with phy_alb"
     609          write(*,*) NF_STRERROR(ierr)
     610        endif
     611      endif
     612
     613      call gather(phy_rug,phy_glo)
     614      if (is_mpi_root.and.is_omp_root) then
     615#ifdef NC_DOUBLE
     616        ierr=NF_PUT_VAR_DOUBLE(nid,id_RUG,phy_glo)
     617#else
     618        ierr=NF_PUT_VAR_REAL(nid,id_RUG,phy_glo)
     619#endif
     620        if(ierr.ne.NF_NOERR) then
     621          write(*,*) "writelim error with phy_rug"
     622          write(*,*) NF_STRERROR(ierr)
     623        endif
     624      endif
     625
     626      call gather(phy_fter,phy_glo)
     627      if (is_mpi_root.and.is_omp_root) then
     628#ifdef NC_DOUBLE
     629        ierr=NF_PUT_VAR_DOUBLE(nid,id_FTER,phy_glo)
     630#else
     631        ierr=NF_PUT_VAR_REAL(nid,id_FTER,phy_glo)
     632#endif
     633        if(ierr.ne.NF_NOERR) then
     634          write(*,*) "writelim error with phy_fter"
     635          write(*,*) NF_STRERROR(ierr)
     636        endif
     637      endif
     638
     639      call gather(phy_foce,phy_glo)
     640      if (is_mpi_root.and.is_omp_root) then
     641#ifdef NC_DOUBLE
     642        ierr=NF_PUT_VAR_DOUBLE(nid,id_FOCE,phy_glo)
     643#else
     644        ierr=NF_PUT_VAR_REAL(nid,id_FOCE,phy_glo)
     645#endif
     646        if(ierr.ne.NF_NOERR) then
     647          write(*,*) "writelim error with phy_foce"
     648          write(*,*) NF_STRERROR(ierr)
     649        endif
     650      endif
     651
     652      call gather(phy_fsic,phy_glo)
     653      if (is_mpi_root.and.is_omp_root) then
     654#ifdef NC_DOUBLE
     655        ierr=NF_PUT_VAR_DOUBLE(nid,id_FSIC,phy_glo)
     656#else
     657        ierr=NF_PUT_VAR_REAL(nid,id_FSIC,phy_glo)
     658#endif
     659        if(ierr.ne.NF_NOERR) then
     660          write(*,*) "writelim error with phy_fsic"
     661          write(*,*) NF_STRERROR(ierr)
     662        endif
     663      endif
     664
     665      call gather(phy_flic,phy_glo)
     666      if (is_mpi_root.and.is_omp_root) then
     667#ifdef NC_DOUBLE
     668        ierr=NF_PUT_VAR_DOUBLE(nid,id_FLIC,phy_glo)
     669#else
     670        ierr=NF_PUT_VAR_REAL(nid,id_FLIC,phy_glo)
     671#endif
     672        if(ierr.ne.NF_NOERR) then
     673          write(*,*) "writelim error with phy_flic"
     674          write(*,*) NF_STRERROR(ierr)
     675        endif
     676      endif
     677
     678!  close file:
     679      if (is_mpi_root.and.is_omp_root) then
     680        ierr = NF_CLOSE(nid)
     681      endif
     682
    570683      end
     684
     685!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    571686
    572687      SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
Note: See TracChangeset for help on using the changeset viewer.