source: LMDZ6/trunk/libf/dyn3dmem/dynetat0_loc.F90 @ 5087

Last change on this file since 5087 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
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 12.3 KB
Line 
1SUBROUTINE dynetat0_loc(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!-------------------------------------------------------------------------------
8  USE parallel_lmdz
9  USE infotrac,    ONLY: nqtot, tracers, niso, iqIsoPha, iH2O, isoName
10  USE strings_mod, ONLY: maxlen, msg, strStack, real2str, int2str, strIdx
11  USE netcdf,      ONLY: NF90_OPEN,  NF90_NOWRITE, NF90_INQUIRE_DIMENSION, NF90_INQ_VARID, &
12                         NF90_CLOSE, NF90_GET_VAR, NF90_INQUIRE_VARIABLE,  NF90_NoErr
13  USE readTracFiles_mod, ONLY: new2oldH2O, newHNO3, oldHNO3, getKey
14  USE control_mod, ONLY: planet_type
15  USE assert_eq_m, ONLY: assert_eq
16  USE comvert_mod, ONLY: pa,preff
17  USE comconst_mod, ONLY: cpp, daysec, dtvr, g, im, jm, kappa, lllm, omeg, rad
18  USE logic_mod, ONLY: fxyhypb, ysinus
19  USE serre_mod, ONLY: clon, clat, grossismx, grossismy
20  USE temps_mod, ONLY: annee_ref, day_ini, day_ref, itau_dyn, start_time
21  USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
22
23  IMPLICIT NONE
24  include "dimensions.h"
25  include "paramet.h"
26  include "comgeom.h"
27  include "description.h"
28  include "iniprint.h"
29!===============================================================================
30! Arguments:
31  CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
32  REAL, INTENT(OUT) ::  vcov(ijb_v:ije_v,llm)      !--- V COVARIANT WIND
33  REAL, INTENT(OUT) ::  ucov(ijb_u:ije_u,llm)      !--- U COVARIANT WIND
34  REAL, INTENT(OUT) ::  teta(ijb_u:ije_u,llm)      !--- POTENTIAL TEMP.
35  REAL, INTENT(OUT) ::     q(ijb_u:ije_u,llm,nqtot)!--- TRACERS
36  REAL, INTENT(OUT) :: masse(ijb_u:ije_u,llm)      !--- MASS PER CELL
37  REAL, INTENT(OUT) ::    ps(ijb_u:ije_u)          !--- GROUND PRESSURE
38  REAL, INTENT(OUT) ::  phis(ijb_u:ije_u)          !--- GEOPOTENTIAL
39!===============================================================================
40! Local variables:
41  CHARACTER(LEN=maxlen) :: mesg, var, modname, oldVar
42  INTEGER, PARAMETER :: length=100
43  INTEGER :: iq, fID, vID, idecal, ierr, iqParent, iName, iZone, iPhase, ix
44  REAL    :: time,tab_cntrl(length)    !--- RUN PARAMS TABLE
45  REAL    :: tnat, alpha_ideal
46  REAL,             ALLOCATABLE :: vcov_glo(:,:),masse_glo(:,:),   ps_glo(:)
47  REAL,             ALLOCATABLE :: ucov_glo(:,:),    q_glo(:,:), phis_glo(:)
48  REAL,             ALLOCATABLE :: teta_glo(:,:)
49  LOGICAL :: lSkip, ll
50  LOGICAL,PARAMETER :: tnat1=.TRUE.
51!-------------------------------------------------------------------------------
52  modname="dynetat0_loc"
53
54!--- Initial state file opening
55  var=fichnom
56  CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
57  CALL get_var1("controle",tab_cntrl)
58
59!!! AS: idecal is a hack to be able to read planeto starts...
60!!!     .... while keeping everything OK for LMDZ EARTH
61  IF(planet_type=="generic") THEN
62    CALL msg('NOTE NOTE NOTE : Planeto-like start files', modname)
63    idecal = 4
64    annee_ref  = 2000
65  ELSE
66    CALL msg('NOTE NOTE NOTE : Earth-like start files', modname)
67    idecal = 5
68    annee_ref  = tab_cntrl(5)
69  END IF
70  im         = tab_cntrl(1)
71  jm         = tab_cntrl(2)
72  lllm       = tab_cntrl(3)
73  day_ref    = tab_cntrl(4)
74  rad        = tab_cntrl(idecal+1)
75  omeg       = tab_cntrl(idecal+2)
76  g          = tab_cntrl(idecal+3)
77  cpp        = tab_cntrl(idecal+4)
78  kappa      = tab_cntrl(idecal+5)
79  daysec     = tab_cntrl(idecal+6)
80  dtvr       = tab_cntrl(idecal+7)
81  etot0      = tab_cntrl(idecal+8)
82  ptot0      = tab_cntrl(idecal+9)
83  ztot0      = tab_cntrl(idecal+10)
84  stot0      = tab_cntrl(idecal+11)
85  ang0       = tab_cntrl(idecal+12)
86  pa         = tab_cntrl(idecal+13)
87  preff      = tab_cntrl(idecal+14)
88!
89  clon       = tab_cntrl(idecal+15)
90  clat       = tab_cntrl(idecal+16)
91  grossismx  = tab_cntrl(idecal+17)
92  grossismy  = tab_cntrl(idecal+18)
93!
94  IF ( tab_cntrl(idecal+19)==1. )  THEN
95    fxyhypb  = .TRUE.
96!   dzoomx   = tab_cntrl(25)
97!   dzoomy   = tab_cntrl(26)
98!   taux     = tab_cntrl(28)
99!   tauy     = tab_cntrl(29)
100  ELSE
101    fxyhypb = .FALSE.
102    ysinus  = tab_cntrl(idecal+22)==1.
103  END IF
104
105  day_ini    = tab_cntrl(30)
106  itau_dyn   = tab_cntrl(31)
107  start_time = tab_cntrl(32)
108
109!-------------------------------------------------------------------------------
110  CALL msg('rad, omeg, g, cpp, kappa = '//TRIM(strStack(real2str([rad,omeg,g,cpp,kappa]))), modname)
111  CALL check_dim(im,iim,'im','im')
112  CALL check_dim(jm,jjm,'jm','jm')
113  CALL check_dim(lllm,llm,'lm','lllm')
114  CALL get_var1("rlonu",rlonu)
115  CALL get_var1("rlatu",rlatu)
116  CALL get_var1("rlonv",rlonv)
117  CALL get_var1("rlatv",rlatv)
118  CALL get_var1("cu"  ,cu)
119  CALL get_var1("cv"  ,cv)
120  CALL get_var1("aire",aire)
121
122  var="temps"
123  IF(NF90_INQ_VARID(fID,var,vID)/=NF90_NoErr) THEN
124    CALL msg('missing field <temps> ; trying with <Time>', modname)
125    var="Time"
126    CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
127  END IF
128  CALL err(NF90_GET_VAR(fID,vID,time),"get",var)
129
130  ALLOCATE(phis_glo(ip1jmp1))
131  CALL get_var1("phisinit",phis_glo)
132  phis (ijb_u:ije_u)  =phis_glo(ijb_u:ije_u);    DEALLOCATE(phis_glo)
133
134  ALLOCATE(ucov_glo(ip1jmp1,llm))
135  CALL get_var2("ucov",ucov_glo)
136  ucov (ijb_u:ije_u,:)=ucov_glo(ijb_u:ije_u,:);  DEALLOCATE(ucov_glo)
137
138  ALLOCATE(vcov_glo(ip1jm,llm))
139  CALL get_var2("vcov",vcov_glo)
140  vcov (ijb_v:ije_v,:)=vcov_glo(ijb_v:ije_v,:);  DEALLOCATE(vcov_glo)
141
142  ALLOCATE(teta_glo(ip1jmp1,llm))
143  CALL get_var2("teta",teta_glo)
144  teta (ijb_u:ije_u,:)=teta_glo(ijb_u:ije_u,:);  DEALLOCATE(teta_glo)
145
146  ALLOCATE(masse_glo(ip1jmp1,llm))
147  CALL get_var2("masse",masse_glo)
148  masse(ijb_u:ije_u,:)=masse_glo(ijb_u:ije_u,:); DEALLOCATE(masse_glo)
149 
150  ALLOCATE(ps_glo(ip1jmp1))
151  CALL get_var1("ps",ps_glo)
152  ps   (ijb_u:ije_u)  =   ps_glo(ijb_u:ije_u);   DEALLOCATE(ps_glo)
153
154!--- Tracers
155  ALLOCATE(q_glo(ip1jmp1,llm))
156  ll = .FALSE.
157#ifdef REPROBUS
158  ll = NF90_INQ_VARID(fID, 'HNO3tot', vID) /= NF90_NoErr                                 !--- DETECT OLD REPRO start.nc FILE
159#endif
160  DO iq=1,nqtot
161    var = tracers(iq)%name
162    oldVar = new2oldH2O(var)
163    lSkip = ll .AND. var == 'HNO3'                                                       !--- FORCE "HNO3_g" READING FOR "HNO3"
164#ifdef REPROBUS
165    ix = strIdx(newHNO3, var); IF(ix /= 0) oldVar = oldHNO3(ix)                          !--- REPROBUS HNO3 exceptions
166#endif
167#ifdef INCA
168    IF(var == 'O3') oldVar = 'OX'                                                        !--- DEAL WITH INCA OZONE EXCEPTION
169#endif
170    !--------------------------------------------------------------------------------------------------------------------------
171    IF(NF90_INQ_VARID(fID, var, vID) == NF90_NoErr .AND. .NOT.lSkip) THEN                !=== REGULAR CASE: AVAILABLE VARIABLE
172      CALL get_var2(var,q_glo); q(ijb_u:ije_u,:,iq)=q_glo(ijb_u:ije_u,:)
173    !--------------------------------------------------------------------------------------------------------------------------
174    ELSE IF(NF90_INQ_VARID(fID, oldVar, vID) == NF90_NoErr) THEN                         !=== TRY WITH ALTERNATE NAME
175      CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized to <'//TRIM(oldVar)//'>', modname)
176      CALL get_var2(oldVar, q_glo); q(ijb_u:ije_u,:,iq)=q_glo(ijb_u:ije_u,:)
177    !--------------------------------------------------------------------------------------------------------------------------
178    ELSE IF(tracers(iq)%iso_iGroup == iH2O .AND. niso > 0) THEN                          !=== WATER ISOTOPES
179      iName    = tracers(iq)%iso_iName
180      iPhase   = tracers(iq)%iso_iPhase
181      iqParent = tracers(iq)%iqParent
182      IF(tracers(iq)%iso_iZone == 0) THEN
183         if (tnat1) then
184                 tnat=1.0
185                 alpha_ideal=1.0
186                 write(*,*) 'attention dans dynetat0: les alpha_ideal sont a 1'
187         else
188          IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
189            CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
190         endif
191         CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized with a simplified Rayleigh distillation law.', modname)
192         q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqParent)*tnat*(q(ijb_u:ije_u,:,iqParent)/30.e-3)**(alpha_ideal-1.)
193         ! Camille 9 mars 2023: point de vigilence: initialisation incohérente
194         ! avec celle de xt_ancien dans la physiq.
195      ELSE
196         CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized to its parent isotope concentration.', modname)
197         ! Camille 9 mars 2023: attention!! seuls les tags qui correspondent à
198         ! izone=izone_init (définie dans isotrac_mod) sont initialisés comme
199         ! les parents. Sinon, c'est nul.
200         ! j'ai fait ça en attendant, mais il faudrait initialiser proprement en
201         ! remplacant 1 par izone_init dans la ligne qui suit.
202         IF(tracers(iq)%iso_iZone == 1) THEN
203           q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqIsoPha(iName,iPhase))
204         ELSE
205           q(ijb_u:ije_u,:,iq) =  0.
206         ENDIF
207      END IF
208    !--------------------------------------------------------------------------------------------------------------------------
209    ELSE                                                                                 !=== MISSING: SET TO 0
210      CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized to zero', modname)
211      q(ijb_u:ije_u,:,iq)=0.
212    !--------------------------------------------------------------------------------------------------------------------------
213    END IF
214  END DO
215  DEALLOCATE(q_glo)
216  CALL err(NF90_CLOSE(fID),"close",fichnom)
217  day_ini=day_ini+INT(time)
218  time=time-INT(time)
219
220
221  CONTAINS
222
223
224SUBROUTINE check_dim(n1,n2,str1,str2)
225  INTEGER,          INTENT(IN) :: n1, n2
226  CHARACTER(LEN=*), INTENT(IN) :: str1, str2
227  CHARACTER(LEN=maxlen) :: s1, s2
228  IF(n1/=n2) CALL abort_gcm(TRIM(modname), 'value of "'//TRIM(str1)//'" = '//TRIM(int2str(n1))// &
229   ' read in starting file differs from gcm value of "'//TRIM(str2)//'" = '//TRIM(int2str(n2)), 1)
230END SUBROUTINE check_dim
231
232
233SUBROUTINE get_var1(var,v)
234  CHARACTER(LEN=*), INTENT(IN)  :: var
235  REAL,             INTENT(OUT) :: v(:)
236  REAL,             ALLOCATABLE :: w2(:,:), w3(:,:,:)
237  INTEGER :: nn(3), dids(3), k, nd, ntot
238
239  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
240  ierr=NF90_INQUIRE_VARIABLE(fID,vID,ndims=nd)
241  IF(nd==1) THEN
242    CALL err(NF90_GET_VAR(fID,vID,v),"get",var); RETURN
243  END IF
244  ierr=NF90_INQUIRE_VARIABLE(fID,vID,dimids=dids)
245  DO k=1,nd; ierr=NF90_INQUIRE_DIMENSION(fID,dids(k),len=nn(k)); END DO
246  ntot=PRODUCT(nn(1:nd))
247  SELECT CASE(nd)
248    CASE(2); ALLOCATE(w2(nn(1),nn(2)))
249      CALL err(NF90_GET_VAR(fID,vID,w2),"get",var)
250      v=RESHAPE(w2,[ntot]); DEALLOCATE(w2)
251    CASE(3); ALLOCATE(w3(nn(1),nn(2),nn(3)))
252      CALL err(NF90_GET_VAR(fID,vID,w3),"get",var)
253      v=RESHAPE(w3,[ntot]); DEALLOCATE(w3)
254  END SELECT
255END SUBROUTINE get_var1
256
257SUBROUTINE get_var2(var,v)
258  CHARACTER(LEN=*), INTENT(IN)  :: var
259  REAL,             INTENT(OUT) :: v(:,:)
260  REAL,             ALLOCATABLE :: w4(:,:,:,:), w3(:,:,:)
261  INTEGER :: nn(4), dids(4), k, nd
262
263
264  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
265  ierr=NF90_INQUIRE_VARIABLE(fID,vID,ndims=nd)
266
267  IF(nd==1) THEN
268    CALL err(NF90_GET_VAR(fID,vID,v),"get",var); RETURN
269  END IF
270  ierr=NF90_INQUIRE_VARIABLE(fID,vID,dimids=dids)
271
272  DO k=1,nd; ierr=NF90_INQUIRE_DIMENSION(fID,dids(k),len=nn(k)); END DO
273
274  SELECT CASE(nd)
275  CASE(3); ALLOCATE(w3(nn(1),nn(2),nn(3)))
276     CALL err(NF90_GET_VAR(fID,vID,w3),"get",var)
277     v=RESHAPE(w3,[nn(1)*nn(2),nn(3)]); DEALLOCATE(w3)
278  CASE(4);  ALLOCATE(w4(nn(1),nn(2),nn(3),nn(4)))
279     CALL err(NF90_GET_VAR(fID,vID,w4),"get",var)
280     v=RESHAPE(w4,[nn(1)*nn(2),nn(3)]); DEALLOCATE(w4)
281  END SELECT
282END SUBROUTINE get_var2
283
284
285SUBROUTINE err(ierr,typ,nam)
286  INTEGER,          INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
287  CHARACTER(LEN=*), INTENT(IN) :: typ    !--- TYPE OF OPERATION
288  CHARACTER(LEN=*), INTENT(IN) :: nam    !--- FIELD/FILE NAME
289  IF(ierr==NF90_NoERR) RETURN
290  SELECT CASE(typ)
291    CASE('inq');   mesg="Field <"//TRIM(nam)//"> is missing"
292    CASE('get');   mesg="Reading failed for <"//TRIM(nam)//">"
293    CASE('open');  mesg="File opening failed for <"//TRIM(nam)//">"
294    CASE('close'); mesg="File closing failed for <"//TRIM(nam)//">"
295  END SELECT
296  CALL ABORT_gcm(TRIM(modname),TRIM(mesg),ierr)
297END SUBROUTINE err
298
299END SUBROUTINE dynetat0_loc
Note: See TracBrowser for help on using the repository browser.