Changeset 4419 for LMDZ6/trunk/libf/dyn3d

Timestamp:

Feb 5, 2023, 9:32:11 PM (22 months ago)

Author:

fhourdin

Message:

Lecture possible du relief dans iniacademic

File:

• LMDZ6/trunk/libf/dyn3d/iniacademic.F90 (modified) (3 diffs)

LMDZ6/trunk/libf/dyn3d/iniacademic.F90

@@ -22 +22 @@
   USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
   USE readTracFiles_mod, ONLY: addPhase
+  use netcdf, only : NF90_NOWRITE,NF90_OPEN,NF90_NOERR,NF90_INQ_VARID
+  use netcdf, only : NF90_CLOSE, NF90_GET_VAR
+
 
   !   Author:    Frederic Hourdin      original: 15/01/93
@@ -63 +66 @@
   real tetajl(jjp1,llm)
   INTEGER i,j,l,lsup,ij, iq, iName, iPhase, iqParent
+
+  integer :: nid_relief,varid,ierr
+  real, dimension(iip1,jjp1) :: relief
 
   REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
@@ -118 +124 @@
   CALL inifilr
 
+
+  !------------------------------------------------------------------
   ! Initialize pressure and mass field if read_start=.false.
+  !------------------------------------------------------------------
+
   IF (.NOT. read_start) THEN
+
+     !------------------------------------------------------------------
+     ! Lecture eventuelle d'un fichier de relief interpollee sur la grille
+     ! du modele.
+     ! On suppose que le fichier relief_in.nc est stoké sur une grille
+     ! iim*jjp1
+     ! Facile a créer à partir de la commande
+     ! cdo remapcon,fichier_output_phys.nc Relief.nc relief_in.nc
+     !------------------------------------------------------------------
+
+     relief=0.
+     ierr = NF90_OPEN ('relief_in.nc', NF90_NOWRITE,nid_relief)
+     if (ierr.EQ.NF90_NOERR) THEN
+         ierr=NF90_INQ_VARID(nid_relief,'RELIEF',varid)
+         if (ierr==NF90_NOERR) THEN
+              ierr=NF90_GET_VAR(nid_relief,varid,relief(1:iim,1:jjp1))
+              relief(iip1,:)=relief(1,:)
+         else
+              CALL abort_gcm ('iniacademic','variable RELIEF pas la',1)
+         endif
+     endif
+     ierr = NF90_CLOSE (nid_relief)
+
+     !------------------------------------------------------------------
+     ! Initialisation du geopotentiel au sol et de la pression
+     !------------------------------------------------------------------
+
+     print*,'relief=',minval(relief),maxval(relief),'g=',g
+     do j=1,jjp1
+        do i=1,iip1
+           phis((j-1)*iip1+i)=g*relief(i,j)
+        enddo
+     enddo
+     print*,'phis=',minval(phis),maxval(phis),'g=',g
+
+     ! ground geopotential
+     !phis(:)=0.
      ps(:)=preff
-
-     ! ground geopotential
-     phis(:)=0.
      CALL pression ( ip1jmp1, ap, bp, ps, p       )