source: codes/icosagcm/trunk/src/diffff @ 150

Last change on this file since 150 was 149, checked in by sdubey, 11 years ago
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 size: 15.6 KB
Line 
1Index: time.f90
2===================================================================
3--- time.f90    (revision 142)
4+++ time.f90    (working copy)
5@@ -10,11 +10,23 @@
6   REAL(rstd),SAVE :: write_period
7   INTEGER,SAVE    :: itau_out, itau_adv, itau_dissip, itau_physics, itaumax
8   
9+  INTEGER,SAVE :: day_step,ndays
10+  REAL(rstd),SAVE :: jD_ref,jH_ref
11+  INTEGER,SAVE :: day_ini,day_end,annee_ref,day_ref
12+  REAL(rstd),SAVE::start_time
13+  CHARACTER(LEN=255) :: time_style
14+  INTEGER,SAVE:: an, mois, jour
15+  REAL(rstd),SAVE:: heure
16+  CHARACTER (LEN=10):: calend
17+
18   PUBLIC create_time_counter_header, update_time_counter, close_time_counter, init_time,  &
19-         dt, write_period, itau_out, itau_adv, itau_dissip, itau_physics, itaumax
20+         dt, write_period, itau_out, itau_adv, itau_dissip, itau_physics, itaumax, &
21+day_step,ndays,jD_ref,jH_ref,day_ini,day_end,annee_ref,day_ref,an, mois, jour,heure, &
22+            calend,time_style
23 
24 
25 
26+
27 CONTAINS
28   
29   SUBROUTINE init_time
30@@ -23,22 +35,18 @@
31   USE mpipara
32   IMPLICIT NONE
33   REAL(rstd) :: run_length
34- 
35+
36+
37+     time_style='dcmip'
38+   CALL getin('time_style',time_style)
39+
40+   IF (TRIM(time_style)=='dcmip')  Then
41     dt=90.
42     CALL getin('dt',dt)
43 
44     itaumax=100
45     CALL getin('itaumax',itaumax)
46 
47-    itau_adv=1
48-    CALL getin('itau_adv',itau_adv)
49-
50-    itau_dissip=1
51-    CALL getin('itau_dissip',itau_dissip)
52-
53-    itau_physics=1
54-    CALL getin('itau_physics',itau_physics)
55-
56     run_length=dt*itaumax
57     CALL getin('run_length',run_length)
58     itaumax=run_length/dt
59@@ -53,6 +61,17 @@
60     write_period=write_period/scale_factor
61     itau_out=FLOOR(.5+write_period/dt)
62     IF (is_mpi_root) PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out
63+    ENDIF
64+
65+    itau_adv=1
66+    CALL getin('itau_adv',itau_adv)
67+
68+    itau_dissip=1
69+    CALL getin('itau_dissip',itau_dissip)
70+
71+    itau_physics=1
72+    CALL getin('itau_physics',itau_physics)
73+
74     
75     CALL create_time_counter_header
76     
77Index: caldyn_gcm.f90
78===================================================================
79--- caldyn_gcm.f90      (revision 142)
80+++ caldyn_gcm.f90      (working copy)
81@@ -10,7 +10,7 @@
82   TYPE(t_field),POINTER :: f_buf_i(:), f_buf_ulon(:), f_buf_ulat(:), f_buf_u3d(:)
83   TYPE(t_field),POINTER :: f_buf_v(:), f_buf_s(:), f_buf_p(:)
84 
85-  PUBLIC init_caldyn, caldyn, write_output_fields
86+  PUBLIC init_caldyn, caldyn, write_output_fields,un2ulonlat
87 
88   INTEGER :: caldyn_hydrostat, caldyn_conserv
89 
90@@ -708,17 +708,17 @@
91     str_pression=ADJUSTL(str_pression)
92     
93     CALL writefield("ps",f_ps)
94-    CALL writefield("dps",f_dps)
95-    CALL writefield("phis",f_phis)
96-    CALL vorticity(f_u,f_buf_v)
97-    CALL writefield("vort",f_buf_v)
98+!    CALL writefield("dps",f_dps)
99+!    CALL writefield("phis",f_phis)
100+!    CALL vorticity(f_u,f_buf_v)
101+!    CALL writefield("vort",f_buf_v)
102 
103-    CALL w_omega(f_ps, f_u, f_buf_i)
104-    CALL writefield('omega', f_buf_i)
105-    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
106-      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
107-      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
108-    ENDIF
109+!    CALL w_omega(f_ps, f_u, f_buf_i)
110+!    CALL writefield('omega', f_buf_i)
111+!    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
112+!      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
113+!      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
114+!    ENDIF
115     
116     ! Temperature
117     CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ;
118@@ -753,10 +753,10 @@
119     
120     ! geopotential
121     CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i)
122-    CALL writefield("p",f_buf_p)
123-    CALL writefield("phi",f_buf_i)
124-    CALL writefield("theta",f_buf1_i) ! potential temperature
125-    CALL writefield("pk",f_buf2_i) ! Exner pressure
126+!    CALL writefield("p",f_buf_p)
127+!    CALL writefield("phi",f_buf_i)
128+     CALL writefield("theta",f_buf1_i) ! potential temperature
129+!    CALL writefield("pk",f_buf2_i) ! Exner pressure
130   
131   
132   END SUBROUTINE write_output_fields
133Index: advect_tracer.f90
134===================================================================
135--- advect_tracer.f90   (revision 142)
136+++ advect_tracer.f90   (working copy)
137@@ -57,7 +57,7 @@
138     CALL transfert_request(f_q,req_i1)
139     CALL transfert_request(f_rhodz,req_i1)
140       
141-    IF (is_mpi_root) PRINT *, 'Advection scheme'
142+!    IF (is_mpi_root) PRINT *, 'Advection scheme'
143 
144 !    DO ind=1,ndomain
145 !       CALL swap_dimensions(ind)
146Index: dissip_gcm.f90
147===================================================================
148--- dissip_gcm.f90      (revision 142)
149+++ dissip_gcm.f90      (working copy)
150@@ -93,7 +93,8 @@
151       rayleigh_friction_type=2
152       IF (is_mpi_root) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2'
153    CASE DEFAULT
154-      IF (is_mpi_root) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), ' in dissip_gcm.f90/init_dissip'
155+      IF (is_mpi_root) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), &
156+        ' in dissip_gcm.f90/init_dissip'
157       STOP
158    END SELECT
159 
160Index: physics.f90
161===================================================================
162--- physics.f90 (revision 142)
163+++ physics.f90 (working copy)
164@@ -7,7 +7,8 @@
165 
166   SUBROUTINE init_physics
167   USE icosa
168-  USE physics_dcmip_mod, init_physics_dcmip=>init_physics
169+  USE physics_dcmip_mod,init_physics_dcmip=>init_physics
170+  USE physics_dry_mod
171   IMPLICIT NONE
172     
173     CALL getin("physics",physics_type)
174@@ -17,36 +18,57 @@
175     
176       CASE ('dcmip')
177         CALL init_physics_dcmip
178+
179+      CASE ('lmd')
180+       CALL init_physics_dry
181       
182       CASE DEFAULT
183-         PRINT*, 'Bad selector for variable physics <',physics_type, &
184+         PRINT*, 'Bad selector for variable physics init <',physics_type, &
185               '> options are <none>, <dcmip>,'
186-         STOP
187+
188     END SELECT
189     
190   END SUBROUTINE init_physics
191 
192-  SUBROUTINE physics(it,f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
193+  SUBROUTINE physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
194   USE icosa
195+  USE physics_dry_mod
196   USE physics_dcmip_mod, physics_dcmip=>physics
197+  USE etat0_mod
198+  USE etat0_heldsz_mod
199   IMPLICIT NONE
200     INTEGER, INTENT(IN)   :: it
201+    REAL(rstd),INTENT(IN)::jD_cur,jH_cur
202     TYPE(t_field),POINTER :: f_phis(:)
203     TYPE(t_field),POINTER :: f_ps(:)
204     TYPE(t_field),POINTER :: f_theta_rhodz(:)
205     TYPE(t_field),POINTER :: f_ue(:)
206     TYPE(t_field),POINTER :: f_q(:)
207+    LOGICAL:: firstcall,lastcall
208     
209     SELECT CASE(TRIM(physics_type))
210       CASE ('none')
211+
212+       SELECT CASE(TRIM(etat0_type))
213+       CASE('heldsz')
214+       CALL transfert_request(f_ps,req_i1)
215+       CALL transfert_request(f_theta_rhodz,req_i1)
216+       CALL transfert_request(f_ue,req_e1)
217+       CALL held_saurez(f_ps,f_theta_rhodz,f_ue)
218+       CASE DEFAULT
219+       PRINT*,"NO PHYSICAL PACAKAGE USED"
220+       END SELECT
221     
222       CASE ('dcmip')
223         CALL physics_dcmip(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
224+
225+      CASE ('dry')
226+        CALL physics_dry(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
227       
228       CASE DEFAULT
229          PRINT*, 'Bad selector for variable physics <',physics_type, &
230               '> options are <none>, <dcmip>,'
231-         STOP
232+  STOP
233     END SELECT
234     
235   END SUBROUTINE physics
236Index: etat0.f90
237===================================================================
238--- etat0.f90   (revision 142)
239+++ etat0.f90   (working copy)
240@@ -1,4 +1,5 @@
241 MODULE etat0_mod
242+    CHARACTER(len=255),SAVE :: etat0_type
243 
244 CONTAINS
245   
246@@ -11,6 +12,10 @@
247     USE etat0_dcmip3_mod, ONLY : etat0_dcmip3=>etat0 
248     USE etat0_dcmip4_mod, ONLY : etat0_dcmip4=>etat0 
249     USE etat0_dcmip5_mod, ONLY : etat0_dcmip5=>etat0 
250+    USE etat0_heldsz_mod, ONLY : etat0_heldsz=>etat0 
251+    USE dynetat0_gcm_mod, ONLY : dynetat0_start=>etat0
252+    USE dynetat0_hz_mod,  ONLY : dynetat0_hz=>etat0
253+
254     IMPLICIT NONE
255     TYPE(t_field),POINTER :: f_ps(:)
256     TYPE(t_field),POINTER :: f_phis(:)
257@@ -18,7 +23,6 @@
258     TYPE(t_field),POINTER :: f_u(:)
259     TYPE(t_field),POINTER :: f_q(:)
260     
261-    CHARACTER(len=255) :: etat0_type
262     etat0_type='jablonowsky06'
263     CALL getin("etat0",etat0_type)
264     
265@@ -27,6 +31,9 @@
266        CALL etat0_jablonowsky06(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
267     CASE ('academic')
268        CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
269+    CASE ('heldsz')
270+       print*,"heldsz test case"
271+       CALL etat0_heldsz(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
272     CASE ('dcmip1')
273        CALL etat0_dcmip1(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
274     CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear')
275@@ -37,6 +44,12 @@
276        CALL etat0_dcmip4(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
277      CASE ('dcmip5')
278        CALL etat0_dcmip5(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
279+     CASE ('readnf_start')
280+         print*,"readnf_start used"   
281+       CALL dynetat0_start(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
282+       CASE ('readnf_hz')
283+         print*,"readnf_hz used"
284+       CALL dynetat0_hz(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
285    CASE DEFAULT
286        PRINT*, 'Bad selector for variable etat0 <',etat0_type, &
287             '> options are <jablonowsky06>, <academic>, <dcmip[1-4]> '
288Index: timeloop_gcm.f90
289===================================================================
290--- timeloop_gcm.f90    (revision 142)
291+++ timeloop_gcm.f90    (working copy)
292@@ -17,7 +17,11 @@
293   USE advect_tracer_mod
294   USE physics_mod
295   USE mpipara
296- 
297+!-------------------------------
298+  USE IOIPSL
299+  USE maxicosa
300+  USE check_conserve_mod
301+!------------------------------
302   IMPLICIT NONE
303   TYPE(t_field),POINTER :: f_phis(:)
304 !  TYPE(t_field),POINTER :: f_theta(:)
305@@ -46,6 +50,10 @@
306   INTEGER :: it,i,j,n, nb_stage, stage, matsuno_period, scheme
307   CHARACTER(len=255) :: scheme_name
308   LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
309+  CHARACTER(len=7) :: first
310+  REAL(rstd),SAVE :: jD_cur, jH_cur
311+  REAL(rstd),SAVE :: start_time
312+
313 !  INTEGER :: itaumax
314 !  REAL(rstd) ::write_period
315 !  INTEGER    :: itau_out
316@@ -84,6 +92,58 @@
317   CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time
318   CALL allocate_field(f_wfluxt,field_t,type_real,llm+1)
319 
320+!----------------------------------------------------
321+  IF (TRIM(time_style)=='lmd')  Then
322+
323+   day_step=180
324+   CALL getin('day_step',day_step)
325+
326+   ndays=1
327+   CALL getin('ndays',ndays)
328+
329+   dt = daysec/REAL(day_step)
330+   itaumax = ndays*day_step
331+
332+   calend = 'earth_360d'
333+   CALL getin('calend', calend)
334+
335+   day_ini = 0
336+   CALL getin('day_ini',day_ini)
337+
338+   day_end = 0
339+   CALL getin('day_end',day_end)
340+
341+   annee_ref = 1998
342+   CALL getin('annee_ref',annee_ref)
343+
344+   start_time = 0
345+   CALL getin('start_time',start_time)
346+
347+   write_period=0
348+   CALL getin('write_period',write_period)
349+     
350+   write_period=write_period/scale_factor
351+   itau_out=FLOOR(write_period/dt)
352+   
353+   PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out
354+
355+  mois = 1 ; heure = 0.
356+  call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
357+  jH_ref = jD_ref - int(jD_ref)
358+  jD_ref = int(jD_ref)
359+ 
360+  CALL ioconf_startdate(INT(jD_ref),jH_ref)
361+  write(*,*)'annee_ref, mois, day_ref, heure, jD_ref'
362+  write(*,*)annee_ref, mois, day_ref, heure, jD_ref
363+  write(*,*)"ndays,day_step,itaumax,dt======>"
364+  write(*,*)ndays,day_step,itaumax,dt
365+  call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
366+  write(*,*)'jD_ref+jH_ref,an, mois, jour, heure'
367+  write(*,*)jD_ref+jH_ref,an, mois, jour, heure
368+  day_end = day_ini + ndays
369+ END IF
370+!----------------------------------------------------
371+
372   scheme_name='runge_kutta'
373   CALL getin('scheme',scheme_name)
374 
375@@ -122,7 +182,6 @@
376 !  CALL allocate_field(f_dum2,field_u,type_real,llm)
377 !  CALL allocate_field(f_dtheta_rhodzm1,field_t,type_real,llm)
378 !  CALL allocate_field(f_dtheta_rhodzm2,field_t,type_real,llm)
379-
380 !  CALL allocate_field(f_theta,field_t,type_real,llm)
381 !  CALL allocate_field(f_dtheta,field_t,type_real,llm)
382 
383@@ -132,7 +191,11 @@
384   CALL init_advect_tracer
385   CALL init_physics
386   
387-  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
388+!============================ SET 0 to all field
389+  CALL initial0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
390+! CALL dynetat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
391+  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
392+!================================================
393   CALL writefield("phis",f_phis,once=.TRUE.)
394   CALL transfert_request(f_q,req_i1)
395 
396@@ -153,17 +216,18 @@
397      rhodz=f_rhodz(ind); ps=f_ps(ind)
398      CALL compute_rhodz(.FALSE., ps, rhodz)   
399   END DO
400+
401   
402   DO it=0,itaumax
403 
404-    IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
405     IF (mod(it,itau_out)==0 ) THEN
406-      CALL writefield("q",f_q)
407+!      IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
408       CALL update_time_counter(dt*it)
409+      CALL compute_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 
410     ENDIF
411-   
412-    CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
413 
414+      CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
415+
416     DO stage=1,nb_stage
417        CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
418             f_phis,f_ps,f_theta_rhodz,f_u, f_q, &
419@@ -177,11 +241,11 @@
420           CALL  leapfrog_matsuno_scheme(stage)
421           
422           !      CASE ('leapfrog')
423-          !        CALL leapfrog_scheme
424+          !      CALL leapfrog_scheme
425           !
426           !      CASE ('adam_bashforth')
427-          !        CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
428-          !        CALL adam_bashforth_scheme
429+          !      CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
430+          !      CALL adam_bashforth_scheme
431        CASE DEFAULT
432           STOP
433        END SELECT
434@@ -192,7 +256,7 @@
435 
436     IF(MOD(it+1,itau_adv)==0) THEN
437        CALL transfert_request(f_wfluxt,req_i1) ! FIXME
438-!       CALL transfert_request(f_hfluxt,req_e1) ! FIXME
439+!      CALL transfert_request(f_hfluxt,req_e1) ! FIXME
440 
441        CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
442        fluxt_zero=.TRUE.
443@@ -207,11 +271,17 @@
444           CALL swap_geometry(ind)
445           rhodz=f_rhodz(ind); ps=f_ps(ind); dps=f_dps(ind);
446           wflux=f_wflux(ind); wfluxt=f_wfluxt(ind)
447-          CALL compute_rhodz(.FALSE., ps, rhodz)   
448+!!!!!          CALL compute_rhodz(.FALSE., ps, rhodz)   
449        END DO
450-
451     END IF
452-
453+!----------------------------------------------------
454+  jD_cur = jD_ref + day_ini - day_ref + it/day_step
455+  jH_cur = jH_ref + start_time + mod(it,day_step)/float(day_step)
456+  jD_cur = jD_cur + int(jH_cur)
457+  jH_cur = jH_cur - int(jH_cur)
458+!      print*,"Just b4 phys"
459+    CALL physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
460+!----------------------------------------------------
461 !    CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
462     ENDDO
463 
464@@ -399,7 +469,7 @@
465        DO j=jj_begin-dd,jj_end+dd
466           DO i=ii_begin-dd,ii_end+dd
467              ij=(j-1)*iim+i
468-             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g
469+             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g
470              IF(comp) THEN
471                 rhodz(ij,l) = m
472              ELSE
473@@ -414,7 +484,7 @@
474           PRINT *, 'Discrepancy between ps and rhodz detected', err
475           STOP
476        ELSE
477-!          PRINT *, 'No discrepancy between ps and rhodz detected'
478+          PRINT *, 'No discrepancy between ps and rhodz detected'
479        END IF
480     END IF
481 
482Index: icosa_gcm.f90
483===================================================================
484--- icosa_gcm.f90       (revision 142)
485+++ icosa_gcm.f90       (working copy)
486@@ -24,6 +24,7 @@
487   CALL compute_domain
488   CALL init_transfert
489   CALL init_writefield
490+
491 !  CALL allocate_field(sum_ne,field_T,type_real)
492 !  CALL allocate_field_glo(sum_ne_glo,field_T,type_real)
493 ! 
Note: See TracBrowser for help on using the repository browser.