Ignore:
Timestamp:
Jun 15, 2015, 8:48:31 PM (9 years ago)
Author:
dcugnet
Message:

In dyn3d/:
etat0dyn_netcdf.F90: "startget_dyn3d" syntax slightly simplified.
dynredem.F90: Shortcut routines (put_var*, cre_var,
dynredem_write_*, dynredem_read_u)

modified to match dyn3dmem version and put in

module dyredem_mod.F90.
dynetat0.F90 -> *.f90: Few simplifications (no usage of NC_DOUBLE
needed => no precompilation)

Add tracers initialization in the isotope case

suppressed by accident.
dynredem_mod.F90: Created to mimic dyn3dmem equivalent.

In dyn3dmem/:
dynetat0_loc.F -> *.f90: Converted into fortran 90 to match the dyn3d
version.
dynredem_loc.F -> *.F90: Converted into fortran 90.
dynredem_mod.F90: Add some shortcut routines to match the dyn3d
version.

In phylmd/:
phyredem.F90: Bug fix: nsw instead of nsoilmx was used as
Tsoil second maximum index.

Bug fix: fevap instead of snow was saved for

"SNOW".
etat0phys_netcdf.F90: "filtreg_mod" module usage suppressed.

Local variable rugo computation removed (not

used).

In dynlonlat_phylonlat/:
grid_atob_m.F90 -> *.f90 DOUBLE PRECISION variables usage removed.

Precompilation o longer needed => .F90 extension.

Location:
LMDZ5/trunk/libf/dyn3d
Files:
1 added
2 edited
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/dyn3d/dynetat0.f90

    r2298 r2299  
    109109  CALL get_var2("cv"   ,cv)
    110110  CALL get_var2("aire" ,aire)
    111   CALL get_var2("phisinit",phis)
    112111  var="temps"
    113112  IF(NF90_INQ_VARID(fID,var,vID)/=NF90_NoErr) THEN
    … …  
    117116  END IF
    118117  CALL err(NF90_GET_VAR(fID,vID,time),"get",var)
     118  CALL get_var2("phisinit",phis)
    119119  CALL get_var3("ucov",ucov)
    120120  CALL get_var3("vcov",vcov)
    … …  
    126126  DO iq=1,nqtot
    127127    var=tname(iq)
    128     IF(NF90_INQ_VARID(fID,var,vID)/=NF90_NoErr) THEN
    129       WRITE(lunout,*)TRIM(modname)//": Tracer <"//TRIM(var)//"> is missing"
    130       WRITE(lunout,*)"         It is hence initialized to zero"
    131       q(:,:,:,iq)=0.
    132     ELSE
    133       CALL err(NF90_GET_VAR(fID,vID,q(:,:,:,iq)),"get",var)
     128    IF(NF90_INQ_VARID(fID,var,vID)==NF90_NoErr) THEN
     129      CALL err(NF90_GET_VAR(fID,vID,q(:,:,:,iq)),"get",var); CYCLE
     130    END IF
     131    WRITE(lunout,*)TRIM(modname)//": Tracer <"//TRIM(var)//"> is missing"
     132    WRITE(lunout,*)"         It is hence initialized to zero"
     133    q(:,:,:,iq)=0.
     134   !--- CRisi: for isotops, theoretical initialization using very simplified
     135   !           Rayleigh distillation las.
     136    IF(ok_isotopes.AND.iso_num(iq)>0) THEN
     137      IF(zone_num(iq)==0) q(:,:,:,iq)=q(:,:,:,iqpere(iq))*tnat(iso_num(iq))    &
     138     &             *(q(:,:,:,iqpere(iq))/30.e-3)**(alpha_ideal(iso_num(iq))-1)
     139      IF(zone_num(iq)==1) q(:,:,:,iq)=q(:,:,:,iqiso(iso_indnum(iq),phase_num(iq)))
    134140    END IF
    135141  END DO
    … …  
    158164SUBROUTINE get_var1(var,v)
    159165  CHARACTER(LEN=*), INTENT(IN)  :: var
    160 #ifdef NC_DOUBLE
    161   DOUBLE PRECISION, INTENT(OUT) :: v(:)
    162 #else
    163166  REAL,             INTENT(OUT) :: v(:)
    164 #endif
    165167  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
    166168  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
    … …  
    170172SUBROUTINE get_var2(var,v)
    171173  CHARACTER(LEN=*), INTENT(IN)  :: var
    172 #ifdef NC_DOUBLE
    173   DOUBLE PRECISION, INTENT(OUT) :: v(:,:)
    174 #else
    175174  REAL,             INTENT(OUT) :: v(:,:)
    176 #endif
    177175  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
    178176  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
    … …  
    182180SUBROUTINE get_var3(var,v)
    183181  CHARACTER(LEN=*), INTENT(IN)  :: var
    184 #ifdef NC_DOUBLE
    185   DOUBLE PRECISION, INTENT(OUT) :: v(:,:,:)
    186 #else
    187182  REAL,             INTENT(OUT) :: v(:,:,:)
    188 #endif
    189183  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
    190184  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)

  • LMDZ5/trunk/libf/dyn3d/dynredem.F90

    r2293 r2299  
    88#endif
    99  USE infotrac
    10   USE netcdf, ONLY:   NF90_CREATE, NF90_DEF_DIM, NF90_REDEF,  NF90_INQ_VARID, &
    11       NF90_CLOBBER,   NF90_CLOSE,  NF90_DEF_VAR, NF90_ENDDEF, NF90_PUT_ATT,   &
    12       NF90_UNLIMITED, NF90_GLOBAL, NF90_FLOAT,   NF90_DOUBLE
    13   USE netcdf95, ONLY: NF95_PUT_VAR
     10  USE netcdf, ONLY: NF90_CREATE, NF90_DEF_DIM, NF90_INQ_VARID, NF90_GLOBAL,    &
     11                    NF90_CLOSE,  NF90_PUT_ATT, NF90_UNLIMITED, NF90_CLOBBER
     12  USE dynredem_mod, ONLY: cre_var, put_var1, put_var2, err, modname, fil
    1413  IMPLICIT NONE
    1514  include "dimensions.h"
    … …  
    2120  include "ener.h"
    2221  include "logic.h"
    23   include "netcdf.inc"
    2422  include "description.h"
    2523  include "serre.h"
    … …  
    3533  INTEGER, PARAMETER :: length=100
    3634  REAL    :: tab_cntrl(length)                     !--- RUN PARAMETERS TABLE
    37   CHARACTER(LEN=20) :: modname
    3835!   For NetCDF:
    3936  CHARACTER(LEN=30) :: unites
    … …  
    4239  INTEGER :: sID, sigID, nID, vID, timID
    4340  INTEGER :: yyears0, jjour0, mmois0
    44   REAL :: zan0, zjulian, hours
    45 !===============================================================================
    46   modname='dynredem0'
     41  REAL    :: zan0, zjulian, hours
     42!===============================================================================
     43  modname='dynredem0'; fil=fichnom
    4744#ifdef CPP_IOIPSL
    4845  CALL ymds2ju(annee_ref, 1, iday_end, 0.0, zjulian)
    … …  
    10299! start_time: start_time of simulation (not necessarily 0.)
    103100  tab_cntrl(32) = start_time
    104 !.........................................................
    105101
    106102!--- File creation
    … …  
    121117
    122118!--- Define and save invariant fields
    123   CALL put_var1("controle","Parametres de controle" ,[indexID],tab_cntrl)
    124   CALL put_var1("rlonu"   ,"Longitudes des points U",[rlonuID],rlonu)
    125   CALL put_var1("rlatu"   ,"Latitudes des points U" ,[rlatuID],rlatu)
    126   CALL put_var1("rlonv"   ,"Longitudes des points V",[rlonvID],rlonv)
    127   CALL put_var1("rlatv"   ,"Latitudes des points V" ,[rlatvID],rlatv)
    128   CALL put_var1("nivsigs" ,"Numero naturel des couches s"    ,[sID]  ,nivsigs)
    129   CALL put_var1("nivsig"  ,"Numero naturel des couches sigma",[sigID],nivsig)
    130   CALL put_var1("ap"      ,"Coefficient A pour hybride"      ,[sigID],ap)
    131   CALL put_var1("bp"      ,"Coefficient B pour hybride"      ,[sigID],bp)
    132   CALL put_var1("presnivs",""                                ,[sID]  ,presnivs)
     119  CALL put_var1(nid,"controle","Parametres de controle" ,[indexID],tab_cntrl)
     120  CALL put_var1(nid,"rlonu"   ,"Longitudes des points U",[rlonuID],rlonu)
     121  CALL put_var1(nid,"rlatu"   ,"Latitudes des points U" ,[rlatuID],rlatu)
     122  CALL put_var1(nid,"rlonv"   ,"Longitudes des points V",[rlonvID],rlonv)
     123  CALL put_var1(nid,"rlatv"   ,"Latitudes des points V" ,[rlatvID],rlatv)
     124  CALL put_var1(nid,"nivsigs" ,"Numero naturel des couches s"    ,[sID]  ,nivsigs)
     125  CALL put_var1(nid,"nivsig"  ,"Numero naturel des couches sigma",[sigID],nivsig)
     126  CALL put_var1(nid,"ap"      ,"Coefficient A pour hybride"      ,[sigID],ap)
     127  CALL put_var1(nid,"bp"      ,"Coefficient B pour hybride"      ,[sigID],bp)
     128  CALL put_var1(nid,"presnivs",""                                ,[sID]  ,presnivs)
    133129! covariant <-> contravariant <-> natural conversion coefficients
    134   CALL put_var2("cu","Coefficient de passage pour U",[rlonuID,rlatuID],cu)
    135   CALL put_var2("cv","Coefficient de passage pour V",[rlonvID,rlatvID],cv)
    136   CALL put_var2("aire","Aires de chaque maille"     ,[rlonvID,rlatuID],aire)
    137   CALL put_var2("phisinit","Geopotentiel au sol"    ,[rlonvID,rlatuID],phis)
     130  CALL put_var2(nid,"cu","Coefficient de passage pour U",[rlonuID,rlatuID],cu)
     131  CALL put_var2(nid,"cv","Coefficient de passage pour V",[rlonvID,rlatvID],cv)
     132  CALL put_var2(nid,"aire","Aires de chaque maille"     ,[rlonvID,rlatuID],aire)
     133  CALL put_var2(nid,"phisinit","Geopotentiel au sol"    ,[rlonvID,rlatuID],phis)
    138134
    139135!--- Define fields saved later
    140136  WRITE(unites,"('days since ',i4,'-',i2.2,'-',i2.2,' 00:00:00')"),&
    141137               yyears0,mmois0,jjour0
    142   CALL put_var0("temps","Temps de simulation",[timID],unites)
    143   CALL put_var0("ucov" ,"Vitesse U"  ,[rlonuID,rlatuID,sID,timID])
    144   CALL put_var0("vcov" ,"Vitesse V"  ,[rlonvID,rlatvID,sID,timID])
    145   CALL put_var0("teta" ,"Temperature",[rlonvID,rlatuID,sID,timID])
     138  CALL cre_var(nid,"temps","Temps de simulation",[timID],unites)
     139  CALL cre_var(nid,"ucov" ,"Vitesse U"  ,[rlonuID,rlatuID,sID,timID])
     140  CALL cre_var(nid,"vcov" ,"Vitesse V"  ,[rlonvID,rlatvID,sID,timID])
     141  CALL cre_var(nid,"teta" ,"Temperature",[rlonvID,rlatuID,sID,timID])
    146142  DO iq=1,nqtot
    147     CALL put_var0(tname(iq),ttext(iq) ,[rlonvID,rlatuID,sID,timID])
     143    CALL cre_var(nid,tname(iq),ttext(iq),[rlonvID,rlatuID,sID,timID])
    148144  END DO
    149   CALL put_var0("masse","Masse d air"    ,[rlonvID,rlatuID,sID,timID])
    150   CALL put_var0("ps"   ,"Pression au sol",[rlonvID,rlatuID    ,timID])
    151   CALL err(NF90_ENDDEF(nid))
     145  CALL cre_var(nid,"masse","Masse d air"    ,[rlonvID,rlatuID,sID,timID])
     146  CALL cre_var(nid,"ps"   ,"Pression au sol",[rlonvID,rlatuID    ,timID])
    152147  CALL err(NF90_CLOSE (nid))
    153148
    … …  
    155150  WRITE(lunout,*)TRIM(modname)//': rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
    156151
    157 
    158 CONTAINS
    159 
    160 
    161 SUBROUTINE put_var0(var,title,did,units)
    162   CHARACTER(LEN=*),           INTENT(IN) :: var, title
    163   INTEGER,                    INTENT(IN) :: did(:)
    164   CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: units
    165 #ifdef NC_DOUBLE
    166   CALL err(NF90_DEF_VAR(nid,var,NF90_DOUBLE,did,vID),var)
    167 #else
    168   CALL err(NF90_DEF_VAR(nid,var,NF90_FLOAT,did,vID),var)
    169 #endif
    170   IF(title/="")      CALL err(NF90_PUT_ATT(nid,vID,"title",title),var)
    171   IF(PRESENT(units)) CALL err(NF90_PUT_ATT(nid,vID,"units",units),var)
    172 END SUBROUTINE put_var0
    173 
    174 
    175 SUBROUTINE put_var1(var,title,did,v,units)
    176   CHARACTER(LEN=*),           INTENT(IN) :: var, title
    177   INTEGER,                    INTENT(IN) :: did(1)
    178 #ifdef NC_DOUBLE
    179   DOUBLE PRECISION,           INTENT(IN) :: v(:)
    180 #else
    181   REAL,                       INTENT(IN) :: v(:)
    182 #endif
    183   CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: units
    184 #ifdef NC_DOUBLE
    185   CALL err(NF90_DEF_VAR(nid,var,NF90_DOUBLE,did,vID),var)
    186 #else
    187   CALL err(NF90_DEF_VAR(nid,var,NF90_FLOAT,did,vID),var)
    188 #endif
    189   IF(title/="")      CALL err(NF90_PUT_ATT(nid,vID,"title",title),var)
    190   IF(PRESENT(units)) CALL err(NF90_PUT_ATT(nid,vID,"units",units),var)
    191   CALL err(NF90_ENDDEF(nid))
    192   CALL NF95_PUT_VAR(nid,vID,v)
    193   CALL err(NF90_REDEF(nid))
    194 END SUBROUTINE put_var1
    195 
    196 
    197 SUBROUTINE put_var2(var,title,did,v,units)
    198   CHARACTER(LEN=*),           INTENT(IN) :: var, title
    199   INTEGER,                    INTENT(IN) :: did(2)
    200 #ifdef NC_DOUBLE
    201   DOUBLE PRECISION,           INTENT(IN) :: v(:,:)
    202 #else
    203   REAL,                       INTENT(IN) :: v(:,:)
    204 #endif
    205   CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: units
    206 #ifdef NC_DOUBLE
    207   CALL err(NF90_DEF_VAR(nid,var,NF90_DOUBLE,did,vID),var)
    208 #else
    209   CALL err(NF90_DEF_VAR(nid,var,NF90_FLOAT,did,vID),var)
    210 #endif
    211   IF(title/="")      CALL err(NF90_PUT_ATT(nid,vID,"title",title),var)
    212   IF(PRESENT(units)) CALL err(NF90_PUT_ATT(nid,vID,"units",units),var)
    213   CALL err(NF90_ENDDEF(nid))
    214   CALL NF95_PUT_VAR(nid,vID,v)
    215   CALL err(NF90_REDEF(nid))
    216 
    217 END SUBROUTINE put_var2
    218 
    219 
    220 SUBROUTINE err(ierr,var)
    221   USE netcdf, ONLY: NF90_STRERROR, NF90_NOERR
    222   IMPLICIT NONE
    223   INTEGER,                    INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
    224   CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: var    !--- VARIABLE NAME
    225   CHARACTER(LEN=256) :: file, msg
    226   IF(ierr==NF90_NoERR) RETURN
    227   msg='Error in "'//TRIM(modname)//'" for file "'//TRIM(fichnom)//'"'
    228   IF(PRESENT(var)) msg=TRIM(msg)//'" and variable "'//TRIM(var)//'"'
    229   WRITE(lunout,*)TRIM(msg)//': '//NF90_STRERROR(ierr)
    230 
    231 END SUBROUTINE err
    232 
    233152END SUBROUTINE dynredem0
    234153!
    … …  
    245164  USE infotrac
    246165  USE control_mod
    247   USE netcdf,   ONLY: NF90_OPEN,  NF90_NOWRITE, NF90_INQ_VARID, NF90_NoErr,    &
    248                       NF90_CLOSE, NF90_WRITE,   NF90_GET_VAR
    249   USE netcdf95, ONLY: NF95_PUT_VAR
    250   USE assert_eq_m, ONLY: assert_eq
     166  USE netcdf,   ONLY: NF90_OPEN,  NF90_NOWRITE, NF90_GET_VAR, NF90_INQ_VARID,  &
     167                      NF90_CLOSE, NF90_WRITE,   NF90_PUT_VAR, NF90_NoErr
     168  USE dynredem_mod, ONLY: dynredem_write_u, dynredem_write_v, dynredem_read_u, &
     169                          err, modname, fil, msg
    251170  IMPLICIT NONE
    252 #include "dimensions.h"
    253 #include "paramet.h"
    254 #include "description.h"
    255 #include "comvert.h"
    256 #include "comgeom.h"
    257 #include "temps.h"
    258 #include "iniprint.h"
     171  include "dimensions.h"
     172  include "paramet.h"
     173  include "description.h"
     174  include "comvert.h"
     175  include "comgeom.h"
     176  include "temps.h"
     177  include "iniprint.h"
    259178!===============================================================================
    260179! Arguments:
    261   CHARACTER(LEN=*), INTENT(IN) :: fichnom           !-- FILE NAME
    262   REAL, INTENT(IN) ::  time                         !-- TIME
    263   REAL, INTENT(IN) ::  vcov(iip1,jjm, llm)          !-- V COVARIANT WIND
    264   REAL, INTENT(IN) ::  ucov(iip1,jjp1,llm)          !-- U COVARIANT WIND
    265   REAL, INTENT(IN) ::  teta(iip1,jjp1,llm)          !-- POTENTIAL TEMPERATURE
    266   REAL, INTENT(IN) ::     q(iip1,jjp1,llm,nqtot)    !-- TRACERS
    267   REAL, INTENT(IN) :: masse(iip1,jjp1,llm)          !-- MASS PER CELL
    268   REAL, INTENT(IN) ::    ps(iip1,jjp1)              !-- GROUND PRESSURE
     180  CHARACTER(LEN=*), INTENT(IN) :: fichnom              !-- FILE NAME
     181  REAL, INTENT(IN)    ::  time                         !-- TIME
     182  REAL, INTENT(IN)    ::  vcov(iip1,jjm, llm)          !-- V COVARIANT WIND
     183  REAL, INTENT(IN)    ::  ucov(iip1,jjp1,llm)          !-- U COVARIANT WIND
     184  REAL, INTENT(IN)    ::  teta(iip1,jjp1,llm)          !-- POTENTIAL TEMPERATURE
     185  REAL, INTENT(INOUT) ::     q(iip1,jjp1,llm,nqtot)    !-- TRACERS
     186  REAL, INTENT(IN)    :: masse(iip1,jjp1,llm)          !-- MASS PER CELL
     187  REAL, INTENT(IN)    ::    ps(iip1,jjp1)              !-- GROUND PRESSURE
    269188!===============================================================================
    270189! Local variables:
    271   INTEGER :: l, iq, nid, vID, nid_trac, vID_trac
     190  INTEGER :: l, iq, nid, vID, ierr, nid_trac, vID_trac
    272191  INTEGER, SAVE :: nb=0
    273192  INTEGER, PARAMETER :: length=100
    274 #ifdef NC_DOUBLE
    275   DOUBLE PRECISION   :: trac_tmp(ip1jmp1,llm)
    276 #else
    277   REAL               :: trac_tmp(ip1jmp1,llm)
    278 #endif
    279193  REAL               :: tab_cntrl(length) ! tableau des parametres du run
    280   CHARACTER(LEN=256) :: modname, var, fil
    281   LOGICAL            :: exist_file
    282 !===============================================================================
    283   modname='dynredem1'
    284   fil=fichnom
     194  CHARACTER(LEN=256) :: var, dum
     195  LOGICAL            :: lread_inca
     196!===============================================================================
     197
     198  modname='dynredem1'; fil=fichnom
    285199  CALL err(NF90_OPEN(fil,NF90_WRITE,nid),"open",fil)
    286200
    287201!--- Write/extend time coordinate
    288202  nb = nb + 1
    289   CALL sav_var1("temps",[time],nb)
     203  var="temps"
     204  CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var)
     205  CALL err(NF90_PUT_VAR(nid,vID,[time]),"put",var)
    290206  WRITE(lunout,*)TRIM(modname)//": Saving for ", nb, time
    291207
    292208!--- Rewrite control table (itaufin undefined in dynredem0)
    293209  var="controle"
    294   CALL get_var1(var,tab_cntrl); tab_cntrl(31)=DBLE(itau_dyn + itaufin)
    295210  CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var)
    296   CALL NF95_PUT_VAR(nid,vID,tab_cntrl)
     211  CALL err(NF90_GET_VAR(nid,vID,tab_cntrl),"get",var)
     212  tab_cntrl(31)=DBLE(itau_dyn + itaufin)
     213  CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var)
     214  CALL err(NF90_PUT_VAR(nid,vID,tab_cntrl),"put",var)
    297215
    298216!--- Save fields
    299   CALL sav_var3("ucov",ucov)
    300   CALL sav_var3("vcov",vcov)
    301   CALL sav_var3("teta",teta)
    302   CALL sav_var3("masse",masse)
    303   CALL sav_var2("ps"   ,ps)
     217  CALL dynredem_write_u(nid,"ucov" ,ucov ,llm)
     218  CALL dynredem_write_v(nid,"vcov" ,vcov ,llm)
     219  CALL dynredem_write_u(nid,"teta" ,teta ,llm)
     220  CALL dynredem_write_u(nid,"masse",masse,llm)
     221  CALL dynredem_write_u(nid,"ps"   ,ps   ,1)
    304222
    305223!--- Tracers in file "start_trac.nc" (added by Anne)
    306   IF (type_trac == 'inca') THEN
    307     fil="start_trac.nc"; INQUIRE(FILE=fil,EXIST=exist_file)
    308     IF(.NOT.exist_file) CALL war(-1,"open",fil)
    309   END IF
    310   DO iq=1,nqtot; var=tname(iq)
    311 
    312   !--- Usual case
    313     IF(type_trac/='inca') THEN
    314       CALL sav_var3(var,q(:,:,:,iq)); CYCLE
     224  lread_inca=.FALSE.; fil="start_trac.nc"
     225  IF(type_trac=='inca') INQUIRE(FILE=fil,EXIST=lread_inca)
     226  IF(lread_inca) CALL err(NF90_OPEN(fil,NF90_NOWRITE,nid_trac),"open")
     227
     228!--- Save tracers
     229  DO iq=1,nqtot; var=tname(iq); ierr=-1
     230    IF(lread_inca) THEN                  !--- Possibly read from "start_trac.nc"
     231      fil="start_trac.nc"
     232      ierr=NF90_INQ_VARID(nid_trac,var,vID_trac)
     233      dum='inq'; IF(ierr==NF90_NoErr) dum='fnd'
     234      WRITE(lunout,*)msg(dum,var)
     235
     236
     237      IF(ierr==NF90_NoErr) CALL dynredem_read_u(nid_trac,var,q(:,:,:,iq),llm)
    315238    END IF
    316 
    317   !--- Special case for INCA tracer read from "start_trac.nc"
    318     IF(NF90_INQ_VARID(nid_trac,var,vID_trac)/=NF90_NoErr) THEN
    319       CALL war(-1,"inq",var,fil)
    320       CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var,fil)
    321       CALL NF95_PUT_VAR(nid,vID,q(:,:,:,iq))
    322     ELSE
    323       WRITE(lunout,*)TRIM(modname)//": <"//TRIM(var)//"> found in "//fil
    324       CALL err(NF90_GET_VAR(nid_trac,vID_trac,trac_tmp),"get",var,fil)
    325     END IF
    326     CALL sav_var3(var,RESHAPE(trac_tmp,SHAPE=[iip1,jjp1,llm]))
     239    fil=fichnom
     240    CALL dynredem_write_u(nid,var,q(:,:,:,iq),llm)
    327241  END DO
    328   CALL err(NF90_CLOSE(nid),"close",fichnom)
    329 
    330 
    331 CONTAINS
    332 
    333 
    334 SUBROUTINE get_var1(var,v)
    335   CHARACTER(LEN=*), INTENT(IN) :: var
    336   REAL,             INTENT(OUT) :: v(:)
    337   CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var)
    338   CALL err(NF90_GET_VAR(nid,vID,v),"get",var)
    339 END SUBROUTINE get_var1
    340 
    341 
    342 SUBROUTINE sav_var1(var,v,start)
    343   CHARACTER(LEN=*),  INTENT(IN) :: var
    344 #ifdef NC_DOUBLE
    345   DOUBLE PRECISION,  INTENT(IN) :: v(:)
    346 #else
    347   REAL,              INTENT(IN) :: v(:)
    348 #endif
    349   INTEGER, OPTIONAL, INTENT(IN) :: start
    350   CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var)
    351   IF(PRESENT(start)) THEN
    352     CALL NF95_PUT_VAR(nid,vID,v,start=[start])
    353   ELSE
    354     CALL NF95_PUT_VAR(nid,vID,v)
    355   END IF
    356 END SUBROUTINE sav_var1
    357 
    358 
    359 SUBROUTINE sav_var2(var,v)
    360   CHARACTER(LEN=*), INTENT(IN) :: var
    361 #ifdef NC_DOUBLE
    362   DOUBLE PRECISION, INTENT(IN) :: v(:,:)
    363 #else
    364   REAL,             INTENT(IN) :: v(:,:)
    365 #endif
    366   CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var)
    367   CALL NF95_PUT_VAR(nid,vID,v)
    368 END SUBROUTINE sav_var2
    369 
    370 
    371 SUBROUTINE sav_var3(var,v)
    372   CHARACTER(LEN=*), INTENT(IN) :: var
    373 #ifdef NC_DOUBLE
    374   DOUBLE PRECISION, INTENT(IN) :: v(:,:,:)
    375 #else
    376   REAL,             INTENT(IN) :: v(:,:,:)
    377 #endif
    378 
    379 print*,'var='//TRIM(var)
    380 print*,SHAPE(v)
    381   CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var)
    382   CALL NF95_PUT_VAR(nid,vID,v)
    383 END SUBROUTINE sav_var3
    384 
    385 
    386 FUNCTION msg(typ,nam,fil)
    387   IMPLICIT NONE
    388   CHARACTER(LEN=256)           :: msg    !--- STANDARDIZED MESSAGE
    389   CHARACTER(LEN=*), INTENT(IN) :: typ    !--- TYPE OF OPERATION
    390   CHARACTER(LEN=*), INTENT(IN) :: nam    !--- FIELD NAME
    391   CHARACTER(LEN=*), INTENT(IN) :: fil    !--- FILE  NAME
    392   SELECT CASE(typ)
    393     CASE('inq');   msg="Missing field <"//TRIM(nam)//">"
    394     CASE('get');   msg="Reading failed for <"//TRIM(nam)//">"
    395     CASE('open');  msg="Opening failed for <"//TRIM(nam)//">"
    396     CASE('close'); msg="Closing failed for <"//TRIM(nam)//">"
    397   END SELECT
    398   msg=TRIM(modname)//": "//TRIM(msg)
    399   IF(typ=="inq".AND.fil/="") msg=TRIM(msg)//" in file <"//TRIM(fil)//">"
    400 
    401 END FUNCTION msg
    402 
    403 
    404 SUBROUTINE err(ierr,typ,nam,fil)
    405   IMPLICIT NONE
    406   INTEGER,                    INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
    407   CHARACTER(LEN=*),           INTENT(IN) :: typ    !--- TYPE OF OPERATION
    408   CHARACTER(LEN=*),           INTENT(IN) :: nam    !--- FIELD NAME
    409   CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: fil    !--- FILE  NAME
    410   CHARACTER(LEN=256) :: file
    411   IF(ierr==NF90_NoERR) RETURN
    412   file=""; IF(PRESENT(fil)) file=fil
    413   CALL ABORT_gcm(modname,msg(typ,nam,file),ierr)
    414 END SUBROUTINE err
    415 
    416 
    417 SUBROUTINE war(ierr,typ,nam,fil)
    418   IMPLICIT NONE
    419   INTEGER,                    INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
    420   CHARACTER(LEN=*),           INTENT(IN) :: typ    !--- TYPE OF OPERATION
    421   CHARACTER(LEN=*),           INTENT(IN) :: nam    !--- FIELD NAME
    422   CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: fil    !--- FILE  NAME
    423   CHARACTER(LEN=256) :: file
    424   IF(ierr==NF90_NoERR) RETURN
    425   file=""; IF(PRESENT(fil)) file=fil
    426   WRITE(lunout,*)msg(typ,nam,file)
    427 END SUBROUTINE war
    428 
     242  CALL err(NF90_CLOSE(nid),"close")
     243  fil="start_trac.nc"
     244  IF(lread_inca) CALL err(NF90_CLOSE(nid_trac),"close")
    429245
    430246END SUBROUTINE dynredem1

  • LMDZ5/trunk/libf/dyn3d/etat0dyn_netcdf.F90

    r2293 r2299  
    1212!  routine (to be called after restget):
    1313!    CALL startget_dyn3d(varname, lon_in,  lat_in, pls, workvar,&
    14 !                 champ, val_exp, lon_in2, lat_in2, ibar)
     14!                          champ, lon_in2, lat_in2, ibar)
    1515!
    1616!    *  Variables should have the following names in the NetCDF files:
    … …  
    8787  USE infotrac
    8888  USE filtreg_mod
    89 !#endif
    9089  IMPLICIT NONE
    9190!-------------------------------------------------------------------------------
    … …  
    120119!*******************************************************************************
    121120  CALL infotrac_init
    122   ALLOCATE(q3d(iip1,jjp1,llm,nqtot))
    123121  CALL inifilr()
    124122
    … …  
    154152! Update uvent, vvent, t3d and tpot
    155153!*******************************************************************************
    156   uvent(:,:,:) = 0.0 ; t3d (:,:,:) = 0.0
    157   vvent(:,:,:) = 0.0 ; tpot(:,:,:) = 0.0
    158   CALL startget_dyn3d('u'   ,rlonu,rlatu,pls,y ,uvent,0.0,rlonv,rlatv,ib)
    159   CALL startget_dyn3d('v'   ,rlonv,rlatv,pls(:,:jjm,:),y(:,:jjm,:),vvent,0.0,  &
     154  uvent(:,:,:) = 0.0 ; vvent(:,:,:) = 0.0 ; t3d (:,:,:) = 0.0
     155  CALL startget_dyn3d('u'   ,rlonu,rlatu,pls,y ,uvent,rlonv,rlatv,ib)
     156  CALL startget_dyn3d('v'   ,rlonv,rlatv,pls(:,:jjm,:),y(:,:jjm,:),vvent,      &
    160157 &                           rlonu,rlatu(:jjm),ib)
    161   CALL startget_dyn3d('t'   ,rlonv,rlatu,pls,y ,t3d ,0.0,rlonu,rlatv,ib)
    162   tpot=t3d
    163   CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,0.0,rlonu,rlatv,ib)
     158  CALL startget_dyn3d('t'   ,rlonv,rlatu,pls,y ,t3d ,rlonu,rlatv,ib)
     159  tpot(:,:,:)=t3d(:,:,:)
     160  CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,rlonu,rlatv,ib)
    164161
    165162  WRITE(lunout,*) 'T3D min,max:',MINVAL(t3d(:,:,:)),MAXVAL(t3d(:,:,:))
    … …  
    174171!  WRITE(lunout,*) 'QSAT :',qsat(10,20,:)
    175172  qd (:,:,:) = 0.0
    176   CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,0.0,rlonu,rlatv,ib)
    177   q3d(:,:,:,:) = 0.0 ; q3d(:,:,:,1) = qd(:,:,:)
    178 
     173  CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,rlonu,rlatv,ib)
     174  ALLOCATE(q3d(iip1,jjp1,llm,nqtot)); q3d(:,:,:,:)=0.0 ; q3d(:,:,:,1)=qd(:,:,:)
    179175  CALL flinclo(fid_dyn)
    180176
    181177#ifdef CPP_PHYS
     178#ifdef CPP_EARTH
    182179! Parameterization of ozone chemistry:
    183180!*******************************************************************************
    … …  
    190187    q3d(:,:,:,i)=q3d(:,:,:,i)*48./ 29.                  !--- Mole->mass fraction         
    191188  END IF
     189
    192190#endif
     191#endif
     192  q3d(iip1,:,:,:)=q3d(1,:,:,:)
    193193
    194194! Intermediate computation
    … …  
    204204    masse(:,jjp1,l)=xps
    205205  END DO
    206   q3d(iip1,:,:,:)=q3d(1,:,:,:)
    207206
    208207! Writing
    … …  
    234233
    235234!#endif
    236 !#endif of #ifdef CPP_EARTH
     235! of ifdef CPP_EARTH
    237236
    238237END SUBROUTINE etat0dyn_netcdf
    … …  
    244243!-------------------------------------------------------------------------------
    245244!
    246 SUBROUTINE startget_dyn3d(var,  lon_in,  lat_in,  pls,  workvar,&
    247                 champ, val_exp, lon_in2, lat_in2, ibar)
     245SUBROUTINE startget_dyn3d(var, lon_in,  lat_in,  pls,  workvar,&
     246                        champ, lon_in2, lat_in2, ibar)
    248247!-------------------------------------------------------------------------------
    249248  IMPLICIT NONE
    … …  
    253252!-------------------------------------------------------------------------------
    254253! Note: An input auxilliary field "workvar" has to be specified in two cases:
    255 !     * for "q":     the saturated humidity.
    256 !     * for "topot": the Exner function.
     254!     * for "q":    the saturated humidity.
     255!     * for "tpot": the Exner function.
    257256!===============================================================================
    258257! Arguments:
    … …  
    263262  REAL,             INTENT(IN)    :: workvar(:, :, :) ! dim (iml, jml, lml)
    264263  REAL,             INTENT(INOUT) :: champ  (:, :, :) ! dim (iml, jml, lml)
    265   REAL,             INTENT(IN)    :: val_exp
    266264  REAL,             INTENT(IN)    :: lon_in2(:)       ! dim (iml)
    267265  REAL,             INTENT(IN)    :: lat_in2(:)       ! dim (jml2)
    … …  
    274272  REAL               :: xppn, xpps
    275273!-------------------------------------------------------------------------------
    276   IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN
    277     iml = assert_eq([SIZE(lon_in),SIZE(pls,1),SIZE(workvar,1),SIZE(champ,1),  &
    278      &               SIZE(lon_in2)],TRIM(modname)//" iml")
    279     jml = assert_eq( SIZE(lat_in),SIZE(pls,2),SIZE(workvar,2),SIZE(champ,2),  &
    280      &                              TRIM(modname)//" jml")
    281     lml = assert_eq(              SIZE(pls,3),SIZE(workvar,3),SIZE(champ,3),  &
    282      &                              TRIM(modname)//" lml")
    283     jml2 = SIZE(lat_in2)
    284 
    285   !--- CHECK IF THE FIELD IS KNOWN
    286     SELECT CASE(var)
    287       CASE('u');    vname='U'
    288       CASE('v');    vname='V'
    289       CASE('t');    vname='TEMP'
    290       CASE('q');    vname='R';    msg='humidity as the saturated humidity'
    291       CASE('tpot'); vname='TEMP'; msg='potential temperature as the Exner function'
    292       CASE DEFAULT;               msg='No rule to extract variable '//TRIM(var)
    293         CALL abort_gcm(modname,TRIM(msg)//' from any data set',1)
    294     END SELECT
    295 
    296   !--- CHECK IF SOMETHING IS MISSING
    297     IF((var=='tpot'.OR.var=='q').AND.MINVAL(workvar)==MAXVAL(workvar)) THEN
    298       msg='Could not compute '//TRIM(msg)//' is missing or constant.'
    299       CALL abort_gcm(modname,TRIM(msg),1)
    300     END IF
    301 
    302   !--- INTERPOLATE 3D FIELD IF NEEDED
    303     IF(var/='tpot') CALL start_inter_3d(TRIM(vname),lon_in,lat_in,lon_in2,    &
    304                                                     lat_in2,pls,champ,ibar)
    305 
    306   !--- COMPUTE THE REQUIRED FILED
    307     SELECT CASE(var)
    308       CASE('u'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO
    309         champ(iml,:,:)=champ(1,:,:)                        !--- Eastward wind
    310 
    311       CASE('v'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO
    312         champ(iml,:,:)=champ(1,:,:)                        !--- Northward wind
    313 
    314       CASE('tpot','q')
    315         IF(var=='tpot') THEN; champ=champ*cpp/workvar      !--- Temperature
    316         ELSE;                 champ=champ*.01*workvar      !--- Relat. humidity
    317           WHERE(champ<0.) champ=1.0E-10
    318         END IF
    319         DO il=1,lml
    320           xppn = SUM(aire(:,1  )*champ(:,1  ,il))/apoln
    321           xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
    322           champ(:,1  ,il) = xppn
    323           champ(:,jml,il) = xpps
    324         END DO
    325     END SELECT
     274  iml=assert_eq([SIZE(lon_in),SIZE(pls,1),SIZE(workvar,1),SIZE(champ,1),       &
     275     &                                    SIZE(lon_in2)], TRIM(modname)//" iml")
     276  jml=assert_eq( SIZE(lat_in),SIZE(pls,2),SIZE(workvar,2),SIZE(champ,2),       &
     277     &                                                    TRIM(modname)//" jml")
     278  lml=assert_eq(              SIZE(pls,3),SIZE(workvar,3),SIZE(champ,3),       &
     279     &                                                    TRIM(modname)//" lml")
     280  jml2=SIZE(lat_in2)
     281
     282!--- CHECK IF THE FIELD IS KNOWN
     283   SELECT CASE(var)
     284    CASE('u');    vname='U'
     285    CASE('v');    vname='V'
     286    CASE('t');    vname='TEMP'
     287    CASE('q');    vname='R';    msg='humidity as the saturated humidity'
     288    CASE('tpot'); msg='potential temperature as the Exner function'
     289    CASE DEFAULT; msg='No rule to extract variable '//TRIM(var)
     290      CALL abort_gcm(modname,TRIM(msg)//' from any data set',1)
     291  END SELECT
     292
     293!--- CHECK IF SOMETHING IS MISSING
     294  IF((var=='tpot'.OR.var=='q').AND.MINVAL(workvar)==MAXVAL(workvar)) THEN
     295    msg='Could not compute '//TRIM(msg)//' is missing or constant.'
     296    CALL abort_gcm(modname,TRIM(msg),1)
    326297  END IF
     298
     299!--- INTERPOLATE 3D FIELD IF NEEDED
     300  IF(var/='tpot') CALL start_inter_3d(TRIM(vname),lon_in,lat_in,lon_in2,      &
     301                                                  lat_in2,pls,champ,ibar)
     302
     303!--- COMPUTE THE REQUIRED FILED
     304  SELECT CASE(var)
     305    CASE('u'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO
     306      champ(iml,:,:)=champ(1,:,:)                   !--- Eastward wind
     307
     308    CASE('v'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO
     309      champ(iml,:,:)=champ(1,:,:)                   !--- Northward wind
     310
     311    CASE('tpot','q')
     312      IF(var=='tpot') THEN; champ=champ*cpp/workvar !--- Potential temperature
     313      ELSE;                 champ=champ*.01*workvar !--- Relative humidity
     314        WHERE(champ<0.) champ=1.0E-10
     315      END IF
     316      DO il=1,lml
     317        xppn = SUM(aire(:,1  )*champ(:,1  ,il))/apoln
     318        xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
     319        champ(:,1  ,il) = xppn
     320        champ(:,jml,il) = xpps
     321      END DO
     322  END SELECT
    327323
    328324END SUBROUTINE startget_dyn3d
    … …  
    768764
    769765!#endif
    770 ! of #ifdef CPP_EARTH
     766! of ifdef CPP_EARTH
    771767
    772768END MODULE etat0dyn
Note: See TracChangeset for help on using the changeset viewer.