Context Navigation

dynetat0_loc.F90 @ 5153

Last change on this file since 5153 was 5084, checked in by Laurent Fairhead, 4 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

Rev	Line
[2299]	1	SUBROUTINE dynetat0_loc(fichnom,vcov,ucov,teta,q,masse,ps,phis,time)
[1632]	2	!
[2299]	3	!-------------------------------------------------------------------------------
	4	! Authors: P. Le Van , L.Fairhead
	5	!-------------------------------------------------------------------------------
	6	! Purpose: Initial state reading.
	7	!-------------------------------------------------------------------------------
	8	USE parallel_lmdz
[4325]	9	USE infotrac, ONLY: nqtot, tracers, niso, iqIsoPha, iH2O, isoName
[4436]	10	USE strings_mod, ONLY: maxlen, msg, strStack, real2str, int2str, strIdx
[5084]	11	USE netcdf, ONLY: NF90_OPEN, NF90_NOWRITE, NF90_INQUIRE_DIMENSION, NF90_INQ_VARID, &
[4120]	12	NF90_CLOSE, NF90_GET_VAR, NF90_INQUIRE_VARIABLE, NF90_NoErr
[4325]	13	USE readTracFiles_mod, ONLY: new2oldH2O, newHNO3, oldHNO3, getKey
[2299]	14	USE control_mod, ONLY: planet_type
	15	USE assert_eq_m, ONLY: assert_eq
[2600]	16	USE comvert_mod, ONLY: pa,preff
[4120]	17	USE comconst_mod, ONLY: cpp, daysec, dtvr, g, im, jm, kappa, lllm, omeg, rad
[2603]	18	USE logic_mod, ONLY: fxyhypb, ysinus
[2598]	19	USE serre_mod, ONLY: clon, clat, grossismx, grossismy
[4120]	20	USE temps_mod, ONLY: annee_ref, day_ini, day_ref, itau_dyn, start_time
[2622]	21	USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
[4120]	22
[2299]	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:
[4120]	41	CHARACTER(LEN=maxlen) :: mesg, var, modname, oldVar
[2299]	42	INTEGER, PARAMETER :: length=100
[4301]	43	INTEGER :: iq, fID, vID, idecal, ierr, iqParent, iName, iZone, iPhase, ix
[4984]	44	REAL :: time,tab_cntrl(length) !--- RUN PARAMS TABLE
	45	REAL :: tnat, alpha_ideal
[2299]	46	REAL, ALLOCATABLE :: vcov_glo(:,:),masse_glo(:,:), ps_glo(:)
	47	REAL, ALLOCATABLE :: ucov_glo(:,:), q_glo(:,:), phis_glo(:)
	48	REAL, ALLOCATABLE :: teta_glo(:,:)
[4301]	49	LOGICAL :: lSkip, ll
[4984]	50	LOGICAL,PARAMETER :: tnat1=.TRUE.
[2299]	51	!-------------------------------------------------------------------------------
	52	modname="dynetat0_loc"
[1632]	53
[2299]	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)
[1632]	58
[2299]	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
[4120]	62	CALL msg('NOTE NOTE NOTE : Planeto-like start files', modname)
[2299]	63	idecal = 4
	64	annee_ref = 2000
	65	ELSE
[4120]	66	CALL msg('NOTE NOTE NOTE : Earth-like start files', modname)
[2299]	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
[1632]	104
[2299]	105	day_ini = tab_cntrl(30)
	106	itau_dyn = tab_cntrl(31)
[2615]	107	start_time = tab_cntrl(32)
[1632]	108
[2299]	109	!-------------------------------------------------------------------------------
[4120]	110	CALL msg('rad, omeg, g, cpp, kappa = '//TRIM(strStack(real2str([rad,omeg,g,cpp,kappa]))), modname)
[2299]	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)
[1632]	121
[2299]	122	var="temps"
	123	IF(NF90_INQ_VARID(fID,var,vID)/=NF90_NoErr) THEN
[4120]	124	CALL msg('missing field <temps> ; trying with <Time>', modname)
	125	var="Time"
[2299]	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)
[1632]	129
[2299]	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)
[1632]	133
[2299]	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)
[1673]	137
[2299]	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)
[1632]	141
[2299]	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)
[1632]	145
[2299]	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)
[1632]	153
[2299]	154	!--- Tracers
	155	ALLOCATE(q_glo(ip1jmp1,llm))
[4301]	156	ll = .FALSE.
[4263]	157	#ifdef REPROBUS
[4301]	158	ll = NF90_INQ_VARID(fID, 'HNO3tot', vID) /= NF90_NoErr !--- DETECT OLD REPRO start.nc FILE
[4263]	159	#endif
[2299]	160	DO iq=1,nqtot
[4120]	161	var = tracers(iq)%name
[4301]	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
[4265]	170	!--------------------------------------------------------------------------------------------------------------------------
[4301]	171	IF(NF90_INQ_VARID(fID, var, vID) == NF90_NoErr .AND. .NOT.lSkip) THEN !=== REGULAR CASE: AVAILABLE VARIABLE
[4265]	172	CALL get_var2(var,q_glo); q(ijb_u:ije_u,:,iq)=q_glo(ijb_u:ije_u,:)
	173	!--------------------------------------------------------------------------------------------------------------------------
[4301]	174	ELSE IF(NF90_INQ_VARID(fID, oldVar, vID) == NF90_NoErr) THEN !=== TRY WITH ALTERNATE NAME
[4265]	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	!--------------------------------------------------------------------------------------------------------------------------
[4120]	178	ELSE IF(tracers(iq)%iso_iGroup == iH2O .AND. niso > 0) THEN !=== WATER ISOTOPES
[4143]	179	iName = tracers(iq)%iso_iName
[4120]	180	iPhase = tracers(iq)%iso_iPhase
	181	iqParent = tracers(iq)%iqParent
	182	IF(tracers(iq)%iso_iZone == 0) THEN
[4984]	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))) &
[4325]	189	CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
[4984]	190	endif
[4120]	191	CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized with a simplified Rayleigh distillation law.', modname)
[4325]	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.)
[4490]	193	! Camille 9 mars 2023: point de vigilence: initialisation incohérente
	194	! avec celle de xt_ancien dans la physiq.
[4120]	195	ELSE
	196	CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized to its parent isotope concentration.', modname)
[4490]	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
[4984]	203	q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqIsoPha(iName,iPhase))
[4490]	204	ELSE
	205	q(ijb_u:ije_u,:,iq) = 0.
	206	ENDIF
[4120]	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	!--------------------------------------------------------------------------------------------------------------------------
[2299]	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)
[1673]	219
	220
[2299]	221	CONTAINS
[1632]	222
	223
[2299]	224	SUBROUTINE check_dim(n1,n2,str1,str2)
	225	INTEGER, INTENT(IN) :: n1, n2
	226	CHARACTER(LEN=*), INTENT(IN) :: str1, str2
[4046]	227	CHARACTER(LEN=maxlen) :: s1, s2
[4301]	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)
[2299]	230	END SUBROUTINE check_dim
[1632]	231
	232
[2299]	233	SUBROUTINE 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
[1632]	238
[2299]	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
	255	END SUBROUTINE get_var1
[1632]	256
[2299]	257	SUBROUTINE get_var2(var,v)
	258	CHARACTER(LEN=*), INTENT(IN) :: var
	259	REAL, INTENT(OUT) :: v(:,:)
[3967]	260	REAL, ALLOCATABLE :: w4(:,:,:,:), w3(:,:,:)
[2299]	261	INTEGER :: nn(4), dids(4), k, nd
[3967]	262
	263
[2299]	264	CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
[3967]	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
[2299]	272	DO k=1,nd; ierr=NF90_INQUIRE_DIMENSION(fID,dids(k),len=nn(k)); END DO
[3967]	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
[2299]	282	END SUBROUTINE get_var2
[1632]	283
	284
[2299]	285	SUBROUTINE 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)
[4120]	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)//">"
[2299]	295	END SELECT
[4120]	296	CALL ABORT_gcm(TRIM(modname),TRIM(mesg),ierr)
[2299]	297	END SUBROUTINE err
[1632]	298
[2299]	299	END SUBROUTINE dynetat0_loc

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format