source: trunk/Roms_agrif/step.F @ 4

Last change on this file since 4 was 4, checked in by pinsard, 17 years ago

modification according to /usr/temp/vech/quiberon/Roms_tools/Roms_Agrif_pisces/

File size: 14.4 KB
Line 
1!
2! $Id: step.F,v 1.18 2005/10/26 15:18:52 pmarches Exp $
3!
4#include "cppdefs.h"
5      subroutine step()
6      implicit none
7#include "param.h"
8#include "scalars.h"
9#include "zoom.h"
10#include "grid.h"
11#include "coupling.h"
12#include "ocean3d.h"
13#include "ocean2d.h"
14
15!$AGRIF_DO_NOT_TREAT
16      integer ntrds,trd,
17     &        my_first,my_last, tile, my_iif
18      common /tile_ranges/ ntrds, trd, my_first,my_last
19!$AGRIF_END_DO_NOT_TREAT
20C$OMP THREADPRIVATE (/tile_ranges/)
21
22      integer :: iifparent
23
24#ifdef AGRIF       
25       IF (.Not.Agrif_Root()) then
26       iifparent =
27     &Agrif_Parent(iif)
28
29         IF (iifparent == (nfast+1)) THEN
30           IF (iif == nfast) THEN
31             iif = mod(iif+1, nfast+2)
32C$OMP PARALLEL
33             call step3D_thread()
34C$OMP END PARALLEL           
35             nbstep3d = nbstep3d + 1
36           ENDIF
37           Return
38         ENDIF
39             
40         IF (iifparent == 0) THEN
41           IF (iif == (nfast+1)) THEN
42             iif = mod(iif+1, nfast+2)
43C$OMP PARALLEL
44             call prestep3D_thread()
45C$OMP END PARALLEL
46           ENDIF
47           Return
48         ENDIF
49                 
50       ENDIF
51#endif
52
53!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
54#ifdef SOLVE3D       
55       iif = mod(iif+1, nfast+2)
56      IF (iif .EQ. 0) then
57C$OMP PARALLEL
58        call prestep3D_thread()
59C$OMP END PARALLEL     
60      endif
61       
62      IF ((iif.NE.0).AND.(iif.NE.(nfast+1))) THEN
63#endif
64C$OMP PARALLEL
65        call step2d_thread()
66C$OMP END PARALLEL
67#ifdef SOLVE3D
68      ENDIF
69
70      IF (iif .EQ. (nfast+1)) then
71C$OMP PARALLEL
72        call step3D_thread()
73C$OMP END PARALLEL     
74        nbstep3d = nbstep3d + 1
75      endif
76#endif     
77!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
78
79#ifdef AGRIF     
80      IF (.Not.Agrif_Root()) THEN             
81        IF (iif == (nfast+1)) THEN
82          iif = 1 + mod(iif, nfast+1)
83C$OMP PARALLEL
84          call prestep3D_thread()
85C$OMP END PARALLEL
86C$OMP PARALLEL
87          call step2d_thread()
88C$OMP END PARALLEL
89        ENDIF
90      ENDIF
91#endif
92
93      return
94      end
95
96      subroutine prestep3D_thread()
97#ifdef AGRIF
98      use Agrif_Util
99#endif
100      implicit none
101#include "param.h"
102#include "scalars.h"
103#include "ncscrum.h"
104#ifdef FLOATS
105# include "floats.h"
106#endif
107#include "ocean2d.h"
108#ifdef STATIONS
109# include "sta.h"
110#endif
111       integer range
112!$AGRIF_DO_NOT_TREAT
113      integer ntrds,trd,
114     &        my_first,my_last, tile, my_iif
115      common /tile_ranges/ ntrds, trd, my_first,my_last
116!$AGRIF_END_DO_NOT_TREAT
117C$OMP THREADPRIVATE (/tile_ranges/)
118#ifdef FLOATS
119       integer chunk_size_flt, Lstr,Lend, flt_str
120       common /floats_step/ flt_str
121#endif
122#ifdef STATIONS
123       integer chunk_size_sta, Mcpstr,Mcpend, sta_str
124       common /sta_step/ sta_str
125#endif
126#ifdef AGRIF
127      integer Agrif_lev_sedim
128#endif
129      integer omp_get_num_threads, omp_get_thread_num 
130
131#ifdef AGRIF
132      Agrif_lev_sedim=Agrif_Nb_Fine_Grids()
133#endif
134
135      ntrds=omp_get_num_threads()
136      trd=omp_get_thread_num()
137C$OMP BARRIER
138      range=(NSUB_X*NSUB_E+ntrds-1)/ntrds
139      my_first=trd*range
140      my_last=min(my_first + range-1, NSUB_X*NSUB_E-1)
141
142C$OMP MASTER
143
144       
145      time=time_start+dt*float(iic-ntstart)
146      tdays=time*sec2day
147
148#ifdef SOLVE3D
149      nstp=1+mod(iic-ntstart,2)
150      nrhs=nstp
151      nnew=3
152#endif
153
154       
155#ifdef FLOATS
156      nfp1=MOD(nfp1+1,NFT+1)  ! Shift time indices
157      nf  =MOD(nf  +1,NFT+1)  ! for floats
158      nfm1=MOD(nfm1+1,NFT+1)
159      nfm2=MOD(nfm2+1,NFT+1)
160      nfm3=MOD(nfm3+1,NFT+1)
161      flt_str=0
162#endif
163#ifdef STATIONS
164      sta_str=0
165#endif
166#if defined AGRIF
167      Agrif_lev_sedim=Agrif_Nb_Fine_Grids()
168#endif
169
170!
171! Model input block: read forcing/climatology data from netCDF files.
172!------ ----- ------ ---- ------------------- ---- ---- ------ ------
173!
174      if (synchro_flag) then
175#ifdef AGRIF
176        IF (Agrif_Root()) THEN
177#endif
178#if defined SOLVE3D && defined TCLIMATOLOGY && !defined ANA_TCLIMA
179            call get_tclima
180#endif
181#if defined SOLVE3D && defined M3CLIMATOLOGY && !defined ANA_M3CLIMA \
182 || (defined M2CLIMATOLOGY && !defined ANA_M2CLIMA)
183            call get_uclima
184#endif
185#if defined ZCLIMATOLOGY && !defined ANA_SSH
186            call get_ssh
187#endif
188#ifdef FRC_BRY
189# ifndef ANA_BRY
190            call get_bry
191#  ifdef BIOLOGY
192            call get_bry_bio
193#  endif
194# endif
195#endif
196#ifdef AGRIF
197        ENDIF
198#endif
199        call get_vbc
200#if defined BBL && !defined ANA_WWAVE
201# ifdef AGRIF
202        IF (Agrif_Fixed().EQ.Agrif_lev_sedim) THEN
203# endif
204          call get_wwave
205# ifdef AGRIF
206        ENDIF
207# endif
208#endif /* BBL && !ANA_WWAVE*/
209        synchro_flag=.false.
210      endif  !<-- synchro_flag
211C$OMP END MASTER
212C$OMP BARRIER
213
214      if (may_day_flag.ne.0) return  !-->  EXIT
215
216#ifdef SOLVE3D
217
218      do tile=my_first,my_last
219# ifdef AGRIF
220        IF (Agrif_Root()) THEN
221# endif
222# ifdef TCLIMATOLOGY
223          call ana_tclima (tile)   ! analytical values are given
224          call set_tclima (tile)   ! if data is missing from clim files
225# endif
226# if defined M2CLIMATOLOGY || defined M3CLIMATOLOGY
227#  if defined ANA_M2CLIMA || defined ANA_M3CLIMA
228          call ana_uclima (tile)
229#  else
230          call set_uclima (tile)
231#  endif
232# endif
233# ifdef ZCLIMATOLOGY
234#  ifdef ANA_SSH
235          call ana_ssh (tile)
236#  else
237          call set_ssh (tile)
238#  endif
239# endif
240# ifdef FRC_BRY
241#  ifdef ANA_BRY
242          call ana_bry(tile)
243#  else
244          call set_bry(tile)
245#  endif
246#  ifdef BIOLOGY
247          call ana_bry_bio (tile)  ! analytical values are given
248          call set_bry_bio (tile)  ! if data is missing from bry files
249#  endif
250# endif
251# ifdef AGRIF
252        ENDIF
253# endif
254        call rho_eos (tile)
255        call set_HUV (tile)
256# if defined BIOLOGY && !defined MPI && !defined PISCES
257        call bio_diag (tile)
258# else
259        call diag(tile)
260# endif
261      enddo
262C$OMP BARRIER
263      do tile=my_first,my_last
264        call set_vbc (tile)
265# ifdef BBL
266#  ifdef AGRIF
267        IF (Agrif_Fixed().EQ.Agrif_lev_sedim) THEN
268#  endif
269#  ifdef ANA_WWAVE
270          call ana_wwave (tile)
271#  else
272          call set_wwave (tile)
273#  endif
274          call bblm (tile)
275#  ifdef AGRIF
276        ENDIF
277#  endif
278# endif
279# if defined SSH_TIDES || defined UV_TIDES
280#  ifdef AGRIF
281        IF (Agrif_Root()) call clm_tides (tile)
282#  else
283        call clm_tides (tile)
284#  endif
285# endif
286      enddo
287C$OMP BARRIER
288
289      if (may_day_flag.ne.0) return  !-->  EXIT
290
291      do tile=my_first,my_last
292# if defined ANA_VMIX
293        call ana_vmix (tile)
294# elif defined LMD_MIXING
295        call lmd_vmix (tile)
296# elif defined BVF_MIXING
297        call bvf_mix (tile)
298# endif
299        call omega (tile)
300# if (defined UV_VIS2 || defined TS_DIF2) && defined SMAGORINSKY
301        call smagorinsky (tile)
302# endif
303      enddo
304C$OMP BARRIER
305      do tile=my_first,my_last
306        call prsgrd (tile)
307        call rhs3d (tile)
308        call pre_step3d (tile)
309# if defined UV_VIS2 || defined UV_VIS4
310# if defined AGRIF && !defined SMAGO_UV
311        if (Agrif_Root()) then
312# endif
313        call uv3dmix (tile) 
314# if defined AGRIF && !defined SMAGO_UV
315        else
316        call uv3dmix_fine(tile)
317        endif
318# endif
319# endif
320# ifdef AVERAGES
321        call set_avg (tile)
322#  if defined DIAGNOSTICS_TS
323        call set_diags_avg (tile)
324#  endif
325#  if defined DIAGNOSTICS_UV
326        call set_diagsM_avg (tile)
327#  endif
328#  ifdef DIAGNOSTICS_BIO
329        call set_bio_diags_avg (tile)
330#  endif
331# endif
332      enddo
333
334C$OMP BARRIER
335C$OMP MASTER
336      nrhs=3
337      nnew=3-nstp
338!
339! Output block: write restart/history files.   
340!
341      call output
342
343C$OMP END MASTER
344C$OMP BARRIER
345#endif  /* SOLVE3D */
346
347      if (may_day_flag.ne.0) return  !-->  EXIT
348
349       
350#if defined AGRIF && defined AGRIF_2WAY
351!
352! Update the outer domain after the last child step
353! in case of 2-way nesting.
354!
355C$OMP BARRIER
356C$OMP SINGLE
357      if (.Not.Agrif_Root()) then
358        call Agrif_update_np1pre
359      endif
360C$OMP END SINGLE
361#endif /*AGRIF && AGRIF_2WAY*/
362
363      return
364      end
365     
366      subroutine step2D_thread()
367#ifdef AGRIF
368      use Agrif_Util
369#endif
370      implicit none
371#include "param.h"
372#include "scalars.h"
373#include "ncscrum.h"
374#include "ocean2d.h"
375#ifdef FLOATS
376# include "floats.h"
377#endif
378#ifdef STATIONS
379# include "sta.h"
380#endif
381       integer range
382!$AGRIF_DO_NOT_TREAT
383      integer ntrds,trd,
384     &        my_first,my_last, tile, my_iif
385      common /tile_ranges/ ntrds, trd, my_first,my_last
386!$AGRIF_END_DO_NOT_TREAT
387C$OMP THREADPRIVATE (/tile_ranges/)
388#ifdef FLOATS
389       integer chunk_size_flt, Lstr,Lend, flt_str
390       common /floats_step/ flt_str
391#endif
392#ifdef STATIONS
393       integer chunk_size_sta, Mcpstr,Mcpend, sta_str
394       common /sta_step/ sta_str
395#endif
396#ifdef AGRIF
397      integer Agrif_lev_sedim
398#endif
399      integer omp_get_num_threads, omp_get_thread_num 
400
401#ifdef AGRIF
402      Agrif_lev_sedim=Agrif_Nb_Fine_Grids()
403#endif
404
405!
406! Solve the 2D primitive equations for the barotropic mode.
407!-------------------------------------------------------------
408! Note that no actual fast-time-step is performed during the
409! auxiliary (nfast+1)th step. It is needed to finalize the fast-time
410! averaging procedure, if any, and to compute the total depth of
411! water column, as well as the new vertical coordinate transformation
412! profiles, z_r, z_w, because the free surface elevation has been
413! changed. All these operations are done during predictor time step.
414! Thus, there is no need for the corrector step during the auxiliary
415! time step.
416!
417
418C$OMP MASTER
419        kstp=knew                   ! This might look a bit silly,
420        knew=kstp+1                 ! since both branches of this
421        if (knew.gt.4) knew=1       ! if statement are identical.
422C$OMP END MASTER
423C$OMP BARRIER
424!
425! Model input block for 2D configurations only !!!
426!
427#ifndef SOLVE3D
428# ifdef AGRIF
429        IF (Agrif_Root()) THEN
430# endif
431          do tile=my_first,my_last
432# if defined M2CLIMATOLOGY
433#  if defined ANA_M2CLIMA
434            call ana_uclima (tile)
435#  else
436            call set_uclima (tile)
437#  endif
438# endif
439# ifdef ZCLIMATOLOGY
440#  ifdef ANA_SSH
441            call ana_ssh (tile)
442#  else
443            call set_ssh (tile)
444#  endif
445# endif
446# if  defined Z_FRC_BRY  || defined M2_FRC_BRY
447#  ifdef ANA_BRY
448            call ana_bry(tile)
449#  else
450            call set_bry(tile)
451#  endif
452# endif
453          enddo
454# ifdef AGRIF
455        ENDIF
456# endif
457C$OMP BARRIER
458        do tile=my_first,my_last
459          call set_vbc (tile)
460# if defined SSH_TIDES || defined UV_TIDES
461#   ifdef AGRIF
462          IF (Agrif_Root()) call clm_tides (tile)
463#   else
464          call clm_tides (tile)
465#   endif
466# endif
467# ifdef AVERAGES
468          call set_avg (tile)
469# endif
470        enddo
471C$OMP BARRIER
472C$OMP MASTER
473        call output
474C$OMP END MASTER
475C$OMP BARRIER
476#endif  /* !defined SOLVE3D */
477       
478                  ! This might look a bit silly,
479                  ! since both branches of this
480                  ! if statement are identical.
481                  ! Nevertheless, it makes sense,
482                  ! since mpc will reverse one of
483                  ! these loops to make zig-zag
484                  ! tile-processing sequence.               
485
486        if (mod(knew,2).eq.0) then
487          do tile=my_first,my_last
488            call     step2d (tile)
489          enddo
490        else
491          do tile=my_first,my_last
492            call     step2d (tile)
493          enddo
494        endif
495
496#if defined AGRIF && defined AGRIF_2WAY
497        if (.not.agrif_root()) then
498            call update2d()
499        endif
500#endif
501
502      return
503      end
504     
505      subroutine step3D_thread()
506#ifdef AGRIF
507      use Agrif_Util
508#endif
509      implicit none
510#include "param.h"
511#include "scalars.h"
512#include "ncscrum.h"
513#ifdef FLOATS
514# include "floats.h"
515#endif
516#ifdef STATIONS
517# include "sta.h"
518#endif
519#include "zoom.h"
520#include "ocean3d.h"
521#include "ocean2d.h"
522#include "coupling.h"
523       integer range
524!$AGRIF_DO_NOT_TREAT
525      integer ntrds,trd,
526     &        my_first,my_last, tile, my_iif
527      common /tile_ranges/ ntrds, trd, my_first,my_last
528      integer i,j,k
529!$AGRIF_END_DO_NOT_TREAT
530C$OMP THREADPRIVATE (/tile_ranges/)
531#ifdef FLOATS
532       integer chunk_size_flt, Lstr,Lend, flt_str
533       common /floats_step/ flt_str
534#endif
535#ifdef STATIONS
536       integer chunk_size_sta, Mcpstr,Mcpend, sta_str
537       common /sta_step/ sta_str
538#endif
539#ifdef AGRIF
540      integer Agrif_lev_sedim
541#endif
542      integer omp_get_num_threads, omp_get_thread_num 
543     
544     
545      real t1,t2,t3
546
547#ifdef AGRIF
548      Agrif_lev_sedim=Agrif_Nb_Fine_Grids()
549#endif
550   
551C$OMP BARRIER
552
553#ifdef SOLVE3D
554
555      do tile=my_first,my_last
556        call set_depth (tile)
557      enddo
558C$OMP BARRIER
559
560      do tile=my_first,my_last
561        call set_HUV2 (tile)
562      enddo
563C$OMP BARRIER
564
565      do tile=my_first,my_last
566        call omega (tile)
567        call rho_eos (tile)
568      enddo
569C$OMP BARRIER
570
571      do tile=my_first,my_last
572        call prsgrd (tile)
573        call rhs3d (tile)
574        call step3d_uv1 (tile)
575      enddo
576C$OMP BARRIER
577
578      do tile=my_first,my_last
579        call step3d_uv2 (tile)
580      enddo
581C$OMP BARRIER
582 
583      do tile=my_first,my_last
584        call omega (tile)
585        call step3d_t (tile)
586# ifdef SEDIMENT
587#  ifdef AGRIF
588        IF (Agrif_Fixed().EQ.Agrif_lev_sedim)
589     &    call sediment (tile)
590#  else
591        call sediment (tile)
592#  endif
593# endif
594
595# if defined TS_DIF2 || defined TS_DIF4
596# if defined AGRIF && !defined SMAGO_TS
597        if (Agrif_Root()) then
598# endif
599        call t3dmix (tile) 
600# if defined AGRIF && !defined SMAGO_TS
601        else
602        call t3dmix_fine(tile)
603        endif
604# endif
605# endif
606      enddo
607C$OMP BARRIER
608
609#endif /* SOLVE3D */
610
611#ifdef FLOATS
612      chunk_size_flt=32
613      do while (flt_str.lt.nfloats) 
614C$OMP CRITICAL
615        Lstr=flt_str+1
616        flt_str=Lstr+chunk_size_flt-1
617C$OMP END CRITICAL
618        Lend=min(Lstr+chunk_size_flt-1,nfloats)
619        call step_floats (Lstr,Lend)
620      enddo
621c        call step_floats (1,nfloats) ! serial version for debugging
622#endif /*FLOATS*/
623
624#ifdef STATIONS
625      chunk_size_sta=32
626      do while (sta_str.lt.nstas)
627C$OMP CRITICAL
628        Mcpstr=sta_str+1
629        sta_str=Mcpstr+chunk_size_sta-1
630C$OMP END CRITICAL
631        Mcpend=min(Mcpstr+chunk_size_sta-1,nstas)
632        call step_sta (Mcpstr,Mcpend)
633      enddo
634#endif /*STATIONS*/
635!
636#if defined AGRIF && defined AGRIF_2WAY
637!
638! Update the outer domain after the last child step
639! in case of 2-way nesting.
640!
641C$OMP BARRIER 
642C$OMP MASTER
643      if ((.Not.Agrif_Root()).and.
644     &    (nbcoarse == Agrif_Irhot())) then
645        call Agrif_update_np1
646      endif
647C$OMP END MASTER
648#endif /*AGRIF && AGRIF_2WAY*/
649
650C$OMP MASTER
651      iic=iic + 1
652
653#ifdef AGRIF     
654      nbcoarse = 1 + mod(nbcoarse, Agrif_IRhot())
655#endif     
656C$OMP END MASTER
657C$OMP BARRIER
658
659      return
660      end           
661
Note: See TracBrowser for help on using the repository browser.