source: LMDZ6/trunk/libf/dyn3d/dynetat0.F90 @ 5182

Last change on this file since 5182 was 5084, checked in by Laurent Fairhead, 12 months ago

Reverting to r4065. Updating fortran standard broke too much stuff. Will do it by smaller chunks
AB, LF
File size: 10.1 KB
RevLine 
[2293]1SUBROUTINE dynetat0(fichnom,vcov,ucov,teta,q,masse,ps,phis,time)
2!
3!-------------------------------------------------------------------------------
4! Authors: P. Le Van , L.Fairhead
5!-------------------------------------------------------------------------------
6! Purpose: Initial state reading.
7!-------------------------------------------------------------------------------
[4325]8  USE infotrac,    ONLY: nqtot, tracers, niso, iqIsoPha, iH2O, isoName
[4193]9  USE strings_mod, ONLY: maxlen, msg, strStack, real2str, int2str
[5084]10  USE netcdf,      ONLY: NF90_OPEN,  NF90_NOWRITE, NF90_INQ_VARID, &
[4120]11                         NF90_CLOSE, NF90_GET_VAR, NF90_NoErr
[4325]12  USE readTracFiles_mod, ONLY: new2oldH2O, newHNO3, oldHNO3, getKey
[2293]13  USE control_mod, ONLY: planet_type
14  USE assert_eq_m, ONLY: assert_eq
[2600]15  USE comvert_mod, ONLY: pa,preff
[2597]16  USE comconst_mod, ONLY: cpp, daysec, dtvr, g, im, jm, kappa, lllm, omeg, rad
[2603]17  USE logic_mod, ONLY: fxyhypb, ysinus
[2598]18  USE serre_mod, ONLY: clon, clat, grossismx, grossismy
[2601]19  USE temps_mod, ONLY: annee_ref, day_ini, day_ref, itau_dyn, start_time
[2622]20  USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
[2600]21
[2293]22  IMPLICIT NONE
23  include "dimensions.h"
24  include "paramet.h"
25  include "comgeom2.h"
26  include "description.h"
27  include "iniprint.h"
28!===============================================================================
29! Arguments:
30  CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
31  REAL, INTENT(OUT) ::  vcov(iip1,jjm, llm)        !--- V COVARIANT WIND
32  REAL, INTENT(OUT) ::  ucov(iip1,jjp1,llm)        !--- U COVARIANT WIND
33  REAL, INTENT(OUT) ::  teta(iip1,jjp1,llm)        !--- POTENTIAL TEMP.
34  REAL, INTENT(OUT) ::     q(iip1,jjp1,llm,nqtot)  !--- TRACERS
35  REAL, INTENT(OUT) :: masse(iip1,jjp1,llm)        !--- MASS PER CELL
36  REAL, INTENT(OUT) ::    ps(iip1,jjp1)            !--- GROUND PRESSURE
37  REAL, INTENT(OUT) ::  phis(iip1,jjp1)            !--- GEOPOTENTIAL
38!===============================================================================
39! Local variables:
[4120]40  CHARACTER(LEN=maxlen) :: mesg, var, modname, oldVar
[2293]41  INTEGER, PARAMETER :: length=100
[4050]42  INTEGER :: iq, fID, vID, idecal, iqParent, iName, iZone, iPhase
[4325]43  REAL    :: time, tnat, alpha_ideal, tab_cntrl(length)    !--- RUN PARAMS TABLE
[4301]44  LOGICAL :: lSkip, ll
[4984]45  LOGICAL,PARAMETER :: tnat1=.TRUE.
[2293]46!-------------------------------------------------------------------------------
47  modname="dynetat0"
48
49!--- Initial state file opening
50  var=fichnom
51  CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
52  CALL get_var1("controle",tab_cntrl)
53
54!!! AS: idecal is a hack to be able to read planeto starts...
55!!!     .... while keeping everything OK for LMDZ EARTH
56  IF(planet_type=="generic") THEN
[4120]57    CALL msg('NOTE NOTE NOTE : Planeto-like start files', modname)
[2293]58    idecal = 4
59    annee_ref  = 2000
60  ELSE
[4120]61    CALL msg('NOTE NOTE NOTE : Earth-like start files', modname)
[2293]62    idecal = 5
63    annee_ref  = tab_cntrl(5)
64  END IF
65  im         = tab_cntrl(1)
66  jm         = tab_cntrl(2)
67  lllm       = tab_cntrl(3)
68  day_ref    = tab_cntrl(4)
69  rad        = tab_cntrl(idecal+1)
70  omeg       = tab_cntrl(idecal+2)
71  g          = tab_cntrl(idecal+3)
72  cpp        = tab_cntrl(idecal+4)
73  kappa      = tab_cntrl(idecal+5)
74  daysec     = tab_cntrl(idecal+6)
75  dtvr       = tab_cntrl(idecal+7)
76  etot0      = tab_cntrl(idecal+8)
77  ptot0      = tab_cntrl(idecal+9)
78  ztot0      = tab_cntrl(idecal+10)
79  stot0      = tab_cntrl(idecal+11)
80  ang0       = tab_cntrl(idecal+12)
81  pa         = tab_cntrl(idecal+13)
82  preff      = tab_cntrl(idecal+14)
83!
84  clon       = tab_cntrl(idecal+15)
85  clat       = tab_cntrl(idecal+16)
86  grossismx  = tab_cntrl(idecal+17)
87  grossismy  = tab_cntrl(idecal+18)
88!
89  IF ( tab_cntrl(idecal+19)==1. )  THEN
90    fxyhypb  = .TRUE.
91!   dzoomx   = tab_cntrl(25)
92!   dzoomy   = tab_cntrl(26)
93!   taux     = tab_cntrl(28)
94!   tauy     = tab_cntrl(29)
95  ELSE
96    fxyhypb = .FALSE.
97    ysinus  = tab_cntrl(idecal+22)==1.
98  END IF
99
100  day_ini    = tab_cntrl(30)
101  itau_dyn   = tab_cntrl(31)
102  start_time = tab_cntrl(32)
103
104!-------------------------------------------------------------------------------
[4120]105  CALL msg('rad, omeg, g, cpp, kappa = '//TRIM(strStack(real2str([rad,omeg,g,cpp,kappa]))), modname)
[2293]106  CALL check_dim(im,iim,'im','im')
107  CALL check_dim(jm,jjm,'jm','jm')
108  CALL check_dim(lllm,llm,'lm','lllm')
109  CALL get_var1("rlonu",rlonu)
110  CALL get_var1("rlatu",rlatu)
111  CALL get_var1("rlonv",rlonv)
112  CALL get_var1("rlatv",rlatv)
113  CALL get_var2("cu"   ,cu)
114  CALL get_var2("cv"   ,cv)
115  CALL get_var2("aire" ,aire)
116  var="temps"
117  IF(NF90_INQ_VARID(fID,var,vID)/=NF90_NoErr) THEN
[4120]118    CALL msg('missing field <temps> ; trying with <Time>', modname)
119    var="Time"
[2293]120    CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
121  END IF
122  CALL err(NF90_GET_VAR(fID,vID,time),"get",var)
[2299]123  CALL get_var2("phisinit",phis)
[2293]124  CALL get_var3("ucov",ucov)
125  CALL get_var3("vcov",vcov)
126  CALL get_var3("teta",teta)
127  CALL get_var3("masse",masse)
128  CALL get_var2("ps",ps)
129
130!--- Tracers
[4268]131  ll=.FALSE.
[4265]132#ifdef REPROBUS
[4301]133  ll = NF90_INQ_VARID(fID, 'HNO3tot', vID) /= NF90_NoErr                                 !--- DETECT OLD REPRO start.nc FILE
[4265]134#endif
[2293]135  DO iq=1,nqtot
[4120]136    var = tracers(iq)%name
[4301]137    oldVar = new2oldH2O(var)
138    lSkip = ll .AND. var == 'HNO3'                                                       !--- FORCE "HNO3_g" READING FOR "HNO3"
139#ifdef REPROBUS
140    ix = strIdx(newHNO3, var); IF(ix /= 0) oldVar = oldHNO3(ix)                          !--- REPROBUS HNO3 exceptions
141#endif
142#ifdef INCA
143    IF(var == 'O3') oldVar = 'OX'                                                        !--- DEAL WITH INCA OZONE EXCEPTION
144#endif
[4120]145    !--------------------------------------------------------------------------------------------------------------------------
[4301]146    IF(NF90_INQ_VARID(fID, var, vID) == NF90_NoErr .AND. .NOT.lSkip) THEN                !=== REGULAR CASE: AVAILABLE VARIABLE
[4120]147      CALL err(NF90_GET_VAR(fID,vID,q(:,:,:,iq)),"get",var)
148    !--------------------------------------------------------------------------------------------------------------------------
[4301]149    ELSE IF(NF90_INQ_VARID(fID, oldVar, vID) == NF90_NoErr) THEN                         !=== TRY WITH ALTERNATE NAME
[4120]150      CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized to <'//TRIM(oldVar)//'>', modname)
151      CALL err(NF90_GET_VAR(fID,vID,q(:,:,:,iq)),"get",oldVar)
152    !--------------------------------------------------------------------------------------------------------------------------
153    ELSE IF(tracers(iq)%iso_iGroup == iH2O .AND. niso > 0) THEN                          !=== WATER ISOTOPES
[4143]154      iName    = tracers(iq)%iso_iName
[4120]155      iPhase   = tracers(iq)%iso_iPhase
156      iqParent = tracers(iq)%iqParent
157      IF(tracers(iq)%iso_iZone == 0) THEN
[4984]158         if (tnat1) then
159                 tnat=1.0
160                 alpha_ideal=1.0
161                 write(*,*) 'attention dans dynetat0: les alpha_ideal sont a 1'
162         else
163          IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
[4325]164            CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
[4984]165         endif
[4120]166         CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized with a simplified Rayleigh distillation law.', modname)
[4325]167         q(:,:,:,iq) = q(:,:,:,iqParent)*tnat*(q(:,:,:,iqParent)/30.e-3)**(alpha_ideal-1.)
[4120]168      ELSE
169         CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized to its parent isotope concentration.', modname)
[4492]170         ! Camille 9 mars 2023: attention!! seuls les tags qui correspondent à
171         ! izone=izone_init (définie dans isotrac_mod) sont initialisés comme
172         ! les parents. Sinon, c'est nul.
173         ! j'ai fait ça en attendant, mais il faudrait initialiser proprement en
174         ! remplacant 1 par izone_init dans la ligne qui suit.
175         IF(tracers(iq)%iso_iZone == 1) THEN
176           q(:,:,:,iq) = q(:,:,:,iqIsoPha(iName,iPhase))
177         ELSE
178           q(:,:,:,iq) = 0.
179         END IF
[4120]180      END IF
181    !--------------------------------------------------------------------------------------------------------------------------
182    ELSE                                                                                 !=== MISSING: SET TO 0
183      CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized to zero', modname)
184      q(:,:,:,iq)=0.
185    !--------------------------------------------------------------------------------------------------------------------------
[2293]186    END IF
187  END DO
188
189  CALL err(NF90_CLOSE(fID),"close",fichnom)
190  day_ini=day_ini+INT(time)
191  time=time-INT(time)
192
193
194  CONTAINS
195
196
197SUBROUTINE check_dim(n1,n2,str1,str2)
198  INTEGER,          INTENT(IN) :: n1, n2
199  CHARACTER(LEN=*), INTENT(IN) :: str1, str2
[4046]200  CHARACTER(LEN=maxlen) :: s1, s2
[4193]201  IF(n1/=n2) CALL abort_gcm(TRIM(modname), 'value of "'//TRIM(str1)//'" = '//TRIM(int2str(n1))// &
202   ' read in starting file differs from gcm value of "'//TRIM(str2)//'" = '//TRIM(int2str(n2)), 1)
[2293]203END SUBROUTINE check_dim
204
205
206SUBROUTINE get_var1(var,v)
207  CHARACTER(LEN=*), INTENT(IN)  :: var
208  REAL,             INTENT(OUT) :: v(:)
209  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
210  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
211END SUBROUTINE get_var1
212
213
214SUBROUTINE get_var2(var,v)
215  CHARACTER(LEN=*), INTENT(IN)  :: var
216  REAL,             INTENT(OUT) :: v(:,:)
217  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
218  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
219END SUBROUTINE get_var2
220
221
222SUBROUTINE get_var3(var,v)
223  CHARACTER(LEN=*), INTENT(IN)  :: var
224  REAL,             INTENT(OUT) :: v(:,:,:)
225  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
226  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
227END SUBROUTINE get_var3
228
229
230SUBROUTINE err(ierr,typ,nam)
231  INTEGER,          INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
232  CHARACTER(LEN=*), INTENT(IN) :: typ    !--- TYPE OF OPERATION
233  CHARACTER(LEN=*), INTENT(IN) :: nam    !--- FIELD/FILE NAME
234  IF(ierr==NF90_NoERR) RETURN
235  SELECT CASE(typ)
[4120]236    CASE('inq');   mesg="Field <"//TRIM(nam)//"> is missing"
237    CASE('get');   mesg="Reading failed for <"//TRIM(nam)//">"
238    CASE('open');  mesg="File opening failed for <"//TRIM(nam)//">"
239    CASE('close'); mesg="File closing failed for <"//TRIM(nam)//">"
[2293]240  END SELECT
[4120]241  CALL ABORT_gcm(TRIM(modname),TRIM(mesg),ierr)
[2293]242END SUBROUTINE err
243
244END SUBROUTINE dynetat0
Note: See TracBrowser for help on using the repository browser.