Changeset 1454 for LMDZ5/trunk/libf/dyn3dpar/iniacademic.F

Timestamp:

Nov 18, 2010, 1:01:24 PM (14 years ago)

Author:

Laurent Fairhead

Message:

Merge of LMDZ5V1.0-dev branch r1453 into LMDZ5 trunk r1434

---

Fusion entre la version r1453 de la branche de développement LMDZ5V1.0-dev
et le tronc LMDZ5 (r1434)

Location:

LMDZ5/trunk

Files:

• . (modified) (1 prop)
• libf/dyn3dpar/iniacademic.F (modified) (6 diffs)

LMDZ5/trunk/libf/dyn3dpar/iniacademic.F

@@ -8 +8 @@
       USE filtreg_mod
       USE infotrac, ONLY : nqtot
-      USE control_mod
+      USE control_mod, ONLY: day_step,planet_type
+#ifdef CPP_IOIPSL
+      USE IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin
+      USE ioipsl_getincom
+#endif
+      USE Write_Field
  
 
@@ -70 +77 @@
       REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
       REAL phi(ip1jmp1,llm)                  ! geopotentiel
-      REAL ddsin,tetarappelj,tetarappell,zsig
+      REAL ddsin,tetastrat,zsig,tetapv,w_pv  ! variables auxiliaires
       real tetajl(jjp1,llm)
       INTEGER i,j,l,lsup,ij
 
+      REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
+      REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
+      LOGICAL ok_geost             ! Initialisation vent geost. ou nul
+      LOGICAL ok_pv                ! Polar Vortex
+      REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex 
+      
       real zz,ran1
       integer idum
@@ -82 +95 @@
 ! 1. Initializations for Earth-like case
 ! --------------------------------------
-      if (planet_type=="earth") then
-c
+c
+        ! initialize planet radius, rotation rate,...
+        call conf_planete
+
         time_0=0.
-        day_ref=0
+        day_ref=1
         annee_ref=0
 
         im         = iim
         jm         = jjm
-        day_ini    = 0
-        omeg       = 4.*asin(1.)/86400.
-        rad    = 6371229.
-        g      = 9.8
-        daysec = 86400.
+        day_ini    = 1
         dtvr    = daysec/REAL(day_step)
         zdtvr=dtvr
-        kappa  = 0.2857143
-        cpp    = 1004.70885
-        preff     = 101325.
-        pa        =  50000.
         etot0      = 0.
         ptot0      = 0.
@@ -120 +127 @@
           if (.not.read_start) then
             phis(:)=0.
-            q(:,:,1)=1.e-10
-            q(:,:,2)=1.e-15
-            q(:,:,3:nqtot)=0.
+            q(:,:,:)=0
             CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
           endif
@@ -129 +134 @@
         if (iflag_phys.eq.2) then
           ! initializations for the academic case
-          ps(:)=1.e5
-          phis(:)=0.
-c---------------------------------------------------------------------
-
-          taurappel=10.*daysec
-
-c---------------------------------------------------------------------
-c   Calcul de la temperature potentielle :
-c   --------------------------------------
-
+          
+!         if (planet_type=="earth") then
+
+          ! 1. local parameters
+          ! by convention, winter is in the southern hemisphere
+          ! Geostrophic wind or no wind?
+          ok_geost=.TRUE.
+          CALL getin('ok_geost',ok_geost)
+          ! Constants for Newtonian relaxation and friction
+          k_f=1.                !friction 
+          CALL getin('k_j',k_f)
+          k_f=1./(daysec*k_f)
+          k_c_s=4.  !cooling surface
+          CALL getin('k_c_s',k_c_s)
+          k_c_s=1./(daysec*k_c_s)
+          k_c_a=40. !cooling free atm
+          CALL getin('k_c_a',k_c_a)
+          k_c_a=1./(daysec*k_c_a)
+          ! Constants for Teta equilibrium profile
+          teta0=315.     ! mean Teta (S.H. 315K)
+          CALL getin('teta0',teta0)
+          ttp=200.       ! Tropopause temperature (S.H. 200K)
+          CALL getin('ttp',ttp)
+          eps=0.         ! Deviation to N-S symmetry(~0-20K)
+          CALL getin('eps',eps)
+          delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
+          CALL getin('delt_y',delt_y)
+          delt_z=10.     ! Vertical Gradient (S.H. 10K)
+          CALL getin('delt_z',delt_z)
+          ! Polar vortex
+          ok_pv=.false.
+          CALL getin('ok_pv',ok_pv)
+          phi_pv=-50.            ! Latitude of edge of vortex
+          CALL getin('phi_pv',phi_pv)
+          phi_pv=phi_pv*pi/180.
+          dphi_pv=5.             ! Width of the edge
+          CALL getin('dphi_pv',dphi_pv)
+          dphi_pv=dphi_pv*pi/180.
+          gam_pv=4.              ! -dT/dz vortex (in K/km)
+          CALL getin('gam_pv',gam_pv)
+          
+          ! 2. Initialize fields towards which to relax
+          ! Friction
+          knewt_g=k_c_a
           DO l=1,llm
-            zsig=ap(l)/preff+bp(l)
-            if (zsig.gt.0.3) then
-             lsup=l
-             tetarappell=1./8.*(-log(zsig)-.5)
-             DO j=1,jjp1
-             ddsin=sin(rlatu(j))-sin(pi/20.)
-             tetajl(j,l)=300.*(1+1./18.*(1.-3.*ddsin*ddsin)+tetarappell)
-             ENDDO
-            else
-c   Choix isotherme au-dessus de 300 mbar
-             do j=1,jjp1
-               tetajl(j,l)=tetajl(j,lsup)*(0.3/zsig)**kappa
-             enddo
-            endif ! of if (zsig.gt.0.3)
+           zsig=presnivs(l)/preff
+           knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
+           kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
+          ENDDO
+          DO j=1,jjp1
+            clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
+          ENDDO
+          
+          ! Potential temperature 
+          DO l=1,llm
+           zsig=presnivs(l)/preff
+           tetastrat=ttp*zsig**(-kappa)
+           tetapv=tetastrat
+           IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
+               tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
+           ENDIF
+           DO j=1,jjp1
+             ! Troposphere
+             ddsin=sin(rlatu(j))
+             tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin
+     &           -delt_z*(1.-ddsin*ddsin)*log(zsig)
+             ! Profil stratospherique isotherme (+vortex)
+             w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
+             tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
+             tetajl(j,l)=MAX(tetajl(j,l),tetastrat)  
+           ENDDO
           ENDDO ! of DO l=1,llm
+ 
+!          CALL writefield('theta_eq',tetajl)
 
           do l=1,llm
@@ -165 +218 @@
           enddo
 
-c       call dump2d(jjp1,llm,tetajl,'TEQ   ')
-
-          CALL pression ( ip1jmp1, ap, bp, ps, p       )
-          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-          CALL massdair(p,masse)
-
-c  intialisation du vent et de la temperature
-          teta(:,:)=tetarappel(:,:)
-          CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
-          call ugeostr(phi,ucov)
-          vcov=0.
-          q(:,:,1   )=1.e-10
-          q(:,:,2   )=1.e-15
-          q(:,:,3:nqtot)=0.
-
-
-c   perturbation aleatoire sur la temperature
-          idum  = -1
-          zz = ran1(idum)
-          idum  = 0
-          do l=1,llm
-            do ij=iip2,ip1jm
-              teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
+
+!         else
+!          write(lunout,*)"iniacademic: planet types other than earth",
+!     &                   " not implemented (yet)."
+!          stop
+!         endif ! of if (planet_type=="earth")
+
+          ! 3. Initialize fields (if necessary)
+          IF (.NOT. read_start) THEN
+            ! surface pressure
+            ps(:)=preff
+            ! ground geopotential
+            phis(:)=0.
+            
+            CALL pression ( ip1jmp1, ap, bp, ps, p       )
+            CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+            CALL massdair(p,masse)
+
+            ! bulk initialization of temperature
+            teta(:,:)=tetarappel(:,:)
+            
+            ! geopotential
+            CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
+            
+            ! winds
+            if (ok_geost) then
+              call ugeostr(phi,ucov)
+            else
+              ucov(:,:)=0.
+            endif
+            vcov(:,:)=0.
+            
+            ! bulk initialization of tracers
+            if (planet_type=="earth") then
+              ! Earth: first two tracers will be water
+              do i=1,nqtot
+                if (i.eq.1) q(:,:,i)=1.e-10
+                if (i.eq.2) q(:,:,i)=1.e-15
+                if (i.gt.2) q(:,:,i)=0.
+              enddo
+            else
+              q(:,:,:)=0
+            endif ! of if (planet_type=="earth")
+
+            ! add random perturbation to temperature
+            idum  = -1
+            zz = ran1(idum)
+            idum  = 0
+            do l=1,llm
+              do ij=iip2,ip1jm
+                teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
+              enddo
             enddo
-          enddo
-
-          do l=1,llm
-            do ij=1,ip1jmp1,iip1
-              teta(ij+iim,l)=teta(ij,l)
+
+            ! maintain periodicity in longitude
+            do l=1,llm
+              do ij=1,ip1jmp1,iip1
+                teta(ij+iim,l)=teta(ij,l)
+              enddo
             enddo
-          enddo
-
-
-
-c     PRINT *,' Appel test_period avec tetarappel '
-c     CALL  test_period ( ucov,vcov,tetarappel,q,p,phis )
-c     PRINT *,' Appel test_period avec teta '
-c     CALL  test_period ( ucov,vcov,teta,q,p,phis )
-
-c   initialisation d'un traceur sur une colonne
-          j=jjp1*3/4
-          i=iip1/2
-          ij=(j-1)*iip1+i
-          q(ij,:,3)=1.
+
+          ENDIF ! of IF (.NOT. read_start)
         endif ! of if (iflag_phys.eq.2)
         
-      else
-        write(lunout,*)"iniacademic: planet types other than earth",
-     &                 " not implemented (yet)."
-        stop
-      endif ! of if (planet_type=="earth")
-      return
       END
 c-----------------------------------------------------------------------