Changeset 149 for codes/icosagcm/trunk/src/timeloop_gcm.f90

Timestamp:

04/03/13 12:05:12 (11 years ago)

Author:

sdubey

Message:

Added few new routines to read NC files and compute diagnostics to r145.

Few routines of dry physics including radiation module, surface process and convective adjustment in new routine phyparam.f90. dynetat to read start files for dynamics.

check_conserve routine to compute conservation of quatities like mass, energy etc.etat0_heldsz.f90 for held-suarez test case initial conditions.

new Key time_style=lmd or dcmip to use day_step, ndays like in LMDZ

File:

• codes/icosagcm/trunk/src/timeloop_gcm.f90 (modified) (13 diffs)

codes/icosagcm/trunk/src/timeloop_gcm.f90

@@ -19 +19 @@
   USE physics_mod
   USE mpipara
+  USE IOIPSL
+  USE maxicosa
+  USE check_conserve_mod
   USE trace
   USE transfert_mod
@@ -43 +46 @@
   REAL(rstd),POINTER :: dtheta_rhodz(:,:),dtheta_rhodzm1(:,:),dtheta_rhodzm2(:,:)
   REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
-
 !  REAL(rstd) :: dt, run_length
   INTEGER :: ind
@@ -49 +51 @@
   CHARACTER(len=255) :: scheme_name
   LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
+  CHARACTER(len=7) :: first 
+  REAL(rstd),SAVE :: jD_cur, jH_cur
+  REAL(rstd),SAVE :: start_time
   LOGICAL, PARAMETER :: check=.FALSE.
 !  INTEGER :: itaumax
@@ -88 +93 @@
   CALL allocate_field(f_wfluxt,field_t,type_real,llm+1)
 
+!----------------------------------------------------
+  IF (TRIM(time_style)=='lmd')  Then
+
+   day_step=180
+   CALL getin('day_step',day_step)
+
+   ndays=1
+   CALL getin('ndays',ndays)
+
+   dt = daysec/REAL(day_step)
+   itaumax = ndays*day_step
+
+   calend = 'earth_360d'
+   CALL getin('calend', calend)
+
+   day_ini = 0
+   CALL getin('day_ini',day_ini)
+
+   day_end = 0
+   CALL getin('day_end',day_end)
+
+   annee_ref = 1998
+   CALL getin('annee_ref',annee_ref)
+
+   start_time = 0
+   CALL getin('start_time',start_time) 
+
+   write_period=0
+   CALL getin('write_period',write_period)
+      
+   write_period=write_period/scale_factor
+   itau_out=FLOOR(write_period/dt)
+   
+   PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out 
+
+  mois = 1 ; heure = 0.
+  call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
+  jH_ref = jD_ref - int(jD_ref) 
+  jD_ref = int(jD_ref) 
+  
+  CALL ioconf_startdate(INT(jD_ref),jH_ref) 
+  write(*,*)'annee_ref, mois, day_ref, heure, jD_ref'
+  write(*,*)annee_ref, mois, day_ref, heure, jD_ref
+  write(*,*)"ndays,day_step,itaumax,dt======>"
+  write(*,*)ndays,day_step,itaumax,dt
+  call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
+  write(*,*)'jD_ref+jH_ref,an, mois, jour, heure'
+  write(*,*)jD_ref+jH_ref,an, mois, jour, heure
+  day_end = day_ini + ndays 
+ END IF 
+!----------------------------------------------------
+
   scheme_name='runge_kutta'
   CALL getin('scheme',scheme_name)
@@ -126 +183 @@
 !  CALL allocate_field(f_dtheta_rhodzm1,field_t,type_real,llm)
 !  CALL allocate_field(f_dtheta_rhodzm2,field_t,type_real,llm)
-
 !  CALL allocate_field(f_theta,field_t,type_real,llm)
 !  CALL allocate_field(f_dtheta,field_t,type_real,llm)
@@ -135 +191 @@
   CALL init_advect_tracer
   CALL init_physics
-
-   
-  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
+!========================================= INITIALIZATION 
+! CALL dynetat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
+  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
   CALL writefield("phis",f_phis,once=.TRUE.)
   CALL transfert_request(f_q,req_i1) 
@@ -158 +214 @@
      CALL compute_rhodz(.FALSE., ps, rhodz)   
   END DO
-  
+
   CALL transfert_request(f_phis,req_i0) 
+  CALL transfert_request(f_phis,req_i1) 
   CALL transfert_request(f_phis,req_i1) 
 
@@ -170 +227 @@
     ENDIF
     
-    IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
     IF (mod(it,itau_out)==0 ) THEN
-      CALL writefield("q",f_q)
+!      IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
       CALL update_time_counter(dt*it)
+      CALL compute_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)  
     ENDIF
-    
-    CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
+
+      CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
 
     DO stage=1,nb_stage
@@ -189 +246 @@
        CASE (mlf)
           CALL  leapfrog_matsuno_scheme(stage)
-          
           !      CASE ('leapfrog')
-          !        CALL leapfrog_scheme
+          !      CALL leapfrog_scheme
           !
           !      CASE ('adam_bashforth')
-          !        CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
-          !        CALL adam_bashforth_scheme
+          !      CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
+          !      CALL adam_bashforth_scheme
        CASE DEFAULT
           STOP
@@ -208 +264 @@
     IF(MOD(it+1,itau_adv)==0) THEN
 !       CALL transfert_request(f_wfluxt,req_i1) ! FIXME
-!       CALL transfert_request(f_hfluxt,req_e1) ! FIXME
+!      CALL transfert_request(f_hfluxt,req_e1) ! FIXME
 
        CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
@@ -227 +283 @@
          END DO
        ENDIF
-
     END IF
-
+!----------------------------------------------------
+  jD_cur = jD_ref + day_ini - day_ref + it/day_step
+  jH_cur = jH_ref + start_time + mod(it,day_step)/float(day_step)
+  jD_cur = jD_cur + int(jH_cur)
+  jH_cur = jH_cur - int(jH_cur)
+!       print*,"Just b4 phys"
+    CALL physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
+!----------------------------------------------------
 !    CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
     ENDDO
@@ -493 +555 @@
           DO i=ii_begin-dd,ii_end+dd
              ij=(j-1)*iim+i
-             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
+             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g
              IF(comp) THEN
                 rhodz(ij,l) = m
@@ -508 +570 @@
           STOP
        ELSE
-!          PRINT *, 'No discrepancy between ps and rhodz detected'
+          PRINT *, 'No discrepancy between ps and rhodz detected'
        END IF
     END IF