Changeset 4419 for LMDZ6/trunk/libf/dyn3dmem/iniacademic_loc.F90

Timestamp:

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

Author:

fhourdin

Message:

Lecture possible du relief dans iniacademic

File:

• LMDZ6/trunk/libf/dyn3dmem/iniacademic_loc.F90 (modified) (4 diffs)

LMDZ6/trunk/libf/dyn3dmem/iniacademic_loc.F90

@@ -23 +23 @@
   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
@@ -67 +70 @@
   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
@@ -126 +133 @@
     ! allocate global fields:
 !    allocate(vcov_glo(ip1jm,llm))
+
     allocate(ucov_glo(ip1jmp1,llm))
     allocate(teta_glo(ip1jmp1,llm))
@@ -135 +143 @@
      ps_glo(:)=preff
 
-     ! ground geopotential
-     phis_glo(:)=0.
+     !------------------------------------------------------------------
+     ! 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_glo((j-1)*iip1+i)=g*relief(i,j)
+        enddo
+     enddo
+     print*,'phis=',minval(phis),maxval(phis),'g=',g
+
      CALL pression ( ip1jmp1, ap, bp, ps_glo, p       )
      if (pressure_exner) then