source: CONFIG/UNIFORM/v7/IPSLCM7/SOURCES/DYNAMICO/trunk/time/timeloop_gcm.F90 @ 6362

Last change on this file since 6362 was 6362, checked in by aclsce, 15 months ago

Modified sponge parameters ICOLMDZ to use same treatment as LMDZ (iflag_sponge=2, mode_sponge=3).
Remind :
iflag_sponge !0 --> for no sponge

!1 --> for sponge over 4 topmost layers
!2 --> for sponge from top to ~1% of top layer pressure

mode_sponge !1 --> u and v relax towards 0

!2 --> u and v relax towards their zonal mean
!3 --> u,v and pot. temp. relax towards their zonal mean
File size: 18.8 KB
Line 
1MODULE timeloop_gcm_mod
2  USE profiling_mod
3  USE icosa
4  USE disvert_mod
5  USE trace
6  USE omp_para
7  USE euler_scheme_mod
8  USE explicit_scheme_mod
9  USE hevi_scheme_mod
10  IMPLICIT NONE
11  PRIVATE
12 
13  INTEGER, PARAMETER :: sync_it=10
14  TYPE(t_message),SAVE :: req_ps0, req_rhodz0, req_mass0, req_theta_rhodz0, req_u0, req_q0, req_W0, req_geopot0
15  LOGICAL, SAVE :: positive_theta
16  INTEGER :: itau_prof, id_timeloop, id_dyn, id_phys, id_dissip, id_adv, id_diags
17  PUBLIC :: init_timeloop, timeloop
18
19CONTAINS
20 
21  SUBROUTINE init_timeloop
22    USE dissip_gcm_mod
23    USE observable_mod
24    USE caldyn_mod
25    USE etat0_mod
26    USE guided_mod
27    USE advect_tracer_mod
28    USE check_conserve_mod
29    USE output_field_mod
30    USE theta2theta_rhodz_mod
31    USE sponge_mod
32
33    CHARACTER(len=255) :: def
34
35    CALL register_id('timeloop',id_timeloop)
36    CALL register_id('dyn',id_dyn)
37    CALL register_id('dissip',id_dissip)
38    CALL register_id('phys',id_phys)
39    CALL register_id('adv',id_adv)
40    CALL register_id('diags',id_diags)
41
42    CALL init_caldyn
43   
44!    IF (xios_output) itau_out=1
45    IF (.NOT. enable_io) itau_out=HUGE(itau_out)
46
47    itau_prof=1000
48    CALL getin('itau_profiling',itau_prof)
49
50    positive_theta=.FALSE.
51    CALL getin('positive_theta',positive_theta)
52    IF(positive_theta .AND. nqtot<1) THEN
53       PRINT *, 'nqtot must be >0 if positive_theta is .TRUE.'
54       STOP
55    END IF
56
57    def='ARK2.3'
58    CALL getin('time_scheme',def)
59    SELECT CASE (TRIM(def))
60    CASE('euler')
61       scheme_family=explicit
62       scheme=euler
63       nb_stage=1
64    CASE ('runge_kutta')
65       scheme_family=explicit
66       scheme=rk4
67       nb_stage=4
68    CASE ('RK2.5')
69       scheme_family=explicit
70       scheme=rk25
71       nb_stage=5
72    CASE ('leapfrog_matsuno')
73       scheme_family=explicit
74       scheme=mlf
75       matsuno_period=5
76       CALL getin('matsuno_period',matsuno_period)
77       nb_stage=matsuno_period+1
78    CASE('ARK2.3')
79       scheme_family=hevi
80       scheme=ark23
81       nb_stage=3
82       CALL set_coefs_ark23(dt)
83    CASE('ARK3.3')
84       scheme_family=hevi
85       scheme=ark33
86       nb_stage=3
87       CALL set_coefs_ark33(dt)
88    CASE ('none')
89       nb_stage=0
90    CASE default
91       PRINT*,'Bad selector for variable scheme : <', TRIM(def),             &
92            ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>,<RK2.5>,<ARK2.3>'
93       STOP
94    END SELECT
95   
96    IF (scheme_family /= hevi) THEN
97       CALL abort_acc("scheme_family /= hevi")
98    END IF
99
100    ! Time-independant orography
101    CALL allocate_field(f_phis,field_t,type_real,name='phis')
102    ! Model state at current time step
103    CALL allocate_field(f_ps,field_t,type_real, name='ps')
104    CALL allocate_field(f_mass,field_t,type_real,llm,name='mass')
105    CALL allocate_field(f_rhodz,field_t,type_real,llm,name='rhodz')
106    CALL allocate_field(f_theta_rhodz,field_t,type_real,llm,nqdyn,name='theta_rhodz')
107    CALL allocate_field(f_u,field_u,type_real,llm,name='u')
108    CALL allocate_field(f_geopot,field_t,type_real,llm+1,name='geopot')
109    CALL allocate_field(f_W,field_t,type_real,llm+1,name='W') ! used only if .not. hydrostatic
110    CALL allocate_field(f_q,field_t,type_real,llm,nqtot,'q')
111    ! Mass fluxes
112    CALL allocate_field(f_hflux,field_u,type_real,llm, ondevice=.TRUE.)    ! instantaneous mass fluxes
113    CALL allocate_field(f_hfluxt,field_u,type_real,llm,ondevice=.TRUE.)   ! mass "fluxes" accumulated in time
114    CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! vertical mass fluxes
115    CALL allocate_field(f_wfluxt,field_t,type_real,llm+1,name='wfluxt',ondevice=.TRUE.)
116   
117    SELECT CASE(scheme_family)
118    CASE(explicit)
119       ! Trends
120       CALL allocate_field(f_dps,field_t,type_real,name='dps')
121       CALL allocate_field(f_dmass,field_t,type_real,llm, name='dmass')
122       CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm,nqdyn,name='dtheta_rhodz')
123       CALL allocate_field(f_du,field_u,type_real,llm,name='du')
124       ! Model state at previous time step (RK/MLF)
125       CALL allocate_field(f_psm1,field_t,type_real,name='psm1')
126       CALL allocate_field(f_massm1,field_t,type_real,llm, name='massm1')
127       CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm,nqdyn,name='theta_rhodzm1')
128       CALL allocate_field(f_um1,field_u,type_real,llm,name='um1')
129    CASE(hevi)
130       ! Trends
131       CALL allocate_fields(nb_stage,f_dps_slow, field_t,type_real,name='dps_slow', ondevice=.TRUE.)
132       CALL allocate_fields(nb_stage,f_dmass_slow, field_t,type_real,llm, name='dmass_slow', ondevice=.TRUE.)
133       CALL allocate_fields(nb_stage,f_dtheta_rhodz_slow, field_t,type_real,llm,nqdyn,name='dtheta_rhodz_fast', ondevice=.TRUE.)
134       CALL allocate_fields(nb_stage,f_du_slow, field_u,type_real,llm,name='du_slow', ondevice=.TRUE.)
135       CALL allocate_fields(nb_stage,f_du_fast, field_u,type_real,llm,name='du_fast', ondevice=.TRUE.)
136       CALL allocate_fields(nb_stage,f_dW_slow, field_t,type_real,llm+1,name='dW_slow')
137       CALL allocate_fields(nb_stage,f_dW_fast, field_t,type_real,llm+1,name='dW_fast')
138       CALL allocate_fields(nb_stage,f_dPhi_slow, field_t,type_real,llm+1,name='dPhi_slow')
139       CALL allocate_fields(nb_stage,f_dPhi_fast, field_t,type_real,llm+1,name='dPhi_fast')
140       f_dps => f_dps_slow(:,1)
141       f_du => f_du_slow(:,1)
142       f_dtheta_rhodz => f_dtheta_rhodz_slow(:,1)
143    END SELECT
144
145    SELECT CASE(scheme)
146    CASE(mlf)
147       ! Model state 2 time steps ago (MLF)
148       CALL allocate_field(f_psm2,field_t,type_real)
149       CALL allocate_field(f_massm2,field_t,type_real,llm)
150       CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm,nqdyn)
151       CALL allocate_field(f_um2,field_u,type_real,llm)
152    END SELECT
153
154    CALL init_theta2theta_rhodz
155    CALL init_dissip
156    CALL init_sponge
157    CALL init_observable
158    CALL init_guided
159    CALL init_advect_tracer
160    CALL init_check_conserve
161
162    CALL etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_geopot,f_W, f_q)
163
164    CALL transfert_request(f_phis,req_i0) 
165    CALL transfert_request(f_phis,req_i1) 
166
167    CALL init_message(f_ps,req_i0,req_ps0)
168    CALL init_message(f_mass,req_i0,req_mass0)
169    CALL init_message(f_rhodz,req_i0,req_rhodz0)
170    CALL init_message(f_theta_rhodz,req_i0,req_theta_rhodz0)
171    CALL init_message(f_u,req_e0_vect,req_u0)
172    CALL init_message(f_q,req_i0,req_q0)
173    CALL init_message(f_geopot,req_i0,req_geopot0)
174    CALL init_message(f_W,req_i0,req_W0)
175
176  END SUBROUTINE init_timeloop
177 
178  SUBROUTINE timeloop
179    USE abort_mod
180    USE dissip_gcm_mod
181    USE sponge_mod
182    USE observable_mod
183    USE etat0_mod
184    USE guided_mod
185    USE caldyn_mod
186    USE advect_tracer_mod
187    USE diagflux_mod
188    USE physics_mod
189    USE mpipara
190    USE transfert_mod
191    USE check_conserve_mod
192    USE xios_mod
193    USE output_field_mod
194    USE write_etat0_mod
195    USE restart_mod
196    USE checksum_mod
197    USE explicit_scheme_mod
198    USE hevi_scheme_mod
199    REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), ps(:)
200
201    REAL(rstd) :: adv_over_out ! ratio itau_adv/itau_out
202    INTEGER :: ind, it,l
203    LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating mass fluxes in time
204    LOGICAL, PARAMETER :: check_rhodz=.FALSE.
205    INTEGER(kind=8) :: start_clock, stop_clock, rate_clock
206    INTEGER :: itau_sync   ! best iteration for synchronisation and ensure 1+1=2
207    INTEGER :: i
208   
209    LOGICAL,SAVE :: first_physic=.TRUE.
210    !$OMP THREADPRIVATE(first_physic)   
211
212    itau_sync=1
213    DO i=2,3*sync_it
214      IF (MOD(86400,INT(i*dt))==0 .AND. ABS((sync_it-itau_sync)*1./sync_it )/sync_it < (sync_it-itau_sync)*1./sync_it) itau_sync=i
215    ENDDO
216    CALL getin("itau_sync",itau_sync)
217    IF (is_master) PRINT*,"Synchronize frontier every itau_sync =",itau_sync
218     
219    CALL switch_omp_distrib_level
220    CALL caldyn_BC(f_phis, f_geopot, f_wflux) ! set constant values in first/last interfaces
221
222    !$OMP BARRIER
223    DO ind=1,ndomain
224       IF (.NOT. assigned_domain(ind)) CYCLE
225       CALL swap_dimensions(ind)
226       CALL swap_geometry(ind)
227       rhodz=f_rhodz(ind); mass=f_mass(ind); ps=f_ps(ind)
228       IF(caldyn_eta==eta_mass) THEN
229          CALL compute_rhodz(.TRUE., ps, rhodz, ondevice=.FALSE.) ! save rhodz for transport scheme before dynamics update ps
230       ELSE
231          DO l=ll_begin,ll_end
232             rhodz(:,l)=mass(:,l)
233          ENDDO
234       END IF
235    END DO
236    !$OMP BARRIER
237
238    fluxt_zero=.TRUE.
239
240    IF(positive_theta) CALL copy_theta_to_q(f_theta_rhodz,f_rhodz,f_q)
241    IF(diagflux_on) THEN
242       adv_over_out = itau_adv*(1./itau_out)
243    ELSE
244       adv_over_out = 0.
245    END IF
246
247    CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,f_q,itau0) 
248
249    Call trace_on
250
251    IF (xios_output) THEN ! we must call update_calendar before any XIOS output
252      IF (is_omp_master) CALL xios_update_calendar(1)
253    END IF
254   
255!    IF (is_write_start) CALL write_etat0(itau0,f_ps, f_phis,f_theta_rhodz,f_u,f_q)
256    IF (is_write_start) CALL write_etat0(itau0,f_ps, f_phis,f_theta_rhodz,f_u,f_q, f_geopot, f_W)
257   
258    CALL write_output_fields_basic(.TRUE., f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)
259
260    !$OMP MASTER
261    CALL SYSTEM_CLOCK(start_clock, rate_clock)
262    !$OMP END MASTER   
263    call update_device_field(f_ps)
264    call update_device_field(f_mass)
265    CALL update_device_field(f_theta_rhodz)
266    CALL update_device_field(f_u)
267    CALL update_device_field(f_q)
268    CALL update_device_field(f_geopot)
269    CALL update_device_field(f_wflux)
270    CALL update_device_field(f_rhodz)
271    call reset_profiling() 
272
273
274    DO it=itau0+1,itau0+itaumax
275
276       CALL print_iteration(it, itau0, itaumax, start_clock, rate_clock)
277
278       CALL enter_profile(id_timeloop)
279
280       IF (xios_output) THEN
281
282          IF(it>itau0+1 .AND. is_omp_master) CALL xios_update_calendar(it-itau0)
283       ELSE
284          CALL update_time_counter(dt*it)
285       END IF
286
287       IF (it==itau0+1 .OR. MOD(it-1,itau_sync)==0) THEN
288          CALL send_message(f_ps,req_ps0)
289          CALL wait_message(req_ps0)
290          CALL send_message(f_mass,req_mass0)
291          CALL wait_message(req_mass0)
292          CALL send_message(f_rhodz,req_rhodz0)
293          CALL wait_message(req_rhodz0)
294          CALL send_message(f_theta_rhodz,req_theta_rhodz0) 
295          CALL wait_message(req_theta_rhodz0)
296          CALL send_message(f_u,req_u0)
297          CALL wait_message(req_u0)
298          CALL send_message(f_q,req_q0) 
299          CALL wait_message(req_q0)
300          IF(.NOT. hydrostatic) THEN
301             CALL send_message(f_geopot,req_geopot0)
302             CALL wait_message(req_geopot0)
303             CALL send_message(f_W,req_W0)
304             CALL wait_message(req_W0)
305          END IF
306       ENDIF
307
308!       CALL checksum(f_ps)
309!       CALL checksum(f_theta_rhodz)
310!       CALL checksum(f_u)
311!       CALL checksum(f_q)
312!       CALL checksum(f_mass)
313!       CALL checksum(f_geopot)
314!       CALL checksum(f_rhodz)
315!       CALL checksum(f_wflux)
316
317       CALL guided(it*dt, f_ps, f_rhodz, f_theta_rhodz,f_u,f_q)
318       
319       CALL enter_profile(id_dyn)
320       SELECT CASE(scheme_family)
321       CASE(explicit)
322          CALL abort_acc("explicit_scheme")
323          CALL explicit_scheme(it, fluxt_zero)
324       CASE(hevi)
325          CALL HEVI_scheme(it, fluxt_zero)
326       END SELECT
327       CALL exit_profile(id_dyn)
328
329!       CALL checksum(f_ps)
330!       CALL checksum(f_theta_rhodz)
331!       CALL checksum(f_u)
332!       CALL checksum(f_q)
333!       CALL checksum(f_mass)
334!       CALL checksum(f_geopot)
335!       CALL checksum(f_rhodz)
336!       CALL checksum(f_wflux)
337
338       CALL enter_profile(id_dissip)
339       IF (MOD(it,itau_dissip)==0) THEN
340         
341          CALL update_border(it*dt, f_ps, f_mass, f_theta_rhodz,f_u,f_q) ! for limited area metric
342         
343          CALL enter_profile(id_diags)
344          CALL check_conserve_detailed(it, AAM_dyn, &
345               f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
346          CALL exit_profile(id_diags)
347
348          CALL dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_u, f_dtheta_rhodz,f_du)
349         
350          CALL euler_scheme(.FALSE.)  ! update only u, theta
351          IF (iflag_sponge > 0) THEN
352             CALL sponge(f_mass,f_u,f_du,f_theta_rhodz,f_dtheta_rhodz)
353             CALL euler_scheme(.FALSE.)  ! update only u, theta
354          END IF
355
356          CALL update_border(it*dt,f_ps, f_mass, f_theta_rhodz,f_u,f_q) ! for limited area metric
357
358          CALL enter_profile(id_diags)
359          CALL check_conserve_detailed(it, AAM_dissip, &
360               f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
361          CALL exit_profile(id_diags)
362       END IF
363       CALL exit_profile(id_dissip)
364!       CALL checksum(f_ps)
365!       CALL checksum(f_theta_rhodz)
366!       CALL checksum(f_u)
367!       CALL checksum(f_q)
368!       CALL checksum(f_hfluxt)
369!       CALL checksum(f_wfluxt)
370!       CALL checksum(f_u)
371!       CALL checksum(f_rhodz)
372
373       CALL enter_profile(id_adv)
374       IF(MOD(it,itau_adv)==0) THEN
375          CALL update_border(it*dt, f_ps, f_mass, f_theta_rhodz,f_u,f_q) ! for limited area metric
376          CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz, & ! update q and rhodz after RK step
377               adv_over_out, f_masst,f_qmasst,f_massfluxt, f_qfluxt)  ! accumulate mass and fluxes if diagflux_on
378          fluxt_zero=.TRUE. ! restart accumulation of hfluxt and qfluxt at next call to explicit_scheme / HEVI_scheme
379          ! At this point advect_tracer has obtained the halos of u and rhodz,
380          ! needed for correct computation of kinetic energy
381          IF(diagflux_on) CALL abort_acc("diagflux_on")
382          IF(diagflux_on) CALL diagflux_energy(adv_over_out, f_phis,f_rhodz,f_theta_rhodz,f_u, f_geopot,f_theta,f_buf_i, f_hfluxt)
383
384          DO ind=1,ndomain
385            IF (.NOT. assigned_domain(ind)) CYCLE
386            CALL swap_dimensions(ind)
387            CALL swap_geometry(ind)
388            rhodz=f_rhodz(ind); mass=f_mass(ind); ps=f_ps(ind)
389            IF(caldyn_eta==eta_mass) THEN
390              CALL compute_rhodz(.TRUE., ps, rhodz, ondevice=.TRUE.) ! save rhodz for transport scheme before dynamics update ps
391            ELSE
392              DO l=ll_begin,ll_end
393                 rhodz(:,l)=mass(:,l)
394              ENDDO
395            END IF
396          END DO
397
398          IF(positive_theta) CALL abort_acc("positive_theta")
399          IF(positive_theta) CALL copy_q_to_theta(f_theta_rhodz,f_rhodz,f_q)
400          CALL update_border(it*dt, f_ps, f_mass, f_theta_rhodz,f_u,f_q) ! for limited area metric
401       END IF
402       CALL exit_profile(id_adv)
403!       CALL checksum(f_ps)
404!       CALL checksum(f_theta_rhodz)
405!       CALL checksum(f_u)
406!       CALL checksum(f_q)
407
408       CALL enter_profile(id_diags)
409!       IF (MOD(it,itau_physics)==0) THEN
410          CALL check_conserve_detailed(it, AAM_dyn, &
411            f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
412          CALL enter_profile(id_phys)
413          CALL physics(it,f_phis, f_geopot, f_ps, f_theta_rhodz, f_u, f_wflux, f_q)
414          CALL exit_profile(id_phys)
415          CALL check_conserve_detailed(it, AAM_phys, &
416               f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
417          !$OMP MASTER
418          IF (first_physic) CALL SYSTEM_CLOCK(start_clock)
419          !$OMP END MASTER   
420          first_physic=.FALSE.
421!       END IF
422
423       IF (MOD(it,itau_check_conserv)==0) THEN
424          CALL update_host_field(f_ps)
425          CALL update_host_field(f_theta_rhodz)
426          CALL update_host_field(f_u)
427          CALL update_host_field(f_dps)
428          CALL update_host_field(f_q)
429          CALL check_conserve_detailed(it, AAM_dyn, &
430               f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
431          CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,f_q,it) 
432       ENDIF
433
434       IF (mod(it,itau_out)==0 ) THEN
435 
436          CALL update_border(it*dt,f_ps, f_mass, f_theta_rhodz,f_u,f_q) ! for limited area metric
437         
438          CALL transfert_request(f_u,req_e1_vect)
439          CALL update_host_field(f_ps)             
440          CALL update_host_field(f_mass)
441          CALL update_host_field(f_theta_rhodz)
442          CALL update_host_field(f_geopot)
443          CALL update_host_field(f_u)
444          CALL update_host_field(f_q)
445          CALL write_output_fields_basic(.FALSE.,f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)
446       ENDIF
447       CALL exit_profile(id_diags)
448
449       CALL exit_profile(id_timeloop)
450    END DO
451   
452    CALL update_host_field(f_ps)
453    CALL update_host_field(f_theta_rhodz)
454    CALL update_host_field(f_u)
455    CALL update_host_field(f_q)
456    CALL update_host_field(f_geopot)
457
458!    CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q)
459    CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q, f_geopot, f_W) 
460
461    CALL update_host_field(f_dps)   
462    CALL check_conserve_detailed(it, AAM_dyn, &
463         f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
464    CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,f_q,it) 
465   
466    !$OMP MASTER
467    CALL SYSTEM_CLOCK(stop_clock)
468   
469    IF (mpi_rank==0) THEN
470       PRINT *,"Time elapsed : ",(stop_clock-start_clock)*1./rate_clock 
471    ENDIF
472    !$OMP END MASTER
473   
474    ! CONTAINS
475  END SUBROUTINE timeloop
476
477  SUBROUTINE print_iteration(it,itau0,itaumax,start_clock,rate_clock)
478    INTEGER :: it, itau0, itaumax, throughput
479    INTEGER(kind=8) :: start_clock, stop_clock, rate_clock
480    REAL :: per_step,total, elapsed
481   
482    IF(is_master) THEN
483       WRITE(*,'(A,I7,A,F14.1)') "It No :",it,"   t :",dt*it
484       IF(MOD(it,10)==0) THEN
485          CALL SYSTEM_CLOCK(stop_clock)
486          elapsed = (stop_clock-start_clock)*1./rate_clock
487          per_step = elapsed/(it-itau0)
488          throughput = INT(dt/per_step)
489          total = per_step*itaumax
490          WRITE(*,'(A,I5,A,F6.2,A,I6)') 'Time spent (s):',INT(elapsed), &
491               '  -- ms/step : ', 1000*per_step, &
492               '  -- Throughput :', throughput
493          WRITE(*,'(A,I5,A,I5)') 'Whole job (min) :', INT(total/60.), &
494               '  -- Completion in (min) : ', INT((total-elapsed)/60.)
495       END IF
496    END IF
497    IF(MOD(it,itau_prof)==0) CALL print_profile
498
499  END SUBROUTINE print_iteration
500
501  SUBROUTINE copy_theta_to_q(f_theta_rhodz,f_rhodz,f_q)
502    TYPE(t_field),POINTER :: f_theta_rhodz(:),f_rhodz(:), f_q(:)
503    REAL(rstd), POINTER :: theta_rhodz(:,:,:), rhodz(:,:), q(:,:,:)
504    INTEGER :: ind, iq
505    DO ind=1, ndomain
506       IF (.NOT. assigned_domain(ind)) CYCLE
507       CALL swap_dimensions(ind)
508       CALL swap_geometry(ind)
509       theta_rhodz=f_theta_rhodz(ind)
510       rhodz=f_rhodz(ind)
511       q=f_q(ind)
512       DO iq=1, nqdyn
513          q(:,:,iq)  = theta_rhodz(:,:,iq)/rhodz(:,:)
514       END DO
515    END DO
516  END SUBROUTINE copy_theta_to_q
517
518  SUBROUTINE copy_q_to_theta(f_theta_rhodz,f_rhodz,f_q)
519    TYPE(t_field),POINTER :: f_theta_rhodz(:),f_rhodz(:), f_q(:)
520    REAL(rstd), POINTER :: theta_rhodz(:,:,:), rhodz(:,:), q(:,:,:)
521    INTEGER :: ind, iq
522    DO ind=1, ndomain
523       IF (.NOT. assigned_domain(ind)) CYCLE
524       CALL swap_dimensions(ind)
525       CALL swap_geometry(ind)
526       theta_rhodz=f_theta_rhodz(ind)
527       rhodz=f_rhodz(ind)
528       q=f_q(ind)
529       DO iq=1,nqdyn
530          theta_rhodz(:,:,iq) = rhodz(:,:)*q(:,:,iq)
531       END DO
532    END DO
533  END SUBROUTINE copy_q_to_theta
534
535END MODULE timeloop_gcm_mod
Note: See TracBrowser for help on using the repository browser.